[R] lme random intercept model vs. random intercept, random slope model yield difference in fixed effects but not difference in the anova.
Joshua Wiley
jwiley.psych at gmail.com
Sun Apr 15 05:55:37 CEST 2012
Hi John,
R-help is not really for statistical questions (see something like
statsexchange) or mixed models in R (there is a SIG mailing list for
those). Just a note for next time.
1) The estimate for the time effect when it is a fixed effect versus
random effect are different things. The former is the average effect
in the population, the latter something like (in loose terms) the
expectation of the distribution of effects (if you add other
predictors, the conditional expectation). Suppose in each subject
there is a negative relationship between value and time, however, some
subjects have both low values of time and value and others have high
values of time and value. In this case, the population effect could
be positive, but the individual effect negative.
2) First of all, just because a parameter is non significant in one
model and is significant in another does not imply that the models
overall are significantly different. Remember the t-value for time is
testing against 0, but the effect in your first model was not 0, is
comparing the fit of the first model to the second, not to a null
model. Second, your second model includes two additional
parameters---the variance of the slope effect and the covariance of
the slope and intercept. So you have 2 additional parameters
estimated and only ~2 change in deviance, which is clearly not
significantly different.
If you have not already, I would strongly encourage you to graph your data.
library(ggplot2)
## base plot
p <- ggplot(repeatdatax, aes(x = time, y = value)) +
geom_point() +
stat_smooth(method = "lm", se = FALSE)
p ## overall
p + facet_wrap(~ subject) ## individual
The first graph shows the overall association, the latter breaks it
down by subject. If this is your full data, considering you only have
more than 1 observation on 50% of your subjects, I would strongly tend
toward parsimony and only include the random intercept.
Cheers,
Josh
On Sat, Apr 14, 2012 at 8:02 PM, John Sorkin
<JSorkin at grecc.umaryland.edu> wrote:
> 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 at 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.
--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/
More information about the R-help
mailing list