# 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

...    ...       ...     ...   ...
```

```
+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[],"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
```

```
+3

source

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.

0

source

All Articles