[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 03:05:08 CEST 2010


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}}




More information about the R-sig-mixed-models mailing list