Recovery factors and their standard error from lme
How could I extract the coefficients (b0 and b1) with their standard errors for each experimental unit (plot) in a linear mixed model such as this one:
with the same dataset (df), and for the model (fitL1): how can I get a dataframe like this ...
plot b0 b0_se b1 b1_se
1 2898.69 53.85 -7.5 4.3
... ... ... ... ...
source to share
The first comment is that this is actually a non-trivial theoretical question: there is a long thread for r-sig-mixed-models that goes into some technical details; you should definitely take a look, although it's a little scary. The main problem is that the estimated values of the coefficients for each group are the sum of the fixed effect parameter and the BLUP / conditional mode for that group, which are different classes of objects (one parameter, one is a conditional average for a random variable), which creates some technical difficulties.
The second point is that (unfortunately) I don't know an easy way to do this in lme
, so my answer is using lmer
(from a package lme4
).
If you are more comfortable doing the lightest thing and ignoring the (possibly undefined) covariance between the fixed effect and BLUP parameters, you can use the code below.
Two alternatives: (1) match your model with a Bayesian hierarchical approach (e.g. batch MCMCglmm
) and calculate the backward prediction standard deviations for each level (2) using parametric loading to compute BLUP / conditional modes, then accept the standard deviations from the loading distributions.
Remember, as usual, this advice does not contain any guarantees.
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
source to share
To get you started there using nlme ...
You can pull out the summary () components using:
summary(fitL1)$tTable[,1] #fixed-effect parameter estimates
summary(fitL1)$tTable[,2] #fixed-effect parameter standard errors
and etc.
You can also multiply them line by line:
summary(fitL1)$tTable[1,1] #the first fixed-effect parameter estimate
summary(fitL1)$tTable[1,2] #the first fixed-effect parameter standard error
to extract individual parameters or standard errors and concatenate them into a data frame using, for example:
df<-data.frame(cbind(summary(fitL1)$tTable[1,1], summary(fitL1)$tTable[1,2]))
names(df)<-c("Estimate","SE")
df
To tweak these parameters for each plot (a random effect I suppose), you can output the random coefficients with:
fitL1$coefficients$random
and add them to the parameter estimates (B0 (intercept), B1, etc.). However, I'm not sure how the standard errors should be adjusted for each plot.
source to share