Skip to content

Growth Curve Modeling in NY Times

9 messages · Roberts, Kyle, Douglas Bates, Afshartous, David +2 more

Original
Douglas Bates
#
On 7/9/07, Roberts, Kyle <kyler at mail.smu.edu> wrote:
I would say that it is definitely lmer-related.  It is for exactly
this type of model that I worked hard to ensure that lmer could handle
large data sets and models with crossed or partially crossed random
effects.
Afshartous, David
# 
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 )
Andrew Robinson
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Simon Blomberg
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      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:

            

      
      
      
    

        

      
      
    

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

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

    

      
      #
    

  


  

      
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 Blomberg
    

    

      
      #
    

  


  

      
On Thu, 2007-07-12 at 06:57 +1000, Andrew Robinson wrote:

            

      
      
      
    

        
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)

            

      
      
      
    

        
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.

            

      
      
      
    

        
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.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Afshartous, David
    

    

      
      #
    

  


  

      
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

            

      
      
      
    

        
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.

            

      
      
      
    

        
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

            

      
      
      
    

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

            

      
      
      
    

        
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.

            

      
      
      
    

        
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?