Почему мне нужно жестко запрограммировать формулу в emmeans() вместо динамического извлечения формулы из динамической линейной модели?

Обновлено: добавлен не-Shiny MRE, хотя у него другая ошибка, чем у Shiny...

Рассмотрим следующий MRE (обратите внимание, что приложение не будет отображаться правильно, этот код предназначен только для запуска его для MRE, чтобы вы могли протестировать его в браузере().) Этот MRE имитирует часть более крупного приложения, указанную пользователем. модель, имеющая значение в ANOVA, и собирается провести соответствующие апостериорные исследования. В этом случае апостериорно при наличии гетероскедастичности, поэтому нам нужно построить модель с разными дисперсиями для каждой комбинации эффектов.

library(shiny)
require(emmeans)
require(nlme)
require(DT)

data<-structure(list(percent = c(2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 
                                 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 
                                 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 
                                 2, 2, 1, 1, 1, 1, 1, 1, 1, 1),
                     cycle = c(4, 4, 4, 4, 4, 4, 4, 
                               4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
                               3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 
                               1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                     density = c(21, 
                                 20, 21, 21, 22, 20, 20, 19, 18, 19, 19, 19, 19, 18, 18, 19, 20, 
                                 22, 22, 21, 22, 22, 20, 21, 20, 20, 19, 19, 19, 19, 19, 19, 23, 
                                 25, 24, 24, 26, 23, 24, 24, 23, 22, 22, 23, 21, 24, 22, 22, 16, 
                                 17, 16, 15, 17, 16, 16, 16, 26, 27, 26, 25, 25, 26, 26, 26), 
                     run = c(8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 4L, 4L, 4L, 4L, 4L, 
                             4L, 4L, 4L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 3L, 3L, 3L, 3L, 
                             3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 
                             2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 
                             1L, 1L, 1L, 1L, 1L, 1L)), class = "data.frame", row.names = c(NA, 
                                                                                           -64L))


# Define UI for application
ui <- fluidPage(

    # Application title
    titlePanel("Blah"),

    # Sidebar with a slider input for number of bins 
    sidebarLayout(
        sidebarPanel(
        ),

        # Show a plot of the generated distribution
        mainPanel(
           #plotOutput("distPlot"),
           dataTableOutput("otherthing"),
           textOutput("post_hoc")
        )
    )
)

# Define server logic required to draw a histogram
server <- function(input, output) {
  
  aov_model <- reactiveVal(NULL)
  
  output$otherthing<-renderDataTable({
    formula<-"density~percent*cycle"#selected by user
    aov_mod <-lm(formula = formula(formula),data = data)
    data$cycle<-factor(data$cycle)
    data$percent<-factor(data$percent)
    aov_model(aov_mod)
    #browser()
    print(aov_mod[["coefficients"]])
  })
  
  output$post_hoc<-renderText({
    data$group<-interaction(data[c(1,2)])
    gls_formula<-formula(aov_model())
    #environment(gls_formula) <- .GlobalEnv
    model<-gls(gls_formula, data = data, weights = varIdent(form = ~1 | group))
    model2<-gls(density~percent*cycle, data = data, weights = varIdent(form = ~1 | group))
    browser()
    results<-emmeans(model,
                     pairwise~percent*cycle,
                     adjust = "tukey",
                     type = "response",
                     data=data)
    
    results2<-emmeans(model2,
                     pairwise~percent*cycle,
                     adjust = "tukey",
                     type = "response",
                     data=data)
    print(results2)
  })

Обратите внимание при разрыве браузера, что результаты контрастов работают нормально, если вы вручную наберете «попарно ~ процент * цикл» для формулы gls, которая позже будет использоваться в emmeans() (выполнение результатов2). Но не будет, если вы используете извлеченную формулу из реактивной линейной модели (выполнение результатов) в emmeans(). В этом случае он выдаст: Ошибка в eval(call$model): объект «gls_formula» не найден и из-за этого еще одна ошибка.

В реальном приложении все это динамически выбирается пользователем и потенциально изменяется приложением, поэтому мне действительно нужна модель gls, использующая динамически извлекаемую формулу, чтобы она работала, а не просто перестраивала выбранные эффекты.

Я не понимаю, почему извлеченная формула, которая, насколько я могу судить, идентична жестко запрограммированной, не работает. Я даже пробовал установить среду на извлеченную, но это, похоже, ничего не меняет.

Ниже приведен пример без Shiny. Выдает другую ошибку, но, возможно, это связано с той же проблемой? Эта ошибка: Ошибка в вызове $model[[2]]: объект типа «символ» не является подмножеством.

require(emmeans)
require(nlme)

data<-structure(list(percent = c(2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 
                                 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 
                                 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 
                                 2, 2, 1, 1, 1, 1, 1, 1, 1, 1),
                     cycle = c(4, 4, 4, 4, 4, 4, 4, 
                               4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
                               3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 
                               1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                     density = c(21, 
                                 20, 21, 21, 22, 20, 20, 19, 18, 19, 19, 19, 19, 18, 18, 19, 20, 
                                 22, 22, 21, 22, 22, 20, 21, 20, 20, 19, 19, 19, 19, 19, 19, 23, 
                                 25, 24, 24, 26, 23, 24, 24, 23, 22, 22, 23, 21, 24, 22, 22, 16, 
                                 17, 16, 15, 17, 16, 16, 16, 26, 27, 26, 25, 25, 26, 26, 26), 
                     run = c(8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 4L, 4L, 4L, 4L, 4L, 
                             4L, 4L, 4L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 3L, 3L, 3L, 3L, 
                             3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 
                             2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 
                             1L, 1L, 1L, 1L, 1L, 1L)), class = "data.frame", row.names = c(NA, 
                                                                                           -64L))



formula<-"density~percent*cycle" #selected by user
effect<-"percent*cycle" #selected by user
data$cycle<-factor(data$cycle)
data$percent<-factor(data$percent)

aov_mod <-lm(formula = formula(formula),data = data)

#later in the app
data$group<-interaction(data[c(1,2)])
gls_formula<-formula(aov_mod)
#environment(gls_formula) <- .GlobalEnv
model<-gls(gls_formula, data = data, weights = varIdent(form = ~1 | group)) #using extracted formula
model2<-gls(density~percent*cycle, data = data, weights = varIdent(form = ~1 | group)) #using hard coded formula
emmeans(model,
        formula(paste("pairwise ~",gsub(pattern = ":",replacement = "*",x = effect))),
        adjust = "tukey",
        type = "response",
        data=data)
emmeans(model2,
        formula(paste("pairwise ~",gsub(pattern = ":",replacement = "*",x = effect))),
        adjust = "tukey",
        type = "response",
        data=data)