[R-sig-ME] random effect variance per treatment group in lmer
Afshartous, David
afshart at exchange.sba.miami.edu
Thu Jul 12 23:30:25 CEST 2007
Simon,
Thanks for your extensive comments. Please see my replies below.
I checked out all lmer calls and so far it seems that none
achieve the desired variance stratification in the desired manner.
Part of the issure may rest in not wanting to have a random
effect for the drug term (see below). If I'm missing something
incredibly simple I apologize in advance to the list.
Dave
> -----Original Message-----
> From: Simon Blomberg [mailto:s.blomberg1 at uq.edu.au]
> Sent: Wednesday, July 11, 2007 8:10 PM
> To: Afshartous, David; r-sig-mixed-models at r-project.org;
> r-sig-mixed-models at r-project.org
> Cc: Andrew Robinson
> Subject: Re: [R-sig-ME] random effect variance per treatment
> group in lmer
>
> 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.
>
upon further thought, this is not precisely the model I want since this
model treats the drug shift from the intercept as random per patient,
and I want this to be a fixed effect only.
However, as the random effect on this shift its own variance, this model
seems to implicitly stratify the random effect variance on the intercept
per drug. I.e., there is patient level variability around an intercept
term (representing the reference level of drug), and there is a separate
patient level variability around the drug slope, representing the shift
to the
next level of drug. but once again, I'd rather not have a random
effect for the drug term.
> 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)
>
BTW, this model estimates okay but has the following problem w/ invoking
coef():
Error in coef(fm.no.cov.2) : unable to align random and fixed effects
>
> > > 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.
>
This lmer call still doesn't model b_ij ~ N(0, tau_i), i.e., more than
one variance. (BTW, I assume that the "drug:Patient" can be replaced
by "Patient" when patients only receive 1 drug, as both versions
produced
identical results for the pseudo data below where that is the case).
> > >
> > > 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.
>
Although I'm still not quite sure this model can be
re-written as such, this model doesn't seem to stratify the
random effect variance as desired. There is a random effect on
the intercept for every patient (once again, "drug:Patient" can be
replaced by "Patient" for pseudo data below), and there is a random
effect
on the intercept for every drug, but the latter's probability
distribution does not have its variance depend on drug level.
> 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.
>
What does one do if the data is from a crossover study and
indeed patients 1,2,3 exist in both treatment 1 and treatment 2?
> 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