[R-sig-eco] Mixed Effects Model with Post-hoc interaction contrasts (lme + glht)

[ONKELINX, Thierry]

Dear Alyse,

The easiest way it to create a new variable that has the interaction. E.g. apa3$TimeReach <- interaction(apa3$Time, apa3$Reach). Then refit  your model with this variable instead of Time and Reach.
lme(APA~ TimeReach, random=~1|Station, method="REML", data=apa3)

The coefficient can be different due to the change in parameterization, but the model fit will be identical.

Now you can use glht with user defined contracts to make the required comparisons. See the examples in ?glht on how to do that.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

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

-----Oorspronkelijk bericht-----
Van: r-sig-ecology-bounces at r-project.org [mailto:r-sig-ecology-bounces at r-project.org] Namens Yeager, Alyse D
Verzonden: maandag 7 april 2014 20:11
Aan: r-sig-ecology at r-project.org
Onderwerp: [R-sig-eco] Mixed Effects Model with Post-hoc interaction contrasts (lme + glht)

Hello all,

I am a Master's student seeking help to run a mixed effects model with a post-hoc (Tukey) analysis on my thesis data.  I am having difficulty setting contrasts for the post-hoc analysis so that it only compares 3 specific interactions (REF vs TRT @PRE), (REF vs TRT @WK4), and (REF vs TRT @WK8).

I need to compare enzyme activity (APA) from periphyton on rocks in a reference (REF) and treatment (TRT) "reach" of a stream.  I have 3 "stations" within each reach (UP, MID, DWN) where 4 replicate rocks were collected, and I have 3 "time" points to observe differences due to the treatment (PRE, WK4, WK8). Thus, 48 samples total.  I have set up a mixed effects model using "reach" and "time" as the fixed effects and "station" as the random effect.

My mixed effects model was run as follows using the "nlme" package:
lme(APA~1+Time+Reach+Time*Reach, random=~1|Station, method="REML", data=apa3)

Here is the output:
Linear mixed-effects model fit by REML
 Data: apa3
        AIC      BIC   logLik
  -302.3262 -284.809 159.1631

Random effects:
 Formula: ~1 | Station
         (Intercept)   Residual
StdDev: 3.833552e-07 0.01938121

Fixed effects: APA ~ 1 + Time + Reach + Time * Reach
                        Value   Std.Error DF   t-value p-value
(Intercept)       0.028333333 0.005594873 64  5.064160  0.0000
TimeWK4          -0.003333333 0.007912346 64 -0.421283  0.6750
TimeWK8          -0.003333333 0.007912346 64 -0.421283  0.6750
ReachTRT          0.003333333 0.007912346 64  0.421283  0.6750
TimeWK4:ReachTRT  0.023333333 0.011189747 64  2.085242  0.0410 TimeWK8:ReachTRT  0.025833333 0.011189747 64  2.308661  0.0242
 Correlation:
                 (Intr) TimWK4 TimWK8 RchTRT TWK4:R
TimeWK4          -0.707
TimeWK8          -0.707  0.500
ReachTRT         -0.707  0.500  0.500
TimeWK4:ReachTRT  0.500 -0.707 -0.354 -0.707 TimeWK8:ReachTRT  0.500 -0.354 -0.707 -0.707  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.14984854 -0.73094850 -0.08599394  0.42996971  2.36483339

Number of Observations: 72
Number of Groups: 3

>From this output, I gathered that there was a significant interaction between Reach and Time.

The issue I am attempting to address is the manner in which I setup a post-hoc test comparing only the two reaches on one specific time point.  We are not interested in comparing, for instance, REF at WK4 to TRT at WK8, but only (REF & TRT @PRE), (REF & TRT @WK4), and (REF & TRT @WK8).  I have been attempting to set up this comparison using glht in the "multcomp" package, but cannot figure out how to prepare a contrast statement with these interactions only.

Anyone who can set up an interaction contrast statement, please help.

Thank you,
Alyse

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.