[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