[R-sig-ME] stratifying autocorrelation parameter in corAR1() in lme

[Afshartous, David]

All,

In repeated measures over time it is possible that the witin patient 
autocorrelation model differs between treatment and placebo; thus one 
might want to stratify the AR1 parameter (phi) accordingly.  I attempted

to model this with corAR1() within the lme call, but the output still
produces only
one phi estimate.  Didn't see this in P&B, the archives, or docs.  
Is there a simple way to do this w/ lme?  (The simulated data below has 
repeated measures for each patient on both drug and placebo (not a
crossover).)

Thanks,
David

ps - the lme call below also stratifies the variance of the random
effect on the
intercept; for details, search archives for "random effect variance per
treatment group in lmer".


##################################################
## simulated data 
set.seed(500)
n.timepoints <- 8
n.subj.per.tx <- 20
sd.d <- 5; 
sd.p <- 2; 
sd.res <- 1.3
drug <- factor(rep(c("D", "P"), each = n.timepoints, times =
n.subj.per.tx)) 
drug.baseline <- rep( c(0,5), each=n.timepoints, times=n.subj.per.tx ) 
## Patient <- rep(1:(n.subj.per.tx*2), each = n.timepoints) 
dat.new$Patient.new = rep(1:20, each=16)    
## add new patient numer ID to emulate design where drug various within
patient
Patient.baseline <- rep( rnorm( n.subj.per.tx*2, sd=c(sd.d, sd.p) ),
each=n.timepoints ) 
time <- factor(paste("Time-", rep(1:n.timepoints, n.subj.per.tx*2),
sep="")) 
time.baseline <-
rep(1:n.timepoints,n.subj.per.tx*2)*as.numeric(drug=="D")
dv <- rnorm( n.subj.per.tx*n.timepoints*2,
mean=time.baseline+Patient.baseline+drug.baseline, sd=sd.res ) 
dat.new <- data.frame(time, drug, dv, Patient) 
xyplot( dv~time|drug, group=Patient, type="l", data=dat.new)
# fit model treats time as a quantitative predictor 
dat.new$Dind <- as.numeric(dat.new$drug == "D") 
dat.new$Pind <- as.numeric(dat.new$drug == "P")
dat.new$time.num = rep(1:n.timepoints, n.subj.per.tx*2)
########################################################################
## lme calls:
fm1 <- lme( dv ~ time.num*drug,  random = ~ 0 + Pind + Dind |
Patient.new,  data=dat.new)
## no autocorrelation 
fm1.AR1 <- lme( dv ~ time.num*drug,  random = ~ 0 + Pind + Dind |
Patient.new,  data=dat.new, correlation = corAR1(0, ~ 1 | Patient.new) )

## AR1 w/o stratification
fm1.AR1.strat <- lme( dv ~ time.num*drug,  random = ~ 0 + Pind + Dind |
Patient.new,  data=dat.new, correlation = corAR1(0, ~ time.num|
Patient.new/drug ) ) 
## attempted call, but doesn't stratify, i.e., only one phi in output


David Afshartous, PhD
University of Miami
School of Business
Rm KE-408
Coral Gables, FL 33124