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

[Andrew Robinson]

Dave,

I don't feel that I am sufficiently well informed about the
conventions in lmer to comment.  It could work that way.  I suggest
that you try some simulations, if you are not convinced by the
solution suggested offline.  

Cheers,

Andrew

On Wed, Jul 11, 2007 at 11:23:39AM -0400, Afshartous, David wrote:
>  
> Simon, Andrew:
> 
> Thanks for the replies.
> I am not interested in stratifying the variance of the innermost
> residuals,
> but rather the variance of the random effects, viz., b_ij (drug i,
> patient j) 
> is a random variable w/ variance depending on i.  
> 
> Possible solution suggested offline for previously supplied pseudo data:
> 
> fm.cov =  lmer(z ~ drug + time + (drug|Patient), data = dat.new ) 
> OR,
> fm.no.cov  =  lmer(z ~ drug + time + (0 + drug|Patient), data = dat.new
> ) 
> 
> Formally, consider:
> 
> Case 1:
> Y_ijk = mu + alpha_i + b_ij + theta_k + espilon_ijk 
> alpha = fixed effect for group, theta = fixed effect for time, 
> b = random effect per patient; b_ij ~ N(0, tau_i)  ## variance of random
> effect depends on treatment
> 
> Case 2: 
> Y_ijk = mu + alpha_i + Indicator_treat_i * b_treatment_ij + 
> 		Indicator_placebo_i * b_placebo_ij + theta_k +
> espilon_ijk
> 
> Indicator_treat_i = 1 if i is in treatment group, 0 otherwise 
> Indicator_placebo_i = 1 if i is in placebo group, 0 otherwise
> 
> where b_treatment_ij and b_placebo_ij are different random effects
> terms, with
> different variances; only one will apply per patient equation as per the
> indicator 
> variables.  The cumbersome notation allows for a covariance since we now
> have "two" random effects. (although it seem nonsensical to want such a 
> covariance)
> 
> Does fm.no.cov estimates Case 1 model and fm.cov estimates Case 2 model?
> 
> Cheers,
> Dave
> 
> 
> 
> 
> 
> 
> 
> 
> -----Original Message-----
> From: Simon Blomberg [mailto:s.blomberg1 at uq.edu.au] 
> Sent: Wednesday, July 11, 2007 1:58 AM
> To: Andrew Robinson
> Cc: Afshartous, David; r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] random effect variance per treatment group in
> lmer
> 
> I think he is asking to stratify the variance of the innermost
> residuals, or at least it's not clear. In lme that can be accomplished
> with weights=varFixed(~1|Patient).
> 
> To stratify at different levels of nesting, say the data is this:
>  dat <- data.frame(inner=rep(1:10, each=5), outer=rep(1:2, each=25),
> x=rnorm(50))
> 
> Then this call to lme does the job:
> 
>  fit <- lme(x ~ 1, random=list(outer=~1, inner=~1), data=dat,
> weights=varComb(varIdent(form=~1|outer), varIdent(form=~1|inner)))
> 
> edited output:
> 
> Combination of variance functions: 
>  Structure: Different standard deviations per stratum
>  Formula: ~1 | outer
>  Parameter estimates:
>         1         2 
> 1.0000000 0.5170794
>  Structure: Different standard deviations per stratum
>  Formula: ~1 | inner
>  Parameter estimates:
>         1         2         3         4         5         6         7
> 8
> 1.0000000 0.3127693 0.4475444 0.7323698 0.3647991 0.5962917 1.4127508
> 1.7664527 
>         9        10 
> 0.9475334 0.3666155 
> 
> Cheers,
> 
> Simon.
> 
> weights=varOn Wed, 2007-07-11 at 15:04 +1000, Andrew Robinson wrote:
> > Hi David,
> > 
> > as far as I am aware, there is no option for stratifying the variance 
> > of random effects in either lme or lmer.  One can stratify the 
> > variance of the innermost residuals in lme, but that is different than
> 
> > what you are asking for.
> > 
> > Cheers,
> > 
> > Andrew
> > 
> > 
> > On Tue, Jul 10, 2007 at 10:23:21AM -0400, Afshartous, David wrote:
> > > 
> > > 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 
> > > )
> > > 
> > > _______________________________________________
> > > R-sig-mixed-models at r-project.org mailing list 
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > 
> --
> Simon Blomberg, BSc (Hons), PhD, MAppStat. 
> Lecturer and Consultant Statistician
> Faculty of Biological and Chemical Sciences The University of Queensland
> St. Lucia Queensland 4072 Australia Room 320 Goddard Building (8)
> T: +61 7 3365 2506
> email: S.Blomberg1_at_uq.edu.au
> 
> Policies:
> 1.  I will NOT analyse your data for you.
> 2.  Your deadline is your problem.
> 
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of
> data. - John Tukey.

-- 
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/