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

[Afshartous, David]

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.