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/
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
