Стандартная ошибка компонента дисперсии на выходе lmer

Мне нужно извлечь компонент дисперсии standard error из вывода lmer .

library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)

Следующее дает оценки компонента дисперсии:

s2 <- VarCorr(model)$Subject[1]

Это НЕ стандартная ошибка дисперсии. И я хочу стандартную ошибку. Как я могу это получить?

РЕДАКТИРОВАТЬ :

Возможно, я не могу объяснить вам, что я имел в виду под стандартной ошибкой компонента дисперсии. Поэтому я редактирую свой пост.

В главе 12 «Эксперименты со случайными факторами» книги Планирование и анализ экспериментов , Дуглас С. Монтгомери , в конце главы пример 12-2 выполнен SAS . В примере 12-2 модель представляет собой двухфакторную факторную модель со случайным эффектом. Выходные данные приведены в таблице 12-17.

введите здесь описание изображения

Я пытаюсь подогнать модель в R на lmer.

library(lme4)
fit <- lmer(y~(1|operator)+(1|part),data=dat)

R-коды для извлечения Estimate, отмеченные цифрой 4 в таблице 12-17:

est_ope=VarCorr(fit)$operator[1]
est_part = VarCorr(fit)$part[1]
sig = summary(fit)$sigma
est_res = sig^2

Теперь я хочу извлечь результаты Std Errors, отмеченные цифрой 5 в таблице 12-17, из вывода lmer.

Огромное спасибо !


r lme4
person user81411    schedule 29.07.2015    source источник
comment
Вы можете опубликовать данные для этого примера? Также было бы полезно знать, для чего вы собираетесь использовать стандартные ошибки (как я указываю ниже, стандартные ошибки оценок дисперсии являются ненадежными показателями неопределенности - интервалы профиля будут лучше)   -  person Ben Bolker    schedule 29.07.2015
comment
Зачем вам нужны стандартные ошибки для параметра, который не симметрично распределен. Вам следует переформулировать вопрос, который вы задаете. Не повторяйте ошибку SAS. Если вам нужен гипотест, используйте функцию анове. Если вы хотите CI, используйте профиль или начальную CI. Есть причины, по которым Imer не дает номер, который вы просите. Хотя Бен расскажет вам, как вы можете его получить. Не позволяйте тому факту, что отчеты SAS или Stata влияют на вас.   -  person pauljohn32    schedule 10.10.2018

Ответы (3)


arrow_upward
8
arrow_downward

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

library("lme4")
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)

(В настоящее время сделать это для оценок REML немного сложнее...)

Извлеките функцию отклонения, параметризованную с точки зрения стандартного отклонения и корреляции, а не с точки зрения факторов Холецкого (обратите внимание, что это внутренняя функция, поэтому нет гарантии, что она будет работать таким же образом в будущем...)

 dd.ML <- lme4:::devfun2(model,useSc=TRUE,signames=FALSE)

Извлеките параметры как стандартные отклонения по исходной шкале:

 vv <- as.data.frame(VarCorr(model)) ## need ML estimates!
 pars <- vv[,"sdcor"]
 ## will need to be careful about order if using this for
 ## a random-slopes model ...

Теперь вычислите матрицу второй производной (Гессиана):

library("numDeriv")
hh1 <- hessian(dd.ML,pars)
vv2 <- 2*solve(hh1)  ## 2* converts from log-likelihood to deviance scale
sqrt(diag(vv2))  ## get standard errors

Это стандартные ошибки стандартных отклонений: удвойте их, чтобы получить стандартные ошибки дисперсий (при преобразовании значения его стандартные ошибки масштабируются в соответствии с производной преобразования).

Я думаю, что это должно сделать это, но вы можете перепроверить это...

person Ben Bolker    schedule 29.07.2015
comment
Не могли бы вы дать мне ссылку на оценку REML. Результат в таблице для стандартной ошибки оценок REML компонента дисперсии. Возможно, поэтому я не получаю ответа. Я правильно указал модель, а также получил REML-оценки компонента дисперсии на VarCorr(model), что соответствует результату в таблице. Было бы очень полезно, если бы вы порекомендовали мне, как вычислить стандартные ошибки стандартных отклонений lmer(Reaction ~ Days + (Days||Subject), sleepstudy) - person user81411; 02.08.2015
comment
Можно ли в lmer рассчитать стандартную ошибку стандартного отклонения оценок REML? - person user81411; 02.08.2015
comment
Если я напишу REML=TRUE и продолжу в том же духе, это не сработает? - person user81411; 03.08.2015
comment
Я думаю, что, вероятно, нет. Я мог бы открыть вопрос по адресу github.com/lme4/lme4/issues... несколько способов сделать это, но нет очень простых решений. - person Ben Bolker; 03.08.2015

arrow_upward
3
arrow_downward

Я не совсем уверен, что вы подразумеваете под «стандартной ошибкой компонента дисперсии». Мое лучшее предположение (на основе вашего кода) заключается в том, что вам нужна стандартная ошибка случайного эффекта. Вы можете получить это с помощью пакета arm:

library(arm)
se.ranef(model)
#$Subject
#    (Intercept)
#308    9.475668
#309    9.475668
#310    9.475668
#330    9.475668
#331    9.475668
#332    9.475668
#333    9.475668
#334    9.475668
#335    9.475668
#337    9.475668
#349    9.475668
#350    9.475668
#351    9.475668
#352    9.475668
#369    9.475668
#370    9.475668
#371    9.475668
#372    9.475668

На самом деле это квадратный корень условной матрицы дисперсии-ковариации случайного эффекта:

sqrt(attr(ranef(model, condVar = TRUE)$Subject, "postVar"))
person Roland    schedule 29.07.2015

arrow_upward

arrow_downward

person    schedule
comment
Запрос as.data.frame для параметра случайных эффектов (ranef) преобразует в длинный формат condsd, соответствующий условному стандартному отклонению. - person Edgar Benitez; 05.02.2021