Я пытаюсь получить таблицу 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
Пожалуйста, посмотрите, как я могу это исправить. Спасибо.
это две разные модели, которые вы примеряете, причем gls что-то сложнее, чем lm. какой из них вы пытаетесь соответствовать? Я могу предвидеть, что функция varIdents()
внутри gls создаст большую проблему. не могли бы вы быть немного более последовательным в примере?
@StupidWolf Спасибо за ответ. Я хочу gls
в конце концов, но я использовал lm
для простоты. Насколько я понимаю, joint_tests
не важно, что мы поставляем, если это поддерживается.
это работает, но это не предназначено для совместимости с муррр и т. д. Что происходит, так это то, что ему нужно снова заглянуть в вашу среду для фрейма данных, а внутри муррр или lapply и т. д., это становится проблематичным. Насколько я понял, вам нужно указать дисперсию с помощью gls?
@StupidWolf Мне нужно указать структуру корреляции. Здесь я предполагаю, что каждому color
разрешено иметь разную дисперсию (вместо объединенной дисперсии для всего набора данных).
Вы согласны с использованием цикла for для получения результатов? Потребуется некоторое время, чтобы разобрать функциюjoint_test, чтобы заставить ее работать внутри муррр.
Я думаю, что это проблема масштаба. Если вместо отображения наjoint_tests написать функцию, соответствующую модели и вызывающуюjoint_tests, то я думаю, что это сработает. Проблема в том, что emmeans() (вызываемая методом join_tests()) нуждается в реконструкции данных и не может найти их в рамках вызова.
Или, может быть, вы можете просто использовать map(joint_tests, data = .)
? Это должно работать, если данные передаются в join_tests.
@RussLenth Спасибо за ваши предложения. Я новичок в циклах и настройке функций, поэтому я еще не ушел далеко. Второе предложение не сработало. Я получил это сообщение об ошибке Error in as_mapper(.f, ...) : argument ".f" is missing, with no default
ОК, я попробовал и получил ту же ошибку, что и вы; но это может сработать, если вместо этого поставить data = diamonds
. Все же лучше написать свою собственную функцию для карты.
Мы можем настроить набор данных с примерами, которые работают:
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)
}
Похоже, вам нужно вернуть «фиты» в конце
По причинам, которые я буду исследовать, 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*
@RonakShah Спасибо за ответ. Я добавил пример.