[R-sig-ME] Incorrect slope, 0.99 vs. correct slope -0.99, due to false convergence
Luca Borger
lborger at uoguelph.ca
Wed Sep 15 04:44:35 CEST 2010
Hello,
it appears this is also because you got only three records for each subject.
If you increase it to just 10, you get a slope exactly -1.0 (and no false
convergence warning):
so, first using 3 records per subject, but lmer instead of lme:
fit12c <-lmer(y ~ time + (1| Subject), data=data.frame(data))
> fit12c
Linear mixed model fit by REML
Formula: y ~ time + (1 | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
474.5 485.6 -233.2 457.7 466.5
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.0000 0.0000
Residual 2.6985 1.6427
Number of obs: 120, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.22396 0.32787 0.68
time 0.99004 0.01296 76.40
Correlation of Fixed Effects:
(Intr)
time -0.889
>
# now, 10 records per subject
#######################################
########################
# Create data #
########################
nsub <- 40
data<-matrix(nrow=10*nsub,ncol=3)
dimnames(data)<-list(NULL,c("Subject","time","y"))
# Define time
time<-rep((1:10)/1,nsub)
data[,"time"]<-time
# Define subject
Subject<-rep((1:nsub),10)
data[,"Subject"]<-sort(Subject)
# Define time value
# Add a bit of noise to time values so we can use random intercept and
random time
data[,"time"]<-data[,"Subject"]+data[,"time"]+rnorm(10*nsub,0,0.001)
# Define y value
data[,"y"]<-data[,"Subject"]+rep((10:1)/1,nsub) +rnorm(10*nsub,0,0.001)
data
fit12c <-lmer(y ~ time + (1| Subject), data=data.frame(data))
> fit12c
Linear mixed model fit by REML
Formula: y ~ time + (1 | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
-3190 -3174 1599 -3213 -3198
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 4.9738e+02 22.3020258
Residual 2.2170e-06 0.0014889
Number of obs: 400, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.200e+01 3.526e+00 15
time -1.000e+00 2.592e-05 -38580
Correlation of Fixed Effects:
(Intr)
time 0.000
>
fit12d <-lmer(y ~ time + (1+time| Subject), data=data.frame(data)) #
gives a false convergence warning
> fit12d
Linear mixed model fit by REML
Formula: y ~ time + (1 + time | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
-2567 -2543 1290 -2579 -2579
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 1.8228e+02 13.5010870
time 2.6926e-01 0.5189022 0.009
Residual 2.4337e-06 0.0015600
Number of obs: 400, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 51.99921 2.13471 24.36
time -0.99997 0.08205 -12.19
Correlation of Fixed Effects:
(Intr)
time 0.009
>
fit12f <-lmer(y ~ time + (0+time| Subject), data=data.frame(data)) # no
false convergence warning, but slope estimate is not exactly -1.00 anymore
> fit12f
Linear mixed model fit by REML
Formula: y ~ time + (0 + time | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
2301 2317 -1146 2292 2293
Random effects:
Groups Name Variance Std.Dev.
Subject time 0.69206 0.8319
Residual 10.02260 3.1658
Number of obs: 400, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 27.8111 0.9359 29.716
time -0.3955 0.1393 -2.838
Correlation of Fixed Effects:
(Intr)
time -0.324
>
fit12g <-lmer(y ~ time + (1| Subject) + (0+time| Subject),
data=data.frame(data)) # false convergence warning, but slope exactly -1
> fit12g
Linear mixed model fit by REML
Formula: y ~ time + (1 | Subject) + (0 + time | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
-3189 -3169 1600 -3214 -3199
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 5.2392e+02 2.2889e+01
Subject time 6.0692e-09 7.7905e-05
Residual 2.1488e-06 1.4659e-03
Number of obs: 400, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.200e+01 3.619e+00 14
time -1.000e+00 2.834e-05 -35291
Correlation of Fixed Effects:
(Intr)
time 0.000
>
#############################################
Just my 2 cents, in case this is of any help.
Cheers,
Luca
---------------------------
Luca Börger, PhD
Postdoctoral Research Fellow
Department of Integrative Biology
University of Guelph
Guelph, Ontario, Canada N1G 2W1
office +1 519 824 4120 ext. 52975
lab +1 519 824 4120 ext. 53594
fax: +1 519 767 1656
email: lborger at uoguelph.ca
www.researcherid.com/rid/C-6003-2008
http://uoguelph.academia.edu/LucaBorger
--------------------------------------------------------------------
----- Original Message -----
From: "John Sorkin" <jsorkin at grecc.umaryland.edu>
To: <r-sig-mixed-models at r-project.org>
Sent: Tuesday, September 14, 2010 9:05 PM
Subject: [R-sig-ME] Incorrect slope, 0.99 vs. correct slope -0.99,due to
false convergence
> Windows XP
> R 2.10.0
>
> I am running a random effects model (random intercept) on a data set in
> which each subject has a negative slope (run script below and see graph).
> When I analyze the data using lme (see end of the script) I get a positive
> slope (0.99 estimate for time) when the true slope should be close
> to -1.0. As has been noted before this is because my data are somewhat
> pathological (the correlation of the intercept and time is -0.889). I
> believe the incorrect slope is due to false convergence. Is there some way
> I can specify starting values to the iterative procedure used in lme so as
> to avoid the false convergence and obtain the correct slope?
> Thanks,
> John
>
>
>
>
> library(nlme)
>
> ########################
> # Create data #
> ########################
> nsub <- 40
> data<-matrix(nrow=3*nsub,ncol=3)
> dimnames(data)<-list(NULL,c("Subject","time","y"))
>
> # Define time
> time<-rep((1:3)/1,nsub)
> data[,"time"]<-time
>
> # Define subject
> Subject<-rep((1:nsub),3)
> data[,"Subject"]<-sort(Subject)
>
> # Define time value
> # Add a bit of noise to time values so we can use random intercept and
> random time
> data[,"time"]<-data[,"Subject"]+data[,"time"]+rnorm(3*nsub,0,0.001)
>
> # Define y value
> data[,"y"]<-data[,"Subject"]+rep((3:1)/1,nsub) +rnorm(3*nsub,0,0.001)
> data
>
> # Plot the data
> plot(data[,"time"],data[,"y"],xlab="Time",ylab="y")
> title("Each subject has a negative slope")
>
> # Regression of y on x ignoring the fact that the data come from 10
> different subjects.
> fit0 <- lm(y~time,data=data.frame(data))
> summary(fit0)
>
> # Compute subject specific regressions.
> for (i in 1:nsub){
> fit0 <- lm(y~time,
> data=data.frame(data[data[,"Subject"]==i,c("time","y")]))
> abline(fit0) # Add sub. specific regression lines to plot.
> coefs[i,] <- coef(fit0)
> }
>
>
> ######################################################################
> ######################################################################
> # This is lme give a slope (time) of 0.99. The true slope should be #
> # approx. -1.0. #
> ######################################################################
> fit12c<-lme(y ~ time, random= ~ 1 | Subject, data=data.frame(data))
> summary(fit12c)
>
>
> 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-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list