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}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
