[R-sig-ME] random effect variance per treatment group in lmer

Afshartous, David afshart at exchange.sba.miami.edu
Sat Jul 14 01:18:31 CEST 2007 

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








 

> -----Original Message-----
> From: Afshartous, David 
> Sent: Friday, July 13, 2007 4:32 PM
> To: Alan Cobo-Lewis
> Cc: r-sig-mixed-models at r-project.org; Afshartous, David; 
> Andrew Robinson
> Subject: RE: random effect variance per treatment group in lmer
> 
> 
> Alan,
> 
> Thanks again for the help.  I defined time.num, Dind, and 
> Pind outside the lmer call as per Doug Bates' suggestion and 
> now it works and I get the same results as your output w/ the 
> suggested lmer call:
> 
> ( fm.het <- lmer( dv ~ time.num*drug + ( 0 + Dind | Patient ) 
> + ( 0 + Pind | Patient ), data=dat.new ) )
> 
> More importantly, this achieves the goal I was seeking 
> originally.  Looks like the key to the issue is not defining 
> a random effect via "(1 | Patient)" 
> as once this is done it doesn't seem possible to stratify; 
> but rather splitting 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.
> 
> Thanks again to you and everyone else for the help!
> 
> Cheers,
> Dave
> 
> 
> 
> 
> 
> 
>  
> 
> > -----Original Message-----
> > From: Alan Cobo-Lewis [mailto:alanc at umit.maine.edu]
> > Sent: Friday, July 13, 2007 12:27 PM
> > To: Afshartous, David
> > Cc: r-sig-mixed-models at r-project.org; " 
> > "Afshartous at basalt.its.maine.edu; Andrew Robinson
> > Subject: Re: random effect variance per treatment group in lmer
> > 
> > 
> > Hmm, could it be a word-wrap issue? I just verified that the code 
> > works on R 2.5 with lme4 and lattice packages installed.
> > See http://www.umaine.edu/visualperception/lme4het
> > 
> > The model I posted assumes that the subj-to-subj variability in 
> > baseline score has higher variability in the "D" condition 
> than in the 
> > "P" condition. You are correct that there's no random 
> effect for drug. 
> > I don't think you *want* there to be a random effect for 
> drug, since 
> > drug doesn't have levels sampled from some population (it's 
> properly a 
> > fixed effect).
> > Instead, the code I posted has the random effect of Patient with a 
> > higher variability in one level of drug than in the other level of 
> > drug.
> > 
> > If you're looking some the drug to have a different time course for 
> > some patients than for other patients I think that'd be a random 
> > effect of Time. (To avoid overparameterization I think you'd almost 
> > surely want to treat time as a quantitative predictor if 
> you're going 
> > to model it as having a random effect.)
> > 
> > One of the reasons that I constructed a full example with a new sim 
> > was to get a big-enough data set together so that the lmer 
> fit fairly 
> > obviously matched the parameters of the sim (to verify the 
> correctness 
> > of the sim and the suggested model formula). Another reason is that 
> > the code of your sim didn't actually have different subj-to-subj 
> > variability in the "D"
> > condition than it did in the "P" condition, so a model with such an 
> > effect was overparameterized for your original simulated 
> data. Could 
> > this be why the lmer call I suggested failed for you?
> > 
> > BTW, nice to see someone from the UM School of Business (I didn't 
> > notice until just now what your email address said). I worked there 
> > about 20 years ago on the Applied AI Reporter when I was an 
> > undergraduate psychology major. But, looking at the SBA web site, I 
> > recognize barely anyone, and those I do recognize probably wouldn't 
> > recognize me.
> > 
> > alan
> > 
> > "Afshartous, David" <afshart at exchange.sba.miami.edu> on 
> Friday, July 
> > 13, 2007 at 10:44 AM -0500 wrote:
> > > 
> > >Alan,
> > >
> > >Thanks for the suggestion.  I noticed the following error
> > msg for that
> > >lmer call:
> > >
> > >> ( fm.het <- lmer( dv ~ rep(1:n.timepoints,
> > n.subj.per.tx*2)*drug + (
> > >> 0
> > >+ as.numeric(drug=="D") | Patient ) + ( 0 + as.numeric(drug=="P") |
> > >Patient ), data=dat.new ) )
> > >Error: length(term) == 3 is not TRUE
> > >
> > >I tried a few changes to the model but the error still 
> exists; I'll 
> > >keep checking.  I assume the rationale for the structure of
> > your lmer
> > >call, where you use as.numeric as opposed to just drug 
> above, is to 
> > >insure that you do not introduce a random effect for drug into the 
> > >model?
> > >
> > >Regards,
> > >Dave
> > >
> > >
> > >
> > >> -----Original Message-----
> > >> From: Alan Cobo-Lewis [mailto:alanc at umit.maine.edu]
> > >> Sent: Wednesday, July 11, 2007 6:40 PM
> > >> To: r-sig-mixed-models at r-project.org
> > >> Cc: " "Afshartous at basalt.its.maine.edu; Afshartous, 
> David; Andrew 
> > >> Robinson
> > >> Subject: Re: random effect variance per treatment group in lmer
> > >> 
> > >> 
> > >> Dave,
> > >> 
> > >> How about using stratifying variance on level of drug using ( 0 +
> > >> as.numeric(drug=="D") | Patient ) + ( 0 +
> > >> as.numeric(drug=="P") | Patient ) Here's some code 
> (whose sim also 
> > >> builds in a fixed effect of time that applies only to the Drug 
> > >> condition).
> > >> 
> > >> 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 ( fm.het <- lmer( dv ~ 
> > >> rep(1:n.timepoints, n.subj.per.tx*2)*drug + ( 0 +
> > >> as.numeric(drug=="D") | Patient ) + ( 0 +
> > >> as.numeric(drug=="P") | Patient ), data=dat.new ) )
> > >> 
> > >> 
> > >> HTH,
> > >> alan
> > >> 
> > >> 
> > >> >From: "Afshartous, David" <afshart at exchange.sba.miami.edu>
> > >>  asked:
> > >> 
> > >> >
> > >> >All,
> > >> >I didn't receive a response to the query below sent to
> > the general
> > >> >R-help mailing list so figured I'd try this mailing list.
> > >> Apologies in
> > >> >advance if this is an overly simplistic question for this
> > >> list; I just
> > >> >started w/ lmer after not using lme for awhile.
> > >> >Cheers,
> > >> >Dave
> > >> >
> > >> >
> > >> >__________________________________________________________
> > >> >
> > >> >All,
> > >> > 
> > >> >How does one specify a model in lmer such that say the
> > random effect
> > >> >for
> > >> >
> > >> >the intercept has a different variance per treatment group?  
> > >> >Thus, in the model equation, we'd have say b_ij represent
> > the random
> > >> >effect for patient j in treatment group i, with variance
> > >> depending on
> > >> >i, i.e,
> > >> >var(b_ij) = tau_i.
> > >> > 
> > >> >Didn't see this in the docs or Pinherio & Bates (section 5.2 is 
> > >> >specific for modelling within group errors).  Sample
> > >> repeated measures
> > >> >code below is for a single random effect variance, where
> > the random
> > >> >effect corresponds to patient.
> > >> >cheers,
> > >> >dave
> > >> > 
> > >> > 
> > >> >z <- rnorm(24, mean=0, sd=1)
> > >> >time <- factor(paste("Time-", rep(1:6, 4), sep="")) Patient
> > >> <- rep(1:4,
> > >> >each = 6) drug <- factor(rep(c("D", "P"), each = 6, times =
> > >> 2)) ## P =
> > >> >placebo, D = Drug dat.new <- data.frame(time, drug, z,
> > >> Patient) fm =
> > >> >lmer(z ~ drug + time + (1 | Patient), data = dat.new )
> > >> 
> > >> 
> > >> 
> > >> --
> > >> Alan B. Cobo-Lewis, Ph.D.		(207) 581-3840 tel
> > >> Department of Psychology		(207) 581-6128 fax
> > >> University of Maine
> > >> Orono, ME 04469-5742     		alanc at maine.edu
> > >> 
> > >> http://www.umaine.edu/visualperception
> > >> 
> > >> 
> > >> 
> > >
> > 
> > 
> > 
> > --
> > Alan B. Cobo-Lewis, Ph.D.		(207) 581-3840 tel
> > Department of Psychology		(207) 581-6128 fax
> > University of Maine
> > Orono, ME 04469-5742     		alanc at maine.edu
> > 
> > http://www.umaine.edu/visualperception
> > 
> > 
> >

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