Сопоставьте «joint_tests» со списком после установки модели «gls»

Я пытаюсь получить таблицу 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

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

@RonakShah Спасибо за ответ. Я добавил пример.

hnguyen 14.12.2020 05:11

это две разные модели, которые вы примеряете, причем gls что-то сложнее, чем lm. какой из них вы пытаетесь соответствовать? Я могу предвидеть, что функция varIdents() внутри gls создаст большую проблему. не могли бы вы быть немного более последовательным в примере?

StupidWolf 14.12.2020 05:11

@StupidWolf Спасибо за ответ. Я хочу gls в конце концов, но я использовал lm для простоты. Насколько я понимаю, joint_tests не важно, что мы поставляем, если это поддерживается.

hnguyen 14.12.2020 05:13

это работает, но это не предназначено для совместимости с муррр и т. д. Что происходит, так это то, что ему нужно снова заглянуть в вашу среду для фрейма данных, а внутри муррр или lapply и т. д., это становится проблематичным. Насколько я понял, вам нужно указать дисперсию с помощью gls?

StupidWolf 14.12.2020 05:14

@StupidWolf Мне нужно указать структуру корреляции. Здесь я предполагаю, что каждому color разрешено иметь разную дисперсию (вместо объединенной дисперсии для всего набора данных).

hnguyen 14.12.2020 05:20

Вы согласны с использованием цикла for для получения результатов? Потребуется некоторое время, чтобы разобрать функциюjoint_test, чтобы заставить ее работать внутри муррр.

StupidWolf 14.12.2020 05:27

Я думаю, что это проблема масштаба. Если вместо отображения наjoint_tests написать функцию, соответствующую модели и вызывающуюjoint_tests, то я думаю, что это сработает. Проблема в том, что emmeans() (вызываемая методом join_tests()) нуждается в реконструкции данных и не может найти их в рамках вызова.

Russ Lenth 14.12.2020 19:02

Или, может быть, вы можете просто использовать map(joint_tests, data = .)? Это должно работать, если данные передаются в join_tests.

Russ Lenth 14.12.2020 19:05

@RussLenth Спасибо за ваши предложения. Я новичок в циклах и настройке функций, поэтому я еще не ушел далеко. Второе предложение не сработало. Я получил это сообщение об ошибке Error in as_mapper(.f, ...) : argument ".f" is missing, with no default

hnguyen 14.12.2020 19:16

ОК, я попробовал и получил ту же ошибку, что и вы; но это может сработать, если вместо этого поставить data = diamonds. Все же лучше написать свою собственную функцию для карты.

Russ Lenth 14.12.2020 23:39
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
0
10
130
2
Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

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

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)
}

Похоже, вам нужно вернуть «фиты» в конце

Russ Lenth 14.12.2020 22:40
Ответ принят как подходящий

По причинам, которые я буду исследовать, 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*

Другие вопросы по теме