добывающих коэффициентов и их стандартные ошибки от LME

Question

Как можно извлечь коэффициенты (b0 и b1) с их соответственно стандартными ошибками для каждой экспериментальной единицы (участок) в виде линейной смешанной модели, такие как этот:

Better fits for a linear model

с этим же набором данных (Df), и для подобранной модели (fitL1): как я мог бы получить кадр данных, как это ...

plot b0  b0_se b1 b1_se 
    1 2898.69 53.85 -7.5 4.3 

    ... ...  ...  ... ... 

Answer 1

Первый комментарий, что это на самом деле является не- тривиальный теоретический вопрос: есть крыса ее long thread on r-sig-mixed-models, который входит в некоторые технические детали; вам обязательно нужно взглянуть, хотя это немного страшно. Основная проблема заключается в том, что оценочные значения коэффициентов для каждой группы представляют собой сумму параметра фиксированного эффекта и BLUP/условный режим для этой группы, которые являются разными классами объектов (один параметр, один является условным средним для случайная величина), что создает некоторые технические трудности.

Второй момент заключается в том, что (к сожалению), я не знаю простой способ сделать это в lme, поэтому мой ответ использует lmer (от lme4 пакета).

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

Возможны две альтернативы: (1) для вашей модели с байесовским иерархическим подходом (например, пакет MCMCglmm) и вычисление стандартных отклонений задних предсказаний для каждого уровня (2) с использованием параметрического бутстрапинга для вычисления BLUPs/условных а затем принять стандартные отклонения от загрузочных распределений.

Пожалуйста, помните, что, как обычно, этот совет поставляется без каких-либо гарантий.

library(lme4) 
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) 
cc <- coef(fm1)$Subject 
## variances of fixed effects 
fixed.vars <- diag(vcov(fm1)) 
## extract variances of conditional modes 
r1 <- ranef(fm1,condVar=TRUE) 
cmode.vars <- t(apply(cv <- attr(r1[[1]],"postVar"),3,diag)) 
seVals <- sqrt(sweep(cmode.vars,2,fixed.vars,"+")) 
res <- cbind(cc,seVals) 
res2 <- setNames(res[,c(1,3,2,4)], 
       c("int","int_se","slope","slope_se")) 
##   int int_se  slope slope_se 
## 308 253.6637 13.86649 19.666258 2.7752 
## 309 211.0065 13.86649 1.847583 2.7752 
## 310 212.4449 13.86649 5.018406 2.7752 
## 330 275.0956 13.86649 5.652955 2.7752 
## 331 273.6653 13.86649 7.397391 2.7752 
## 332 260.4446 13.86649 10.195115 2.7752 

Answer 2

Чтобы вы там, используя частично nlme ...

Вы можете тянуть компоненты резюме() с помощью:

summary(fitL1)$tTable[,1] #fixed-effect parameter estimates 
summary(fitL1)$tTable[,2] #fixed-effect parameter standard errors 

т.д.

Вы можете в дальнейшем подмножество то по строкам:

summary(fitL1)$tTable[1,1] #the first fixed-effect parameter estimate 
summary(fitL1)$tTable[1,2] #the first fixed-effect parameter standard error 

извлекать отдельные параметры или стандартные ошибки и объединить их в кадр данных, используя, например:

df<-data.frame(cbind(summary(fitL1)$tTable[1,1], summary(fitL1)$tTable[1,2])) 
names(df)<-c("Estimate","SE") 
df 

Для настройки этих параметров для каждого участка (случайный эффект, я полагаю), вы можете вытащить случайные коэффициенты с:

fitL1$coefficients$random 

и добавить их в оценках параметров (B0 (перехватывать), B1 и т.д.). Однако я не уверен, как стандартные ошибки должны быть скорректированы для каждого сюжета.