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. 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.