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

Afshartous, David afshart at exchange.sba.miami.edu
Fri Jul 13 22:31:40 CEST 2007 

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