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

Alan Cobo-Lewis alanc at umit.maine.edu
Fri Jul 13 18:36:11 CEST 2007 

"Douglas Bates" <bates at stat.wisc.edu> on Friday, July 13, 2007 at 12:16 PM -0500 wrote:

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

rep(1:n.timepoints, n.subj.per.tx*2) was a cheesy way of turning time back into a quantitative predictor.

Since rep(1:n.timepoints, n.subj.per.tx*2)*drug wasn't wrapped in I(), this expression in the context of a model formula indicates that there's a fixed linear effect of time, a main effect of drug, and an interaction term (in other words letting the
fixed effect of time be different for the "D" condition than it is for the "P" condition).

"Douglas Bates" <bates at stat.wisc.edu> on Friday, July 13, 2007 at 12:16 PM -0500 wrote:

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

Yes, that would be clearer! As would adding a variable to the data frame for quantitative time! This would allow a fairly simple lmer call like

fm.het <- lmer( dv ~ time.num*drug + (0+Dind|Patient) + (0+Pind|Patient), ... )
for a parameterization in terms of an interaction
or
fm.het <- lmer( dv ~ drug/time.num + (0+Dind|Patient) + (0+Pind|Patient), ... )
for a parameterization in terms of the time course being nested within level of drug

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



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