12
$\begingroup$

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)
$\endgroup$
27
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for the reply, that cleared up how the individual coefficients are calculated. One further question though, when the above paragraph stated that they adjusted for residuals is this the same as including the random effects in the coefficients or is that a separate issue? $\endgroup$ – Andrews Oct 30 '14 at 9:27
  • $\begingroup$ You state that coef(summary(fm1))[, "Estimate"] produces the mean of all subjects. I take that to imply the intercept is the mean value of "Reaction", but 'mean(sleepstudy$Reaction) = 298.5079' . This is the interpretation of the intercept in a linear regression and I'm trying to understand the coefficients in mixed effects models and am struggling to get things (like this to add up). Why isn't the intercept 298.5? $\endgroup$ – svannoy Oct 13 '16 at 21:08
  • $\begingroup$ @svannoy Note that the intercept is the estimated value of the dependent variable when all covariates are equal to zero. Here, Days is not centered at the mean, and therefore the intercept of the model is lower than the mean value of Reaction. Furthermore, the regression line is the best linear fit. Hence, the observations are not necessarily on this line. $\endgroup$ – Sven Hohenstein Oct 14 '16 at 7:32
  • $\begingroup$ @SvenHohenstein thanks, for clarifying that. Is it accurate then to say that in OLS regression with just a single categorical predictor variable and deviation coding (contr.sum()) the intercept is precisely the grand mean; whereas in RC regression with a random intercept model the intercept is the estimated grand mean? $\endgroup$ – svannoy Oct 16 '16 at 11:59
  • $\begingroup$ @svannoy Yes, with deviation coding the intercept reflects the Grand Mean in OLS regression. In a linear mixed model, the intercept is indeed the estimated Grand Mean. $\endgroup$ – Sven Hohenstein Oct 16 '16 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.