I am currently working with linear mixed models (LMM) and want to perform some residual analyses. Since the residuals from linear mixed models "are correlated and do not necessarily have constant variance", they are inadequate to check the assumptions underlying the model (Fitzmaurice et al. 2011, p. 266). It has thus been suggested (Waternaux et al. 1989) to use Cholesky transformed residuals instead to investigate the model assumptions, see also Houseman et al. (2004) and Santos Nobre & da Motta Singer (2007).
That is, suppose we have a LMM of the form Yi = Xiβ + Zibi + Ɛi, where β is a vector of fixed effect parameters, bi is a vector of normally distributed random effects, and where Xi and Zi are design matrices.
Then, the residuals ei = Yi - Xiβ are correlated, as are their observed counterpart using estimated β-coefficients instead of the unknown true values. We denote those residuals as ri. Decorrelating the residuals is done through finding the Cholesky decomposition of the covariance matrix of ri. That is, letting Cov(ei) = Σi, Let Lii' be the Cholesky decomposition of the estimate of Σi based upon ri. The Cholesky transformed residuals are then defined as r*i = Li-1ri. These are the residuals I am trying to get from the lme4 (Bates et al. 2015) package but so far without success.
However, using data from the homepage of Fitzmaurice et al. (2011), I have managed to recreate the residuals from a case study in the book, using the nlme package (as opposed to the lme4 package). Code is provided below.
require(foreign) # To read dta-files
require(nlme) # The nlme package used to estimate the LMM.
require(mgcv) # Needed for the extract.lme.cov function
fat <- read.dta("https://content.sph.harvard.edu/fitzmaur/ala2e/fat.dta") # We read the data
# Add spline function for time, with a knot at age of menarch, denoted time = 0 in the data.
fat$timeknot <- fat$time * (fat$time > 0)
# Estimate the LMM with nlme
nlme.model <- lme(pbf ~ time + timeknot, data=fat, random = ~ 1 + time + timeknot | id)
# Calculate the residuals
raw.residuals <- residuals(nlme.model, level=0)
# Cholesky residuals
est.cov <- extract.lme.cov(nlme.model, fat) # We extract the blocked covariance function of the residuals
Li <- t(chol(est.cov)) # We find the Cholesky transformation of the residuals. (The transform is to get the lower trangular matrix.)
cholesky.residuals <- solve(Li) %*% raw.residuals # We then calculate the transformed residuals.
hist(cholesky.residuals) # We plot the residuals
My question is thus, is there an easy way to extract the estimated covariance matrix of the residuals, using the lme4 package? Or to get the Cholesky transformed residuals directly?
Thankful for any help or pointers. With regards,
/Phil
References:
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48.
Fitzmaurice, G. M., Laird, N. M., & Ware, J. H. (2011). Applied longitudinal analysis, Second Edition. John Wiley & Sons.
Houseman, E. A., Ryan, L. M., & Coull, B. A. (2004). Cholesky residuals for assessing normal errors in a linear model with correlated outcomes. Journal of the American Statistical Association, 99(466), 383-394.
Santos Nobre, J., & da Motta Singer, J. (2007). Residual analysis for linear mixed models. Biometrical Journal: Journal of Mathematical Methods in Biosciences, 49(6), 863-875.
Waternaux, C., Laird, N. M., & Ware, J. H. (1989). Methods for analysis of longitudinal data: blood-lead concentrations and cognitive development. Journal of the American Statistical Association, 84(405), 33-41.
