[R] lme random intercept model vs. random intercept, random slope model yield difference in fixed effects but not difference in the anova.
John Sorkin
JSorkin at grecc.umaryland.edu
Sun Apr 15 05:02:35 CEST 2012
I am running two mixed effects regressions. One (fitRandomIntercept) has a random intercept, the second (fitRandomInterceptSlope) has a random intercept and a random slope. In the random slope regression, the fixed effect for time is not significant. In the random intercept random slope model, the fixed effect for time is significant. Despite the difference in the results obtained from the two models. a comparison of the two models, performed via anova(fitRandomIntercept,fitRandomInterceptSlope), shows that there is no significant difference between the two models. I don't understand how this can happen, and I don't know which model I should report. The random intercept random slope model makes physiologic sense, but the principle of parsimony would suggest I report the random intercept model given that it is simpler than the random intercept random slope model, and there is no significant difference between the two models.
Can someone help me understand (1) why one model has a significant slope where as other does not, and (2) given the difference in the two models why the ANOVA comparison of the two model is not significant.
Thanks,
John
Log and code follows:
> # Define data.
> line <- c(1,2,6,11,12,16,17,18,19,21,22,23,24,25,26,31,32,33,34,35,36,41,42,43,46,47,48,49,51,52,56,57,61,66,67,71,72,73,77,82,87,92,97,107,112,117)
> subject<- c(1,1,2,3,3,4,4,4,4,5,5,5,5,5,6,7,7,7,7,7,8,9,9,9,10,10,10,10,11,11,12,12,13,14,14,15,15,15,16,17,18,19,20,22,23,24)
> time <- c(1,3,1,1,6,1,3,7,4,1,3,7,3,35,1,1,3,10,2,25,1,1,3,9,1,3,9,2,1,6,1,3,1,1,3,1,3,11,7,7,7,7,7,7,7,6)
> value <- c(22,4,39,47,5,34,3,33,21,42,9,86,56,39,57,71,8,57,62,47,79,60,10,68,47,6,46,48,57,11,85,12,34,30,1,42,7,33,1,1,1,1,1,1,1,2)
> # Add data to dataframe.
> repeatdatax <- data.frame(line=line,subject=subject,time=time,value=value)
> # Print the data.
> repeatdatax
line subject time value
1 1 1 1 22
2 2 1 3 4
3 6 2 1 39
4 11 3 1 47
5 12 3 6 5
6 16 4 1 34
7 17 4 3 3
8 18 4 7 33
9 19 4 4 21
10 21 5 1 42
11 22 5 3 9
12 23 5 7 86
13 24 5 3 56
14 25 5 35 39
15 26 6 1 57
16 31 7 1 71
17 32 7 3 8
18 33 7 10 57
19 34 7 2 62
20 35 7 25 47
21 36 8 1 79
22 41 9 1 60
23 42 9 3 10
24 43 9 9 68
25 46 10 1 47
26 47 10 3 6
27 48 10 9 46
28 49 10 2 48
29 51 11 1 57
30 52 11 6 11
31 56 12 1 85
32 57 12 3 12
33 61 13 1 34
34 66 14 1 30
35 67 14 3 1
36 71 15 1 42
37 72 15 3 7
38 73 15 11 33
39 77 16 7 1
40 82 17 7 1
41 87 18 7 1
42 92 19 7 1
43 97 20 7 1
44 107 22 7 1
45 112 23 7 1
46 117 24 6 2
>
> # Run random effects regression, with random intercept.
> library(nlme)
>
> #random intercept
> fitRandomIntercept <- lme(value~time,random=~1 |subject,data=repeatdatax)
> summary(fitRandomIntercept)
Linear mixed-effects model fit by REML
Data: repeatdatax
AIC BIC logLik
432.7534 439.8902 -212.3767
Random effects:
Formula: ~1 | subject
(Intercept) Residual
StdDev: 5.78855 25.97209
Fixed effects: value ~ time
Value Std.Error DF t-value p-value
(Intercept) 31.70262 5.158094 22 6.146189 0.0000
time -0.26246 0.632612 22 -0.414888 0.6822
Correlation:
(Intr)
time -0.611
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.0972719 -0.9923743 0.1216571 0.6734183 2.0290540
Number of Observations: 46
Number of Groups: 23
>
> #random intercept and slope
> fitRandomInterceptSlope <- lme(value~time,random=~1+time|subject,data=repeatdatax)
> summary(fitRandomInterceptSlope)
Linear mixed-effects model fit by REML
Data: repeatdatax
AIC BIC logLik
434.7684 445.4735 -211.3842
Random effects:
Formula: ~1 + time | subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05423581 (Intr)
time 2.05242164 -0.477
Residual 23.70228346
Fixed effects: value ~ time
Value Std.Error DF t-value p-value
(Intercept) 38.85068 5.205499 22 7.463392 0.0000
time -2.45621 1.081599 22 -2.270903 0.0333
Correlation:
(Intr)
time -0.648
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.30668304 -0.63797284 -0.09662859 0.69620396 2.05363165
Number of Observations: 46
Number of Groups: 23
>
> # Compare random intercept model to random intercept and slope model.
> anova(fitRandomIntercept,fitRandomInterceptSlope)
Model df AIC BIC logLik Test L.Ratio
fitRandomIntercept 1 4 432.7534 439.8902 -212.3767
fitRandomInterceptSlope 2 6 434.7684 445.4735 -211.3842 1 vs 2 1.985043
p-value
fitRandomIntercept
fitRandomInterceptSlope 0.3706
>
My code :
# Define data.
line <- c(1,2,6,11,12,16,17,18,19,21,22,23,24,25,26,31,32,33,34,35,36,41,42,43,46,47,48,49,51,52,56,57,61,66,67,71,72,73,77,82,87,92,97,107,112,117)
subject<- c(1,1,2,3,3,4,4,4,4,5,5,5,5,5,6,7,7,7,7,7,8,9,9,9,10,10,10,10,11,11,12,12,13,14,14,15,15,15,16,17,18,19,20,22,23,24)
time <- c(1,3,1,1,6,1,3,7,4,1,3,7,3,35,1,1,3,10,2,25,1,1,3,9,1,3,9,2,1,6,1,3,1,1,3,1,3,11,7,7,7,7,7,7,7,6)
value <- c(22,4,39,47,5,34,3,33,21,42,9,86,56,39,57,71,8,57,62,47,79,60,10,68,47,6,46,48,57,11,85,12,34,30,1,42,7,33,1,1,1,1,1,1,1,2)
# Add data to dataframe.
repeatdatax <- data.frame(line=line,subject=subject,time=time,value=value)
# Print the data.
repeatdatax
# Run random effects regression, with random intercept.
library(nlme)
#random intercept
fitRandomIntercept <- lme(value~time,random=~1 |subject,data=repeatdatax)
summary(fitRandomIntercept)
#random intercept and slope
fitRandomInterceptSlope <- lme(value~time,random=~1+time|subject,data=repeatdatax)
summary(fitRandomInterceptSlope)
# Compare random intercept model to random intercept and slope model.
anova(fitRandomIntercept,fitRandomInterceptSlope)
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}
More information about the R-help
mailing list