Skip to content
Prev 249 / 20628 Next

random effect variance per treatment group in lmer

All,

Here is a brief summary of the problem, model, and solution 
as the multiple e-mails have probably made it difficult to read.

Problem: 
In the context of growth curve data with drug/placebo
conditions, the subject to subject variability is greater in 
the drug group than the placebo group.  It is desired to model
this differential variability by having a random effect on the
intercept, where the variance of this random effect depends on drug.
(Time may either be treated as discrete or continuous.)

Model:

fm.het <- lmer( dv ~ time.num*drug + ( 0 + Dind | Patient )  + ( 0 +
Pind | Patient ), data=dat.new ) )

For a data set w/ the aforementioned variability (generated via code
below), this lmer call fits a model where time is continuous and there 
is an interaction between time and drug, i.e., a slope shift.  Thus,
there 
are fixed effects for intercept, time trend, drug (representing a shift
from drug reference level when using treatment contrasts), and slope
shift.
Dind and Pind are indicator variables for drug and placebo.  Note that
one does not define a random effect via "(1 | Patient)" as once this is
done 
it doesn't seem possible to stratify; but rather one splits the "1" into

the two pieces as per the drug dichotomy, and getting separate random 
effects variances for each while mainting drug as a fixed effect.
Examining
ranef(fm.het) shows that this is accomplished.

Solution:
The estimated model from the simulated data clearly accomplished 
the desired stratification:
Random effects:
 Groups   Name Variance Std.Dev.
 Patient  Dind 29.0183  5.3869  
 Patient  Pind  3.9571  1.9893  
 Residual       1.7108  1.3080 


Note that the second term (3.957) does not represent variability around
the fixed effect for the shift from drug to placebo, but rather
variability 
around the fixed effect for the placebo intercept itself (obtained by 
summing the appropriate fixed effects in the output; one can always
remove
the intercept from the lmer call to inrease the correspondence between
the
fixed effect and random effect ouput).

Thanks again to all for the extensive help,

David

code below for the simulated data:
## from Alan Cobo-Lewis e-mail:
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) 
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)
( fm.het <- lmer( dv ~ time.num*drug + ( 0 + Dind | Patient ) + ( 0 +
Pind | Patient ), data=dat.new ) )