I would like to extract the slopes for each individual in a mixed effect model, as outlined in the following paragraph
Mixed effects models were used to characterize individual paths of change in the cognitive summary measures, including terms for age, sex, and years of education as fixed effects (Laird and Ware, 1982; Wilson et al., 2000, 2002c)…. Residual, individual cognitive decline slope terms were extracted from the mixed models, after adjustment for the effects of age, sex, and education. Person-specific, adjusted residual slopes were then used as a quantitative outcome phenotype for the genetic association analyses. These estimates equate to the difference between an individual’s slope and the predicted slope of an individual of the same age, sex, and education level.
De Jager, P. L., Shulman, J. M., Chibnik, L. B., Keenan, B. T., Raj, T., Wilson, R. S., et al. (2012). A genome-wide scan for common variants affecting the rate of age-related cognitive decline. Neurobiology of Aging, 33(5), 1017.e1–1017.e15.
I have looked at using the
coef
function to extract the coefficients for each individual, but I am unsure if this is the correct approach to be using.Can anyone provide some advice on how to do this?
#example R code library(lme4) attach(sleepstudy) fml <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) beta <- coef(fml)$Subject colnames(beta) <- c("Intercept", "Slope") beta summary(beta) summary(fm1)
Answer
The model:
library(lme4)
data(sleepstudy)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
The function coef
is the right approach for extracting individual
differences.
> coef(fm1)$Subject
(Intercept) Days
308 253.6637 19.6662581
309 211.0065 1.8475834
310 212.4449 5.0184067
330 275.0956 5.6529540
331 273.6653 7.3973908
332 260.4446 10.1951151
333 268.2455 10.2436611
334 244.1725 11.5418622
335 251.0714 -0.2848735
337 286.2955 19.0955694
349 226.1950 11.6407008
350 238.3351 17.0814915
351 255.9829 7.4520286
352 272.2687 14.0032989
369 254.6806 11.3395025
370 225.7922 15.2897513
371 252.2121 9.4791308
372 263.7196 11.7513155
These values are a combination of the fixed effects and the variance components
(random effects). You can use summary
and coef
to obtain the coefficients
of the fixed effects.
> coef(summary(fm1))[ , "Estimate"]
(Intercept) Days
251.40510 10.46729
The intercept is 251.4 and the slope (associated with Days
) is 10.4.
These coeffcients are the mean of all subjects. To obtain the random effects,
you can use ranef
.
> ranef(fm1)$Subject
(Intercept) Days
308 2.2585637 9.1989722
309 -40.3985802 -8.6197026
310 -38.9602496 -5.4488792
330 23.6905025 -4.8143320
331 22.2602062 -3.0698952
332 9.0395271 -0.2721709
333 16.8404333 -0.2236248
334 -7.2325803 1.0745763
335 -0.3336936 -10.7521594
337 34.8903534 8.6282835
349 -25.2101138 1.1734148
350 -13.0699598 6.6142055
351 4.5778364 -3.0152574
352 20.8635944 3.5360130
369 3.2754532 0.8722166
370 -25.6128737 4.8224653
371 0.8070401 -0.9881551
372 12.3145406 1.2840295
These values are the variance components of the subjects. Every row
corresponds to one subject. Inherently the mean of each column is zero since
the values correspond to the differences in relation to the fixed effects.
> colMeans(ranef(fm1)$Subject)
(Intercept) Days
4.092529e-13 -2.000283e-13
Note that these values are equal to zero, deviations are due to imprecision
of floating point number representation.
The result of coef(fm1)$Subject
incoporates the fixed effects into
the random effects, i.e., the fixed effect coefficients are added to
the random effects. The results are individual intercepts and slopes.
Attribution
Source : Link , Question Author : Andrews , Answer Author : Sven Hohenstein