Как получить коэффициенты и их доверительные интервалы в моделях со смешанными эффектами?

Question

В 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)... что вы думаете?

Answer 1

Не уверен, когда он был добавлен, но теперь функция confint() реализована в lme4. Например, работает следующий пример:

library(lme4)
m = lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
confint(m)

Answer 2

Существуют два новых пакета, lmerTest и lsmeans, которые могут рассчитывать 95% доверительные интервалы для lmer а также glmer выход. Может быть, вы можете посмотреть на них? И coefplot2, я думаю, может сделать это тоже (хотя, как указывает Бен ниже, не так изощренно, из стандартных ошибок в статистике Вальда, в отличие от приближения Кенварда-Роджера и / или Satterthwaite df, используемого в lmerTest а также lsmeans)... просто обидно, что в упаковке до сих пор нет встроенных заговорщиков lsmeans (как есть в упаковке effects() что, кстати, также возвращает 95% доверительные пределы lmer а также glmer объекты, но делает это путем переоснащения модели без каких-либо случайных факторов, что, очевидно, не правильно).

Answer 3

Я собираюсь добавить сюда немного. Если 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"), что значительно медленнее, но в некоторых случаях точнее.

Answer 4

Предполагая нормальную аппроксимацию для фиксированных эффектов (что было бы возможно и для ограничения), мы можем получить 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

Answer 5

Я предлагаю вам использовать старый добрый lme (в пакете nlme). У него есть ограничение, и, если вам нужно ограничение контрастов, есть ряд вариантов (оцениваемых в моделях, контраст в контрастах, glht в мульткомпьютере).

Почему p-значения и ограничение отсутствуют в lmer: см. Http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html.

Answer 6

Чтобы найти коэффициент, вы можете просто использовать функцию итога 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)

Answer 7

Я бы предложил tab_model() функция от sjPlotпакет в качестве альтернативы. Чистый и читаемый вывод, готовый к уценке. Ссылка здесь и примеры здесь.

Для более склонных к зрению plot_model() из того же пакета тоже может пригодиться.

Альтернативное решение - через parameters пакет с использованием model_parameters() функция.