Блог

Я пытаюсь получить таблицу ANOVA типа 3 emmeans::joint_tests() из списка со следующим кодом. Я не совсем понимаю сообщение об ошибке. Код, который меня обучает, взят из http://pages.stat.wisc.edu /~yandell/R_for_data_sciences/curate/tidyverse.html

 library(dplyr)
library(nlme)
library(emmeans)
data("diamonds")
diamonds %>%
  split(.$cut) %>%
  map(~ gls(price ~ x   y   z,  
  weights = varIdent(form = ~ 1|color),
  data = .))%>%
map(summary)
 

Сообщение об ошибке , по — видимому , предполагает , что я каким — то образом сохраняю свои отдельные модели , а затем применяю joint_tests их . Чего я не понимаю, так это почему рабочий процесс работает для summary , но не для joint_tests . Когда мы подгоняем отдельные модели, это summary(model) или joint_tests(model) это печатает сводную таблицу или таблицу ANOVA соответственно.

 data("diamonds")
diamonds %>%
  split(.$cut) %>%
  map(~ gls(price ~ x   y   z,  
  weights = varIdent(form = ~ 1|color),
  data = .))%>%
map(joint_tests)

Error in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : Perhaps a 'data' or 'params' argument is needed
 

Использование map(~ joint_tests) дало странный список, подобный этому

 $Fair
function (object, by = NULL, show0df = FALSE, cov.reduce = range, 
    ...) 
{
    if (!inherits(object, "emmGrid")) {
        args = .zap.args(object = object, cov.reduce = cov.reduce, 
            ..., omit = "submodel")
        object = do.call(ref_grid, args)
    }
    facs = setdiff(names(object@levels), by)
    do.test = function(these, facs, result, ...) {
        if ((k <- length(these)) > 0) {
            emm = emmeans(object, these, by = by, ...)
            tst = test(contrast(emm, interaction = "consec"), 
                joint = TRUE, status = TRUE)
            tst = cbind(ord = k, `model term` = paste(these, 
                collapse = ":"), tst)
            result = rbind(result, tst)
            last = max(match(these, facs))
        }
        else last = 0
        if (last < (n <- length(facs))) 
            for (i in last   seq_len(n - last)) result = do.test(c(these, 
                facs[i]), facs, result, ...)
        result
    }
    result = suppressMessages(do.test(character(0), facs, NULL, 
        ...))
    result = result[order(result[[1]]), -1, drop = FALSE]
    if (!show0df) 
        result = result[result$df1 > 0, , drop = FALSE]
    class(result) = c("summary_emm", "data.frame")
    attr(result, "estName") = "F.ratio"
    attr(result, "by.vars") = by
    if (any(result$note != "")) {
        msg = character(0)
        if (any(result$note %in% c(" d", " d e"))) 
            msg = .dep.msg
        if (any(result$note %in% c("   e", " d e"))) 
            msg = c(msg, .est.msg)
        attr(result, "mesg") = msg
    }
    else result$note = NULL
    result
}
<bytecode: 0x7ff68eb4a0a8>
<environment: namespace:emmeans>

$Good
function (object, by = NULL, show0df = FALSE, cov.reduce = range, 
    ...) 
{
 

Вот как я обошелся joint_tests без списка.

 diamond.gls <-  gls(price ~ x   y   z,  
  weights = varIdent(form = ~ 1|color),
  data = diamonds)
joint_tests(diamond.gls)
model term df1      df2  F.ratio p.value
 x            1 14311.72 4898.859 <.0001 
 y            1 12964.08   46.231 <.0001 
 z            1  8380.71   24.576 <.0001
 

Пожалуйста, посмотрите, как я могу это исправить. Спасибо.

По причинам, которые я буду исследовать, joint_tests() требуется data аргумент, когда это gls модель, по крайней мере, при вызове в теле функции. Чтобы преодолеть это, нам нужно написать функцию, которая соответствует модели и выполняется joint_tests() . Вот параллельная иллюстрация:

 mod_jt = function(dat) {
  mod = gls(breaks ~ tension, data = dat)
  joint_tests(mod, data = dat)
}

warpbreaks %>% split(.$wool) %>% map(mod_jt) 
 

… и мы получаем результаты:

 $A
 model term df1 df2 F.ratio p.value
 tension      2  24   7.288 0.0034 


$B
 model term df1 df2 F.ratio p.value
 tension      2  24   4.059 0.0303 
 

Я думаю, что код, который у вас был, будет работать с lm моделью, по крайней мере, с новейшей версией emmeans*

Мы можем настроить набор данных с примерами, которые работают:

 dat = droplevels(subset(diamonds,cut %in% c("Ideal","Premium","Good")))
dat$X = cut(dat$z,c(-0.5,4,11))
dat$clarity = factor(dat$clarity,ordered=FALSE)
dat$color = factor(dat$color,ordered=FALSE)
fit = gls(price ~ X*clarity, weights = varIdent(form = ~ 1|color),
data=subset(dat,cut=="Ideal"))
joint_tests(fit)

 model term df1      df2   F.ratio p.value
 X            1 15002.85 12145.835 <.0001 
 clarity      7 11834.99   351.899 <.0001 
 X:clarity    7 11834.99   352.344 <.0001 
 

Чтобы это работало нормально для подмножества, нам нужно заставить его работать. Причина, по которой вы сталкиваетесь с ошибкой, заключается в том, что joint_tests() в вашей среде снова просматривается data.frame, а внутри функции map () это невозможно.

Один из простых способов — использовать цикл for и сохранить результаты в списке:

 fits = list()

for(i in unique(dat$cut)){

f  = gls(price ~ X*clarity,  
                weights = varIdent(form = ~ 1|color),
                data = subset(dat,cut==i))
res = joint_tests(f)
fits[[i]] = list(f=f,res=res)
}