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

56

Я использую lme4 в R, чтобы соответствовать смешанной модели

lmer(value~status+(1|experiment)))

где значение непрерывно, статус и эксперимент являются факторами, и я получаю

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

Как я могу узнать, что влияние статуса является значительным? R сообщает только а не p-значения .Tп

ECII
источник
1
Если исходить из ответов на этот вопрос, то задаешься вопросом, что на самом деле здесь интересует OP: тестирование коэффициентов по отношению к нулевому (ванильный тест проводится в регулярной линейной регрессии по отношению к нулевому H 0 : β = β нуль ) или тестирование по минимизация дисперсии ( F- тест мы получаем от множества типов ANOVA). Эти два стремятся к разным вещам. Просвещающий ответ, хотя и не о моделях со смешанными эффектами, находится здесь . TЧАС0:βзнак равноβзначение NULLF
Firebug

Ответы:

61

По этой теме много информации в FAQ по GLMM . Тем не менее, в вашем конкретном случае, я бы предложил использовать

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)

потому что вам не нужно ничего из того, что lmerпредлагает (более высокая скорость, обработка скрещенных случайных эффектов, GLMM ...). lmeдолжно дать вам точно такой же коэффициент и дисперсии оценок , но и вычислить ДФ и р-значения для вас (что делать имеет смысл в «классический» дизайн , такие , как вы , кажется, есть). Возможно, вы также захотите рассмотреть случайный термин ~status|experiment( учитывающий изменение эффектов состояния между блоками или эквивалентно включающий взаимодействие «статус за экспериментом»). Постеры выше также верны, что ваша tстатистика настолько велика, что ваше p-значение определенно будет <0,05, но я могу представить, что вам нужны «реальные» p-значения.

Бен Болкер
источник
3
Я не знаю об этом ответе. lmerможет так же легко сообщать о тех же самых p-значениях, но не по уважительным причинам. Я предполагаю, что это комментарий, что здесь есть какие-то "реальные" p-значения, которые меня раздражают. Вы можете утверждать, что вы можете найти одну возможную отсечку, и что любая разумная отсечка пропущена. Но вы не можете утверждать, что есть реальное значение p.
Джон
11
Для классического дизайна (сбалансированного, вложенного и т. Д.), Я думаю, я действительно могу утверждать, что существует реальное значение p, то есть вероятность получения оценки бета наблюдаемой величины или больше, если нулевая гипотеза (бета = 0) были ложными ... lme4 не предоставляет эти знаменатели df, я полагаю, потому что в целом сложнее обнаружить из структуры модели lme4, когда указанная модель является такой, где будет работать некоторая эвристика для вычисления классического знаменателя df ...
Бен Болкер
попробуйте summary(m1)вместо этого (я использую это с пакетом nlme)
jena
36

Вы можете использовать пакет lmerTest . Вы просто устанавливаете / загружаете его, и модели lmer расширяются. Так, например,

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)

даст вам результаты с р-значениями. Если p-значения правильные, это немного оспаривается, но если вы хотите их получить, это способ их получить.

pbx101
источник
28

Если вы можете справиться с отказом от p-значений ( и вы должны это сделать ), вы можете вычислить отношение правдоподобия, которое представляло бы вес доказательств влияния статуса, с помощью:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
    formula = value ~ (1|experiment) + status
    , REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
    formula = value ~ (1|experiment)
    , REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))
Майк Лоуренс
источник
16
Обратите внимание, что отношения правдоподобия являются асимптотическими, то есть не учитывают неопределенность в оценке остаточной дисперсии ...
Бен Болкер
5
Я заинтересован в вашей последней строке. Какова интерпретация результата? Есть ли источники, на которые я могу взглянуть на это?
mguzmann
13

Проблема в том, что вычисление p-значений для этих моделей не является тривиальным, см. Здесь обсуждение, поэтому авторы lme4пакета намеренно решили не включать p-значения в выходные данные. Вы можете найти метод расчета их, но они не обязательно будут правильными.

Мишель
источник
9

Подумай, что ты спрашиваешь. Если вы просто хотите узнать, проходит ли общее значение p для эффекта состояния какое-то произвольное значение отсечения, например 0,05, то это просто. Во-первых, вы хотите узнать общий эффект. Вы можете получить это от anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)

Теперь у вас есть значение F. Вы можете взять это и посмотреть в некоторых F таблицах. Просто выберите минимально возможный деном. степени свободы. Отсечка там будет около 20. Ваш F может быть больше, чем это, но я могу ошибаться. Даже если это не так, посмотрите на количество степеней свободы по сравнению с обычным расчетом ANOVA, используя количество экспериментов, которые у вас есть. Придерживаясь этого значения, вы уменьшаете примерно до 5 для отсечки. Теперь вы легко проходите это в своем кабинете. «Истинный» df для вашей модели будет выше, чем это, потому что вы моделируете каждую точку данных, а не агрегированные значения, которые будет моделировать ANOVA.

Если вы на самом деле хотите получить точное значение p, такой вещи не будет, если вы не готовы сделать теоретическое утверждение об этом. Если вы читаете Pinheiro & Bates (2001, и, возможно, еще несколько книг по этой теме ... см. Другие ссылки в этих ответах), и у вас нет аргумента для конкретного df, тогда вы можете использовать это. Но вы все равно не ищете точное значение p. Я упоминаю об этом, потому что поэтому вы не должны сообщать точное значение p, только то, что ваша отсечка пройдена.

Вы должны действительно рассмотреть ответ Майка Лоуренса, потому что сама идея просто придерживаться точки прохода для p-значений в качестве окончательной и наиболее важной информации для извлечения из ваших данных, как правило, ошибочна (но может и не быть в вашем случае, так как мы не на самом деле достаточно информации, чтобы знать). Майк использует любопытную версию вычисления LR, которая может быть интересной, но может быть трудно найти много документации по ней. Если вы посмотрите на выбор и интерпретацию модели с помощью AIC, она вам может понравиться.

Джон
источник
9

Редактировать: этот метод больше не поддерживается в более новых версиях lme4. Используйте пакет lmerTest, как предложено в этом ответе pbx101 .

Существует запись в списке R автором lme4 для почему не отображаются значения р. Вместо этого он предлагает использовать образцы MCMC, что вы делаете с помощью pvals.fnc из пакета languageR:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)

См. Http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf для примера и подробностей.

Джефф
источник
3
lme4 больше не поддерживает это. Этот пост может быть обновлен, чтобы избавить людей от необходимости выяснять это, как я только что сделал.
timothy.s.lau
5

Вы заинтересованы в том, чтобы узнать, оказывает ли комбинированный эффект statusсущественное влияние на value? Если это так, вы можете использовать Anovaфункцию в carпакете (не путать с anovaфункцией в базе R).

dat <- data.frame(
  experiment = sample(c("A","B","C","D"), 264, replace=TRUE), 
  status = sample(c("D","R","A"), 264, replace=TRUE), 
  value = runif(264)   
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)

Посмотрите ?Anovaпосле загрузки carпакета.

smillig
источник
Любая идея, как car::Anova()избежать липких проблем, связанных с вычислением p-значений, которые связывает Мишель?
Майк Лоуренс
Я не знаю, но я полагаю, что он избегает проблем, игнорируя их! Перечитав оригинальный пост, я чувствую, что, возможно, неправильно понял вопрос. Если ОП хочет точные значения p для параметров фиксированных эффектов, он или она в беде. Но если ОП просто хочет знать, значимы ли они, я думаю, что значения t больше, чем какая-либо неопределенность в том, как будет рассчитываться точное значение p. (Другими словами, они значимы.)
smillig
1
Я думаю, что было бы неплохо перенаправить на расчет ANOVA, чтобы выяснить общий эффект статистики, но я не уверен, что сглаживание p-значений - это хорошо. Обычная anovaкоманда даст вам F.
Джон
Я думаю, что это немного липче, чем очевидно. Запуск ANOVA действителен, когда вы хотите минимизировать дисперсию, но из формулировки вопроса я думаю, что OP хочет установить предельный эффект переменных, то есть проверить коэффициенты против нуля.
Firebug
0

Функция pvals.fncбольше не поддерживается lme4. Используя пакет lmerTest package, можно использовать другой метод для вычисления p-значения, такой как приближения Кенварда-Роджера.

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")
Стефано Каччиаторе
источник
0

Простая загрузка пакета afex выведет p-значения в выводе функции lmer из пакета lme4 (вам не нужно использовать afex; просто загрузите его):

library(lme4)  #for mixed model
library(afex)  #for p-values

Это автоматически добавит столбец p-значения к выводу lmer (yourmodel) для фиксированных эффектов.

Марьям Нассери
источник