Skip to content

Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?

5 messages · roberto toro, Adam D. I. Kramer, Mark Difford +1 more

roberto toro
# 
Thanks for answering Mark!

I tried with the coding of the interaction you suggested:
But is it normal that the DF are 2303? DF is 2303 even for the estimate of
LobeO that has only 662 values (331 for Tissue=white and 331 for Tissue=grey).
I'm not sure either that Sex, Lobe and Tissue are correctly handled.... why are
there different estimates called Sex:LobeO, Sex:LobeP, etc, and not just
Sex:Lobe as with aov()?. Why there's Tissuew, but not Sex1, for example?

Thanks again!
roberto

ps1. How would you code this with lmer()?
ps2. this is part of the output of mod<-lme:
Linear mixed-effects model fit by REML
 Data: vlt
       AIC      BIC    logLik
  57528.35 57639.98 -28745.17

Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    11294.65

 Formula: ~1 | tfac %in% Subject
        (Intercept) Residual
StdDev:    10569.03 4587.472

Fixed effects: Volume ~ Sex * Lobe * Tissue
                       Value Std.Error   DF    t-value p-value
(Intercept)        245224.61  1511.124 2303  162.27963  0.0000
Sex                  2800.01  1866.312  329    1.50029  0.1345
LobeO             -180794.83  1526.084 2303 -118.46975  0.0000
LobeP             -131609.27  1526.084 2303  -86.23984  0.0000
LobeT              -73189.97  1526.084 2303  -47.95932  0.0000
Tissuew            -72461.05  1526.084 2303  -47.48168  0.0000
Sex:LobeO            -663.27  1884.789 2303   -0.35191  0.7249
Sex:LobeP           -2146.08  1884.789 2303   -1.13863  0.2550
Sex:LobeT            1379.49  1884.789 2303    0.73191  0.4643
Sex:Tissuew          5387.65  1884.789 2303    2.85849  0.0043
LobeO:Tissuew       43296.99  2158.209 2303   20.06154  0.0000
LobeP:Tissuew       50952.21  2158.209 2303   23.60856  0.0000
LobeT:Tissuew      -15959.31  2158.209 2303   -7.39470  0.0000
Sex:LobeO:Tissuew   -5228.66  2665.494 2303   -1.96161  0.0499
Sex:LobeP:Tissuew   -1482.83  2665.494 2303   -0.55631  0.5781
Sex:LobeT:Tissuew   -6037.49  2665.494 2303   -2.26506  0.0236
Mark Difford
# 
Hi Roberto,

It's difficult to comment further on specifics without access to your data
set. A general point is that the output from summary(aov.object) is not
directly comparable with summary(lme.object). The latter gives you a summary
of a fitted linear regression model, not an analysis of variance model, and
what you "see" will depend on what contrasts were in place when the model
was fitted.

If you haven't changed these then they will be so-called treatment
contrasts. What you are seeing for Lobe (which plainly is coded as a factor)
in the output from summary(lme.object) are the regression coefficients for
each level of Lobe relative to its reference/treatment/baseline level, which
is your (Intercept). If you fitted your model with, say, Helmert or
sum-to-zero contrasts then these values would change.

To see what your current reference level is do levels(dataset$Lobe). See
?levels.

What you want to look at to begin with is: anova(lme.object).

HTH, Mark.
roberto toro wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Adam D. I. Kramer
    

    

      
      #
    

  


  

      
On Sun, 14 Sep 2008, roberto toro wrote:

        

            

      
      
      
    

        
lme is basically doing a regression, not an ANOVA as you're used to it. You
may want anova(mod) instead of summary(mod) to see aggregated effects. Or,
you could define contrasts among your levels by assigning to
contrasts(vlt$Lobe), for example.

Also, in the above model, you're only looking at modeling a separate average
volume for each subject-within-tfac; if I read you correctly, you actually
want to model a lobe and tissue effect for each subject for each tfac, in
which case you would want something like what was in my last post.

--Adam

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Mark Difford
    

    

      
      #
    

  


  

      
Hi Roberto,

The other thing you can do --- if you don't wish to step across to lmer(),
where you will be able to exactly replicate the crossed-factor error
structure --- is stay with aov(... + Error()), but fit the factor you are
interested in last. Assume it is Sex. Then fit your model as

aov.model <- aov(Volume ~ Lobe * Tissue * Sex + Error(Subject/(Lobe *
Tissue))

This should give you a so-called "Type-II" test for Sex. You may verify this
by fitting the model without the Error term and using Anova() from the car
package (which does Type-II/III tests) to look at the SS and F values.

I say should, because the only concern I have is whether this procedure is
affected by the presence of an Error term in the model. Establishing this is
beyond my capabilities.

Regards, Mark.

        
roberto toro wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Benoit Boulinguiez
    

    

      
      #
    

  


  

      
Hi all,

I try to use R to compare models of adsorption onto activated carbon
materials.
Three models are compared:

(1)Langmuir
qe~a*b*Ce/(1+b*Ce)

(2)Freundlich
qe=b*Ce^a

(3)Langmuir-Freundlich
qe~a*b*Ce^z/(1+b*Ce^z)

Experimentally I measure Ce and qe, thus I obtain a dozen of paired values
(Ce,qe) which represents an isotherm of adsorption.
I use the nls() function to assess the parameters of the models and
confint() to get their confidence intervals.

To estimate the fit, I plot the residuals for each models and compare the
SSR values. However, as the model (3) has 3 parameters where models (1) and
(2) have 2, I use the standard deviation to estimate models' fitting each
other.

May I use a sort of R2 determined with the cor() function between the data
and the fitted values, to compare the models as well?


Regards/Cordialement


Benoit Boulinguiez