[R-sig-ME] Incorrect slope, 0.99 vs. correct slope -0.99, due to false convergence
John Sorkin
jsorkin at grecc.umaryland.edu
Wed Sep 15 13:03:46 CEST 2010
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\...{{dropped:16}}
More information about the R-sig-mixed-models
mailing list