[R-sig-ME] random effect variance per treatment group in lmer
Simon Blomberg
s.blomberg1 at uq.edu.au
Thu Jul 12 02:09:59 CEST 2007
On Thu, 2007-07-12 at 06:57 +1000, Andrew Robinson wrote:
> 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 )
The above model specifies a random intercept with one random effect per
Patient, and a random slope term for drug, with 1 random effect per
patient. Covariances of the random effects for intercept and drug are
estimated.
This is the model with zero covariances for the random effects, with
Patient as the single level of grouping:
fm.cov <- lmer(z ~ drug + time + (1|Patient) + (0 + 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
If your notation is correct, then this is the lmer call:
fm <- lmer(z ~ drug + time + (1|drug:Patient), data=dat.new)
So you get different random effects on the intercept for each drug *
Patient combination. you can estimate one variance of these random
effects.
> >
> > Case 2:
> > Y_ijk = mu + alpha_i + Indicator_treat_i * b_treatment_ij +
> > Indicator_placebo_i * b_placebo_ij + theta_k +
> > espilon_ijk
Hold on, I think the above model can be rewritten as:
Y_ijk = mu + alpha_i + Indicator_i * b1_i + Indicator_ij * b2_ij +
theta_k + epsilon_ijk
fm <- lmer(z ~ drug + time + (1|drug) + (1|drug:Patient), data=dat.new)
Here we have 2 levels of grouping of random effects on the intercept: at
the drug level (b1), and at the drug*patient level (or equivalently,
Patient within drug level (b2)). So two variances are estimated: for b1
and b2. So to get the total random effect for each patient, just sum the
appropriate random effects across the grouping levels.
The only trick with lmer (compared to lme) is that the Patient j's
should have unique identifiers. Don't have Patients 1,2,3 for within
treatment 1 and 1,2,3 for patients within treatment 2. Use 1,2,3 for
treatment 1 and 4,5,6 for treatment 2 etc.
I hope I have now understood your problem correctly!
Simon.
> >
> > 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.
>
--
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.
More information about the R-sig-mixed-models
mailing list