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

[Douglas Bates]

On 7/13/07, Afshartous, David <afshart at exchange.sba.miami.edu> 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 expect that the problem is due to "AsIs" terms in the model frame
not being identified correctly in later function calls that build
model matrices.  Parsing model formulas and building the model frame
and model matrices in general settings is a bit of a black art.

The easiest way out is to add variables to the dat.new data frame. Say

dat.new$Dind <- as.numeric(dat.new$drug == "D")
dat.new$Pind <- as.numeric(dat.new$drug == "P")

(fm.het <- lmer(dv ~ ... + (0+Dind|Patient)+(0+Pind|Patient))

I'm afraid I don't understand what the expression

rep(1:n.timepoints, n.subj.per.tx*2)*drug

would mean.  It seems that the variable drug is a factor and it
doesn't make sense to multiply a factor by an integer variable.

I hope this helps,
Doug Bates

> 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
> >
> >
> >
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>