[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 13:31:16 CEST 2010
Hello,
I forgot the starting values question ;-)
How about (but please check this is actually appropriate):
## no starting values
fit12b<-lmer(y ~ time + (1| Subject),
data=data.frame(data),verbose = TRUE)
0: 514.30292: 0.942809
1: 466.47314: 0.00000
2: 466.47314: 0.00000
> summary(fit12b)
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.6 466.5
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.0000 0.0000
Residual 2.6983 1.6427
Number of obs: 120, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.22386 0.32786 0.68
time 0.99005 0.01296 76.40
Correlation of Fixed Effects:
(Intr)
time -0.889
>
#### with starting value
fit12c<-lmer(y ~ time + (1|Subject), start = (50),
data=data.frame(data), verbose = TRUE)
0: 378.91223: 50.0000
1: 369.93186: 53.0381
2: 303.05576: 81.8034
3: 221.14500: 138.040
4: 146.49123: 221.786
5: 69.694894: 360.911
6: -5.9883317: 583.179
7: -81.547428: 942.393
8: -156.62461: 1522.07
9: -229.86189: 2445.53
10: -298.29322: 3868.55
11: -359.72808: 6047.84
12: -404.16415: 8866.56
13: -431.83889: 12188.6
14: -443.87889: 14943.8
15: -453.39962: 23592.1
16: -453.41087: 23528.7
17: -453.49005: 22563.6
18: -453.49030: 22592.6
19: -453.49032: 22591.
20: -453.49032: 22591.
21: -453.49032: 22591.
> summary(fit12c)
Linear mixed model fit by REML
Formula: y ~ time + (1 | Subject)
Data: data.frame(data)
AIC BIC logLik deviance REMLdev
-445.5 -434.3 226.7 -465.3 -453.5
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 5.4560e+02 23.3581476
Residual 1.0691e-06 0.0010340
Number of obs: 120, groups: Subject, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 45.0002761 3.6932472 12
time -1.0000095 0.0001156 -8651
Correlation of Fixed Effects:
(Intr)
time -0.001
>
HTH
Cheers,
Luca
----- Original Message -----
From: "John Sorkin" <jsorkin at grecc.umaryland.edu>
To: <r-sig-mixed-models at r-project.org>; "Luca Borger" <lborger at uoguelph.ca>
Sent: Wednesday, September 15, 2010 7:03 AM
Subject: Re: [R-sig-ME] Incorrect slope, 0.99 vs. correct slope -0.99,due to
false convergence
> Luca,
> First my thanks for your thoughts.
> You are correct that increasing the number of measurements each subject
> has from three to ten eliminates the false convergence. While your
> solution works, it changes the problem and does not answer my original
> question. I have only three observations per subject, and I would like
> to know if it is possible to modify the initial parameters used in the
> iterative solution used by lme so that the false convergence does not
> occur. The essence of my problem is, I believe, that the variance of the
> intercepts is much larger than the residual variance leading to bad
> initial parameter estimates and false convergence. I believe that if I
> can modify the starting parameters I can get true convergence.
> Thank you,
> John
>
> John Sorkin M.D., Ph.D.
> Chief, Biostatistics and Informatics
> Baltimore VA Medical Center GRECC,
> University of Maryland School of Medicine Claude D. Pepper OAIC,
> University of Maryland Clinical Nutrition Research Unit, and
> Baltimore VA Center Stroke of Excellence
>
> 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)
> jsorkin at grecc.umaryland.edu
>
>>>> "Luca Borger" <lborger at uoguelph.ca> 9/14/2010 10:44 PM >>>
> 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
>>
>
>
> Confidentiality Statement:
> This email message, including any attachments, is for ...{{dropped:7}}
More information about the R-sig-mixed-models
mailing list