High correlation among random effects for longitudinal model
Dear Joshua, I'm chiming in late, so it's possible that someone already pointed this out and I didn't notice. A better way to specify a polynomial in R is to use poly() in the model formula. By default, this produces orthogonal polynomial regressors (at least in the fixed effects) but the same fit to the data. For example,
time <- 1:5 X <- poly(time, 2)
X
1 2 [1,] -0.6324555 0.5345225 [2,] -0.3162278 -0.2672612 [3,] 0.0000000 -0.5345225 [4,] 0.3162278 -0.2672612 [5,] 0.6324555 0.5345225 attr(,"coefs") attr(,"coefs")$alpha [1] 3 3 attr(,"coefs")$norm2 [1] 1 5 10 14 attr(,"degree") [1] 1 2 attr(,"class") [1] "poly" "matrix"
colSums(X)
1 2 0.000000e+00 1.110223e-16
crossprod(X)
1 2 1 1.000000e+00 -1.110223e-16 2 -1.110223e-16 1.000000e+00 My guess is that this will also reduce the correlations among the random effects. If you really must have raw polynomials, then poly(time, 2, raw=TRUE) offers the advantage that model-structure-aware functions can understand that the linear and quadratic regressors are part of the same term in the model. Whether high correlations among the random effects are really a problem, my guess is that they aren't, because lme() uses a log-Cholesky factorization of the random-effects covariance matrix anyway. In some contexts, high correlations might produce numerical instability, but, as I said, probably not here. Ben would know. I hope this helps, John ----------------------------- John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada Web: socialsciences.mcmaster.ca/jfox/
-----Original Message-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org]
On Behalf Of Joshua Rosenberg
Sent: Monday, April 2, 2018 5:33 PM
To: Ben Bolker <bbolker at gmail.com>
Cc: R SIG Mixed Models <r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] High correlation among random effects for longitudinal
model
Dear Stuart and Ben,
Thank you, this worked to significantly reduce the correlations between the
intercept and the linear and quadratic terms (though still quite high between the
linear and quadratic term):
Random effects:
Formula: ~time + I(time^2) | student_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 18.671959 (Intr) time
time 11.029842 -0.262
I(time^2) 8.359834 -0.506 0.959
Residual 29.006598
Could I ask if that correlation between the linear (time) and quadratic
I(time^2) terms is cause for concern - and if so, how to think about
(potentially) addressing this?
Josh
On Sun, Apr 1, 2018 at 12:34 PM Ben Bolker <bbolker at gmail.com> wrote:
On Sun, Apr 1, 2018 at 12:20 PM, Stuart Luppescu <lupp at uchicago.edu> wrote:
On Sun, 2018-04-01 at 12:55 +0000, Joshua Rosenberg wrote:
lme(outcome ~ time + I(time^2),
random = ~ time + I(time^2),
correlation = corAR1(form = ~ time | individual_ID),
data = d_grouped)
I have a question / concerns about the random effects, as they are
highly correlated (intercept and linear term = -.95; intercept and
quadratic term = .96; linear term and quadratic term = -.995):
I think this is an ordinary occurrence for the intercept and time
trend to be negatively correlated. The way to avoid this is to
center the time variable at a point in the middle of the series, so,
instead of setting the values of time to {0, 1, 2, 3, 4} use {-2, -1, 0, 1, 2}.
Agreed. This is closely related, but not identical to, the phenomenon where the *fixed effects* are highly correlated.
-- Stuart Luppescu Chief Psychometrician (ret.) UChicago Consortium on School Research http://consortium.uchicago.edu
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Joshua Rosenberg, Ph.D. Candidate Educational Psychology ?&? Educational Technology Michigan State University http://jmichaelrosenberg.com [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models