Как получить коэффициенты и их доверительные интервалы в моделях со смешанными эффектами?
В lm
а также glm
модели, я использую функции coef
а также confint
для достижения цели:
m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
coef(m)
confint(m)
Теперь я добавил случайный эффект к модели - использовал модели смешанных эффектов, используя lmer
функция из пакета lme4. Но тогда функции coef
а также confint
не работай больше на меня!
> mix1 = lmer(resp ~ 0 + var1 + var1:var2 + (1|var3))
# var1, var3 categorical, var2 continuous
> coef(mix1)
Error in coef(mix1) : unable to align random and fixed effects
> confint(mix1)
Error: $ operator not defined for this S4 class
Я пытался Google и использовать документы, но безрезультатно. Пожалуйста, укажите мне в правильном направлении.
РЕДАКТИРОВАТЬ: Я также подумал, подходит ли этот вопрос больше к https://stats.stackexchange.com/ но я считаю, что это скорее технический, чем статистический, поэтому я пришел к выводу, что он подходит здесь (SO)... что вы думаете?
7 ответов
Не уверен, когда он был добавлен, но теперь функция confint() реализована в lme4. Например, работает следующий пример:
library(lme4)
m = lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
confint(m)
Существуют два новых пакета, lmerTest и lsmeans, которые могут рассчитывать 95% доверительные интервалы для lmer
а также glmer
выход. Может быть, вы можете посмотреть на них? И coefplot2, я думаю, может сделать это тоже (хотя, как указывает Бен ниже, не так изощренно, из стандартных ошибок в статистике Вальда, в отличие от приближения Кенварда-Роджера и / или Satterthwaite df, используемого в lmerTest
а также lsmeans
)... просто обидно, что в упаковке до сих пор нет встроенных заговорщиков lsmeans
(как есть в упаковке effects()
что, кстати, также возвращает 95% доверительные пределы lmer
а также glmer
объекты, но делает это путем переоснащения модели без каких-либо случайных факторов, что, очевидно, не правильно).
Я собираюсь добавить сюда немного. Если
m
это подогнанный
(g)lmer
модель (большинство из них тоже работают):
-
fixef(m)
- канонический способ извлечения коэффициентов из смешанных моделей (это соглашение началось сnlme
и перешел наlme4
) - вы можете получить полную таблицу коэффициентов с
coef(summary(m))
; если вы загрузили до подгонки модели, или преобразовали модель после подгонки (а затем загрузки) черезcoef(summary(as(m,"merModLmerTest")))
, то таблица коэффициентов будет включать p-значения. (Таблица коэффициентов представляет собой матрицу; вы можете извлечь столбцы, например,ctab[,"Estimate"]
,ctab[,"Pr(>|t|)"]
, или преобразовать матрицу в фрейм данных и использовать$
-индексирование.) - Как указано выше, вы можете получить доверительные интервалы профиля правдоподобия с помощью
confint(m)
; они могут потребовать больших вычислительных ресурсов. Если вы используетеconfint(m, method="Wald")
вы получите стандартные доверительные интервалы +/- 1.96SE. (lme
используетintervals(m)
вместоconfint()
.)
Если вы предпочитаете использовать
broom.mixed
:
-
tidy(m,effects="fixed")
дает вам таблицу с оценками, стандартными ошибками и т. д. -
tidy(as(m,"merModLmerTest"), effects="fixed")
(или подгонка сlmerTest
в первую очередь) включает p-значения - добавление дает (Wald) CI
- добавление
conf.method="profile"
(вместе сconf.int=TRUE
) дает КЭ профиля правдоподобия
Вы также можете получить доверительные интервалы с помощью параметрической начальной загрузки (
method="boot"
), что значительно медленнее, но в некоторых случаях точнее.
Предполагая нормальную аппроксимацию для фиксированных эффектов (что было бы возможно и для ограничения), мы можем получить 95% доверительные интервалы
оценка + 1,96* стандартная ошибка.
Следующее не относится к компонентам дисперсии / случайным эффектам.
library("lme4")
mylm <- lmer(Reaction ~ Days + (Days|Subject), data =sleepstudy)
# standard error of coefficient
days_se <- sqrt(diag(vcov(mylm)))[2]
# estimated coefficient
days_coef <- fixef(mylm)[2]
upperCI <- days_coef + 1.96*days_se
lowerCI <- days_coef - 1.96*days_se
Я предлагаю вам использовать старый добрый lme (в пакете nlme). У него есть ограничение, и, если вам нужно ограничение контрастов, есть ряд вариантов (оцениваемых в моделях, контраст в контрастах, glht в мульткомпьютере).
Почему p-значения и ограничение отсутствуют в lmer: см. Http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html.
Чтобы найти коэффициент, вы можете просто использовать функцию итога lme4
m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
m_summary <- summary(m)
иметь все коэффициенты:
m_summary$coefficient
Если вы хотите доверительный интервал, умножьте стандартную ошибку на 1,96:
CI <- m_summary$coefficient[,"Std. Error"]*1.96
print(CI)
Я бы предложил
tab_model()
функция от
sjPlot
пакет в качестве альтернативы. Чистый и читаемый вывод, готовый к уценке. Ссылка здесь и примеры здесь.
Для более склонных к зрению
plot_model()
из того же пакета тоже может пригодиться.
Альтернативное решение - через
parameters
пакет с использованием
model_parameters()
функция.