Hello,
This is a follow up question to my previous one
http://tolstoy.newcastle.edu.au/R/e4/help/08/02/3600.html
I am attempting to model relationship satisfaction (MAT) scores
(measurements at 5 time points), using participant (spouseID) and
couple id (ID) as grouping variables, and time (years) and conflict
(MCI.c) as predictors. I have been instructed to include random
effects for the slopes of both predictors as well as the intercepts,
and then to drop non-significant random effects from the model. The
instructor and the rest of the class is using HLM 6.0, which gives p-
values for random effects, and the procedure is simply to run a model,
note which random effects are not significant, and drop them from the
model. I was hoping I could to something analogous by using the anova
function to compare models with and without a particular random
effect, but I get dramatically different results than those obtained
with HLM 6.0.
For example, I wanted to determine if I should include a random effect
for the variable "MCI.c" (at the couple level), so I created two
models, one with and one without, and compared them:
> m.3 <- lmer(MAT ~ 1 + years + MCI.c + (1 + years | spouseID) + (1 +
years + MCI.c | ID), data=Data, method = "ML")
> m.1 <- lmer(MAT ~ 1 + years + MCI.c + (1 + years + MCI.c |
spouseID) + (1 + years + MCI.c | ID), data=Data, method = "ML")
> anova(m.1, m.3)
Data: Data
Models:
m.3: MAT ~ 1 + years + MCI.c + (1 + years | spouseID) + (1 + years +
m.1: MCI.c | ID)
m.3: MAT ~ 1 + years + MCI.c + (1 + years + MCI.c | spouseID) + (1 +
m.1: years + MCI.c | ID)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
m.3 12 5777.8 5832.7 -2876.9
m.1 15 5780.9 5849.5 -2875.4 2.9428 3 0.4005
The corresponding output from HLM 6.0 reads
Random Effect Standard Variance df Chi-
square P-value
Deviation Component
------------------------------------------------------------------------------
INTRCPT1, R0 6.80961 46.37075 60
112.80914 0.000
YEARS slope, R1 1.49329 2.22991 60
59.38729 >.500
MCI slope, R2 5.45608 29.76881 60
90.57615 0.007
------------------------------------------------------------------------------
To me, this seems to indicate that HLM 6.0 is suggesting that the
random effect should be included in the model, while R is suggesting
that it need not be. This is not (quite) a "why do I get different
results with X" post, but rather an "I'm worried that I might be
doing
something wrong" post. Does what I've done look reasonable? Is there a
better way to go about it?
Thank you very much for reading this, and for any advice.
-Ista
On 4/22/08, Ista Zahn <istazahn at gmail.com> wrote:> Hello, > This is a follow up question to my previous one http://tolstoy.newcastle.edu.au/R/e4/help/08/02/3600.html> I am attempting to model relationship satisfaction (MAT) scores > (measurements at 5 time points), using participant (spouseID) and > couple id (ID) as grouping variables, and time (years) and conflict > (MCI.c) as predictors. I have been instructed to include random > effects for the slopes of both predictors as well as the intercepts, > and then to drop non-significant random effects from the model. The > instructor and the rest of the class is using HLM 6.0, which gives p- > values for random effects, and the procedure is simply to run a model, > note which random effects are not significant, and drop them from the > model. I was hoping I could to something analogous by using the anova > function to compare models with and without a particular random > effect, but I get dramatically different results than those obtained > with HLM 6.0.> For example, I wanted to determine if I should include a random effect > for the variable "MCI.c" (at the couple level), so I created two > models, one with and one without, and compared them:> > m.3 <- lmer(MAT ~ 1 + years + MCI.c + (1 + years | spouseID) + (1 + > years + MCI.c | ID), data=Data, method = "ML") > > m.1 <- lmer(MAT ~ 1 + years + MCI.c + (1 + years + MCI.c | > spouseID) + (1 + years + MCI.c | ID), data=Data, method = "ML") > > anova(m.1, m.3) > Data: Data > Models: > m.3: MAT ~ 1 + years + MCI.c + (1 + years | spouseID) + (1 + years + > m.1: MCI.c | ID) > m.3: MAT ~ 1 + years + MCI.c + (1 + years + MCI.c | spouseID) + (1 + > m.1: years + MCI.c | ID) > Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) > m.3 12 5777.8 5832.7 -2876.9 > m.1 15 5780.9 5849.5 -2875.4 2.9428 3 0.4005> The corresponding output from HLM 6.0 reads> Random Effect Standard Variance df Chi- square P-value > Deviation Component> ------------------------------------------------------------------------------ > INTRCPT1, R0 6.80961 46.37075 60 112.80914 0.000 > YEARS slope, R1 1.49329 2.22991 60 59.38729 >.500 > MCI slope, R2 5.45608 29.76881 60 90.57615 0.007 > ------------------------------------------------------------------------------> To me, this seems to indicate that HLM 6.0 is suggesting that the > random effect should be included in the model, while R is suggesting > that it need not be. This is not (quite) a "why do I get different > results with X" post, but rather an "I'm worried that I might be doing > something wrong" post. Does what I've done look reasonable? Is there a > better way to go about itThe first thing I would try to determine is whether the model fit by HLM is equivalent to the model you have fit with lmer. The part of the HLM model output you have shown lists only variance components. It does not provide covariances or correlations. The lmer model is fitting a 3 by 3 symmetric positive definite variance-covariance matrix with a total of 6 parameters - 3 variances and 3 covariances. It may be that HLM is fitting a simpler model in which the covariances are all zero. The next question would be exactly how HLM is calculating that p-value. I wonder where the 60 degrees of freedom comes from. Do you happen to have 60 couples in the study?