Как получить «общее» значение p и величину эффекта для категориального фактора в смешанной модели (lme4)?

28

Я хотел бы получить p-значение и величину эффекта независимой категориальной переменной (с несколькими уровнями) - то есть «в целом», а не для каждого уровня в отдельности, как нормальный вывод из lme4R. Это так же, как то, о чем люди сообщают при запуске ANOVA.

Как я могу получить это?

user3288202
источник
Какую статистику вы хотите именно? Вы можете использовать anova()функцию, чтобы получить таблицу anova с линейными смешанными моделями, как с линейными моделями.
Смиллиг
Я пробовал anova (), но он дает мне значения Df, Sum Sq, Mean Sq и F. Я не вижу величины эффекта и значения p. У вас есть идеи по этому поводу?
user3288202
1
Под размером эффекта вы подразумеваете что-то вроде эквивалента ? Что касается p-значений, существует долгая и существенная дискуссия вокруг их оценки и их применения в . Посмотрите обсуждение в этом вопросе для более подробной информации. р2lme4
Смиллиг
Спасибо за ссылку, Смилиг. Означает ли это, что из-за проблемы с вычислением значения p величина эффекта фактора в целом также является проблемой?
user3288202
Они не связаны напрямую. Однако следует помнить, что линейная смешанная модель не ведет себя точно так же, как линейная модель без случайных эффектов, поэтому мера, которая может быть подходящей для линейной модели, не обязательно обобщается на смешанные модели.
Смиллиг

Ответы:

48

Обе концепции, о которых вы упоминаете (p-значения и величины эффектов линейных смешанных моделей), имеют присущие им проблемы. Что касается величины эффекта , цитируя Дуга Бейтса, оригинального автора lme4,

Предполагая, что кто-то хочет определить меру , я думаю, что можно привести аргумент для обработки штрафованной остаточной суммы квадратов из линейной смешанной модели так же, как мы рассматриваем остаточную сумму квадратов из линейной модели. Или можно использовать только остаточную сумму квадратов без штрафа или минимальную остаточную сумму квадратов, получаемую из данного набора терминов, что соответствует матрице с бесконечной точностью. Я не знаю, правда. Это зависит от того, что вы пытаетесь охарактеризовать.р2

Для получения дополнительной информации вы можете посмотреть эту ветку , эту ветку и это сообщение . По сути, проблема заключается в том, что не существует согласованного метода включения и декомпозиции дисперсии от случайных эффектов в модели. Однако есть несколько стандартов, которые используются. Если вы посмотрите на Wiki, настроенный для / в списке рассылки r-sig-mixed-models , есть пара подходов.

Один из предложенных методов рассматривает корреляцию между подобранными и наблюдаемыми значениями. Это может быть реализовано в R, как предложил Джаррет Бернс в одном из этих потоков:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

Например, скажем, мы оцениваем следующую линейную смешанную модель:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

Мы можем рассчитать величину эффекта, используя функцию, определенную выше:

r2.corr.mer(fm1)
# [1] 0.0160841

Подобная альтернатива рекомендуется в статье Ронхуи Сюй , именуемой , и может быть вычислена в R просто:Ω02

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

Что касается p-значений , это гораздо более спорный вопрос (по крайней мере, в R / lme4сообществе). Смотрите обсуждения в вопросах здесь , здесь и здесь, среди многих других. Снова ссылаясь на страницу Wiki, есть несколько подходов для проверки гипотез о влиянии в линейных смешанных моделях. Перечислено от «худшего к лучшему» (согласно авторам вики-страницы, в которую, я полагаю, входят Даг Бейтс, а также Бен Болкер, который вносит здесь большой вклад)

  • Wald Z-тесты
  • Для сбалансированных вложенных LMM, где можно вычислить df: t-тесты Вальда
  • Проверка отношения правдоподобия, либо путем настройки модели так, чтобы параметр можно было изолировать / отбросить (с помощью anovaили drop1), либо с помощью вычисления профилей правдоподобия
  • MCMC или параметрические доверительные интервалы начальной загрузки

Они рекомендуют подход выборки по методу Монте-Карло в цепочке Маркова, а также перечисляют ряд возможностей для реализации этого из псевдо и полностью байесовских подходов, перечисленных ниже.

Псевдо-байесовский:

  • Выборочная выборка, как правило (1) с использованием плоских априорных значений и (2) начиная с MLE, возможно, с использованием приблизительной оценки дисперсии-ковариации для выбора распределения кандидатов
  • Через mcmcsamp(если доступно для вашей проблемы: т.е. LMM с простыми случайными эффектами - не GLMM или сложными случайными эффектами)
    Через pvals.fncв languageRпакете, оболочка для mcmcsamp)
  • В AD Model Builder, возможно, через glmmADMBпакет (используйте mcmc=TRUEопцию) или R2admbпакет (напишите свое собственное определение модели в AD Model Builder) или вне R
  • С помощью simфункции из armпакета (имитирует апостериорный только для коэффициентов бета (фиксированный эффект)

Полностью байесовские подходы:

  • Через MCMCglmmпакет
  • Использование glmmBUGS(оболочка WinBUGS / интерфейс R )
  • Использование JAGS / WinBUGS / OpenBUGS и т. Д. Через пакеты rjags/ r2jags/ R2WinBUGS/BRugs

Для иллюстрации того, как это может выглядеть, ниже приведена MCMCglmmоценка с использованием MCMCglmmпакета, который, как вы увидите, дает результаты, аналогичные приведенной выше модели, и имеет некие байесовские p-значения:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

Надеюсь, это поможет несколько. Я думаю, что лучший совет для тех, кто начинает с линейных смешанных моделей и пытается оценить их в R, - это прочитать Вики-часто задаваемые вопросы, откуда взялась большая часть этой информации. Это отличный ресурс для всевозможных тем со смешанными эффектами, от базовых до продвинутых, от моделирования до черчения.

smillig
источник
Спасибо большое, улыбнись. Поэтому я не могу сообщить размер эффекта для общих параметров.
user3288202
р2
3
+6, впечатляюще понятный, всеобъемлющий и полностью аннотированный.
gung - Восстановить Монику
1
Кроме того, вы можете взглянуть на пакет afex и особенно на смешанную функцию. смотрите здесь
beginR
6

Что касается вычисления значений значимости ( p ), Luke (2016). Оценка значимости в линейных моделях смешанных эффектов в R сообщает, что оптимальным методом является приближение Кенварда-Роджера или Саттервейта для степеней свободы (доступно в R с такими пакетами, как lmerTestили afex).

Аннотация

Модели смешанных эффектов все чаще используются при анализе экспериментальных данных. Однако в пакете lme4 в R стандарты оценки значимости фиксированных эффектов в этих моделях (т. Е. Получения значений p) несколько расплывчаты. Для этого есть веские причины, но поскольку исследователи, использующие эти модели, во многих случаях обязаны сообщать значения p, необходим некоторый метод оценки значимости результатов модели. В этой статье приводятся результаты моделирования, показывающие, что два наиболее распространенных метода оценки значимости, использующие тесты отношения правдоподобия и применение распределения z к значениям Вальда t из выходных данных модели (t-as-z), являются несколько антиконсервативными, особенно для небольших выборок. Другие методы оценки значимости,Результаты этого моделирования показывают, что коэффициенты ошибок типа 1 наиболее близки к 0,05, если модели подбираются с использованием REML, а значения p выводятся с использованием приближений Кенварда-Роджера или Саттертвейта, поскольку оба этих приближения дают приемлемые коэффициенты ошибок типа 1 даже для меньших образцы.

(выделение добавлено)

Пабло Бернабеу
источник
4
+1 Спасибо, что поделились этой ссылкой. Я просто кратко прокомментирую, что приближение Кенварда-Роджера доступно в lmerTestпакете.
говорит амеба, восстанови Монику
5

Я использую lmerTestпакет. Это удобно включает оценку значения p в anova()выходных данных для моих анализов MLM, но не дает размер эффекта по причинам, приведенным в других постах здесь.

Bruna
источник
1
В моем случае я предпочитаю парное сравнение с использованием lsmeans, поскольку оно дает мне все пары контрастов, включая значения p. Если я использую lmerTest, мне придется шесть раз запустить модель с разными базовыми линиями, чтобы увидеть все пары контрастов.
user3288202