[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