[R] Multiple comparisons with a mixed effects model

[Sean Henderson]

Hello,

I have used R in the past to conduct multiple comparisons on standard linear models, but am a bit confused as to how to go about doing it with a mixed effects model.

I am conducting a bioindication study using carabid beetles in which I have four treatment types (forest harvest types with varying levels of canopy structure retention), and am using canopy closure percent as a covariate in my ANCOVA model. The response variable is species diversity (richness).

I am wanting to use a mixed effects model because there are only two replicates of each treatment type, and within each of those blocks there were 3 sites established - far enough apart that they should be independent based on literature but there is still potentially site-to-site variation that I'd like to account for.

My question is, how would I go about doing a protected t-test on the treatments to see if there are significant differences between them? Is it even possible to do this? I have tried doing it the same way I've done with linear models in the past, but instead of a t-value I get a z-value? 

Below is the code for what I have done.

Thank you very much for any help!

Sean


>
fm1 <- lme(Diversity ~ Canopy + Treatment, data=data, random = ~ 1 | Site)
>
anova(fm1)
            numDF denDF   F-value p-value
(Intercept)     1    18 310.32845  <.0001
Canopy          1    18  28.38084  <.0001
Treatment       3    18   4.46571  0.0164
>
summary(fm1)
Linear
mixed-effects model fit by REML
 Data: data 
       AIC     BIC    logLik
  86.00919 92.2418 -36.00459
 
Random
effects:
 Formula: ~1 | Site
        (Intercept) Residual
StdDev:    1.086125 0.407297
 
Fixed
effects: Diversity ~ Canopy + Treatment 
                        Value Std.Error DF    t-value p-value
(Intercept)         2.9941329 1.3470256 18  2.2227735  0.0393
Canopy              0.0312745 0.0157384 18  1.9871455  0.0623
TreatmentGap        0.4884026 0.6742472 18  0.7243672  0.4782
TreatmentOpen      -1.9764646 0.9016932 18 -2.1919479  0.0418
TreatmentSemi-open
-1.8719043 0.7405578 18 -2.5276952  0.0211
 Correlation: 
                   (Intr) Canopy TrtmnG TrtmnO
Canopy             -0.936                     
TreatmentGap       -0.355  0.116              
TreatmentOpen      -0.772  0.627  0.441       
TreatmentSemi-open
-0.624  0.427  0.499  0.603
 
Standardized
Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-0.5770339
-0.1937237 -0.1051985  0.1503428  0.9151522 
 
Number
of Observations: 23
Number
of Groups: 23 
>
comparisons <- glht(fm1, linfct=mcp(Treatment="Tukey"))
>
summary(comparisons, test=adjusted("none"))
 
         Simultaneous Tests for General Linear
Hypotheses
 
Multiple
Comparisons of Means: Tukey Contrasts
 
 
Fit:
lme.formula(fixed = Diversity ~ Canopy + Treatment, data = data, 
    random = ~1 | Site)
 
Linear
Hypotheses:
                        Estimate Std. Error z
value Pr(>|z|)    
Gap
- Closed == 0         0.4884     0.6742   0.724 0.468840    
Open
- Closed == 0       -1.9765     0.9017  -2.192 0.028383 *  
Semi-open
- Closed == 0  -1.8719     0.7406  -2.528 0.011481 *  
Open
- Gap == 0          -2.4649     0.8549  -2.883 0.003937 ** 
Semi-open
- Gap == 0     -2.3603     0.7108  -3.321 0.000898 ***
Semi-open
- Open == 0     0.1046     0.7453   0.140 0.888435    
---
Signif.
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05
‘.’ 0.1 ‘ ’ 1 
(Adjusted
p values reported -- none method)



      __________________________________________________________________
Make your browsing faster, safer, and easier with the new Internet Explorer® 8. Optimized for Yahoo! Get it Now for Free! at http://downloads.yahoo.com/ca/internetexplorer/