Получение оценок средней величины эффекта и доверительных интервалов при отсутствии комбинации уровней факторов
Я подогнал модель с помощью пакета метафор с модераторским взаимодействием, для которой нет данных для комбинации уровней факторов.
library(orchaRd)
data(fish)
warm_dat <- fish
warm_dat$factor2<-as.factor(seq(c("A","B","C"),133) )
mod_fit <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
random = list(~1 | group_ID, ~1 | es_ID),
mods = ~ trait.type*factor2,
method = "REML", test = "t",
control=list(optimizer = "optim", optmethod = "Nelder-Mead"), data = warm_dat)
Я хотел использовать пакет orchaRd для получения средних размеров эффекта и доверительных интервалов, но получил ошибку.
results1 <- mod_results(mod_fit, group = "group_ID", mod = "trait.type", by = "Measurement_cat")
results2 <- mod_results(mod_fit, group = "group_ID", mod = "trait.type")
Ошибка в ref_grid(result,...): Что-то пошло не так: Несогласные элементы в опорной сетке.
Я просмотрел код, используемый в mod_results, и обнаружил, что он основан на emmeans. Я попытался адаптировать код, вручную добавив недостающую комбинацию уровней факторов в эталонную сетку, но полученные результаты оказались идентичны оценкам, полученным с помощью сводной функции.
coef1<-mod_fit$b
coef1[12]<-NA
names(coef1)<-c(row.names(mod_fit$b),"trait.typelife-history:factor2C")
vcov1<-stats::vcov(mod_fit)
vcov1<- rbind(vcov1,rep(NA, times=11))
vcov1<- cbind(vcov1,rep(NA, times=12))
row.names(vcov1)<-c(row.names(vcov(mod_fit)),"trait.typelife-history:factor2C")
colnames(vcov1)<-row.names(vcov1)
grid <- emmeans::qdrg(formula = stats::formula(mod_fit),
data = warm_dat, coef = coef1, vcov = vcov1,
df = mod_fit$k - 1)
mm <- emmeans::emmeans(grid, specs =~trait.type * factor2, df = as.numeric(mod_fit$ddf[[1]]))
mm_pi <- pred_interval_esmeans(mod_fit, mm, mod = trait.type)
Любые предложения относительно того, как получить оценки и доверительные интервалы в этом случае, будут высоко оценены.
Точнее, некоторые разработчики пакетов просто выбрасывают некоторые столбцы матрицы X; но для правильной обработки невозможности оценки вы должны сохранить эти столбцы и признать, что «выбрасывание» столбца на самом деле является способом ограничить его коэффициент регрессии равным нулю. Если вы посмотрите на имена столбцов model_matrix(formula(mod_fit), data=warm_dat) и сравните их с именами коэффициентов вашей модели, возможно, вам удастся собрать новый полный набор коэффициентов и установить для исключенных коэффициентов значение NA. Тогда это может сработать.
Да, я согласен, что в долгосрочной перспективе было бы лучше, если бы glmmTMB заполнили NA значения отсутствующих/неидентифицируемых параметров. Не то, чтобы мы дошли до этого - можно было бы открыть проблему...
@BenBolker Не проблема, теперь, когда у нас есть работающая поддержка emmeans для моделей glmmTMB (а также merMod). Он использует код, аналогичный тому, который я предлагаю здесь, для заполнения недостающих значений там, где они необходимы. Так что не планирую создавать проблему...
Ответы 1
После моего второго комментария я думаю, что вы можете сделать все наоборот, поскольку вы устанавливаете некоторые необходимые вам оценки NA, а не создаете более длинный список коэффициентов, включающий NA. Попробуйте что-то вроде этого:
X = model.matrix(stats::formula(mod_fit), data = warm_dat)
b = mod_fit$b
full.coef = X[1, ] * NA # a bunch of NAs with the right names
full.coef[names(b)] = b # so now full.coef has NAs only where excluded from b
grid <- emmeans::qdrg(formula = stats::formula(mod_fit),
data = warm_dat, coef = full.coef, vcov = vcov(mod_fit),
df = mod_fit$k - 1)
# etc.
Обратите внимание, что часть vcov должна быть неполной версией с удаленными строками и столбцами.
Обновлять
Код в ОП работает не совсем правильно. Итак, для factor2 я использовал
warm_dat$factor2 = factor(rep(LETTERS[1:3], 137))[1:410]
Предложенное мной исправление почти работает. Но дополнительная проблема заключается в том, что перехват имеет имя intrcpt вместо обычного (Intercept), возвращаемого model.matrix. (Что такого плохого в том, что разработчики пакетов следуют установленным соглашениям???) В любом случае это требует дополнительной строки кода для установки перехвата.
Ниже приведен результат, который я получаю. Обратите внимание, что второй элемент правильно помечен как неоценимый.
data(fish, package = "orchaRd")
warm_dat <- fish
warm_dat$factor2 = factor(rep(LETTERS[1:3], 137))[1:410]
mod_fit <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
random = list(~1 | group_ID, ~1 | es_ID),
mods = ~ trait.type*factor2,
method = "REML", test = "t",
control=list(optimizer = "optim", optmethod = "Nelder-Mead"), data = warm_dat)
## Warning: Redundant predictors dropped from the model.
library(emmeans)
b = mod_fit$b
full.coef = NA * model.matrix(stats::formula(mod_fit), data = warm_dat)[1, ]
full.coef[rownames(b)] = b # so now full.coef has NAs only where excluded from b
full.coef[1] = b[1] # needed because of unconventional naming of the intercept
RG <- emmeans::qdrg(formula = stats::formula(mod_fit),
data = warm_dat, coef = full.coef, vcov = vcov(mod_fit),
df = mod_fit$k - 1)
confint(RG)
## trait.type factor2 prediction SE df lower.CL upper.CL
## behaviour A 0.953537 0.2992 409 0.3654 1.5417
## life-history A nonEst NA NA NA NA
## morphology A 0.001582 0.0224 409 -0.0425 0.0456
## physiology A 0.038965 0.0454 409 -0.0502 0.1282
## behaviour B -0.028576 0.1827 409 -0.3878 0.3306
## life-history B 0.444830 0.1296 409 0.1901 0.6995
## morphology B 0.004905 0.0226 409 -0.0394 0.0492
## physiology B -0.000301 0.0467 409 -0.0922 0.0916
## behaviour C 0.057822 0.1889 409 -0.3135 0.4291
## life-history C 0.972144 0.1695 409 0.6389 1.3054
## morphology C -0.004610 0.0234 409 -0.0507 0.0415
## physiology C 0.004195 0.0404 409 -0.0752 0.0836
##
## Confidence level used: 0.95
Created on 2024-07-26 with reprex v2.1.0
Спасибо, что нашли время предложить решение. Это сработает, если я поменяю names(b) на rownames(b). Однако я не могу получить краткое описание объекта emmeans: mm <- emmeans::emmeans(grid, specs =~trait.type * factor2, df = as.numeric(mod_fit$ddf[[1]])) summary(mm)
@cosalofa, этот fish набор данных доступен в пакете R или где-то общедоступен? Было бы полезно иметь воспроизводимый пример.
Я вижу наборы данных fish в некоторых установленных мной пакетах, но ни один из них не имеет имен переменных, соответствующих тем, которые используются в примере: pp <- help.search("^fish$", agrep = FALSE, ignore.case = FALSE)$matches$Package; for (p in pp) { cat(p, "\n"); data("fish", package = p); print(names(fish)); rm("fish") }
Извиняюсь, забыл добавить пакет, в котором он находится (orchaRd). Я изменил исходное сообщение, чтобы включить это.
Я все еще не могу воспроизвести анализ, потому что строка warm_dat$factor2<-as.factor(seq(c("A","B","C"),133) ) выдает ошибку — недопустимое использование seq(). И я мог бы догадаться rep(c("A","B","C"),133)", за исключением того, что warm_dat имеет 410 строк и 3*133 = 399. Итак, не могли бы вы дать нам фактический код, который вы использовали, и сначала протестировать его?
Извините, я скопировал код из R в первый раз, и он сработал. Должно быть, я заранее обработал набор данных. Я вижу, что вы обновили свое решение. Теперь все работает отлично, большое спасибо!
@cosalofa Я также добавил ссылку на эту публикацию в более ранний выпуск, который я разместил на сайте метафора, и, возможно, разработчик соответствующим образом улучшит эту функциональность.
Если у вас есть недостающие ячейки, это создает проблему невозможности оценки, и вполне возможно, что эта ситуация не обрабатывается правильно пакетом метафора. Возможно, будет полезно попытаться сообщить об этой проблеме на веб-сайт этого пакета. Кроме того, вы можете получить некоторое представление, попробовав другой пакет (возможно, glmmTMB или nlme?), который может соответствовать чему-то достаточно близкому к тому, что вам нужно.