Search code examples
rrankmixed-modelstukey

Mixed Model [in place of Repeated Measures ANOVA], but need RANKED Tukey results PER sampling date (NOT overall)

I have data from 6 different treatments, with sampling repeated 8 times over 3 years (no missing data points). My individual bins from each treatment are equally split in 7 randomly distributed blocks. To analyze, I am using a Mixed Model (nlme package).

Example of the data : 

bin     treatment   Bloc        date        CONTAM
b1      TR_A            1       t0      4.753038458
b2      TR_A            2       t0      4.709589136
b3      TR_A            3       t0      4.72668357
b4      TR_A            4       t0      4.647430928
b5      TR_A            5       t0      4.670129151
b6      TR_A            6       t0      4.647430928
b7      TR_A            7       t0      4.811256762
b8      TR_B            1       t0      4.551238194
b9      TR_B            2       t0      4.662660293
b10     TR_B            3       t0      4.753038458
b11     TR_B            4       t0      4.69554541
b12     TR_B            5       t0      4.69554541

Packages used:

nlme ; lattice ; nortest ; multcomp

Here is the script I am currently using:

mod.lme=lme(CONTAM~Treatment*date,random=~1|bloc/bin,data=data)
summary(mod.lme)
anova(mod.lme)
summary(glht(mod.lme, linfct=mcp(Treatment = "Tukey")), test = adjusted(type = "bonferroni"))

This works just fine (ANOVA<0.001) but is not giving me the information I need.

-> I obtain an overall Tukey, but since I am dealing with degradation data I expect the treatments to be similar at the beginning and end, but DIFFERENT in the middle.

----> Therefore, I am seeking a test that will give me the RANKED (a, ab, bc...) Tukey results PER sampling date, while taking in consideration that this is a repeated measures model.

Any ideas? :)

Thanks!

FYI, I've attempted the solutions from this question: post hoc test for a two way mixed model anova

1

library(GAD)
snk.test(mod.lme, term="Treatment*date", among="Treatment", within="date")

My inconclusive result : # Error in object$model[, 2:(length(object$x) + 1)] : incorrect number of dimensions

2

This second one gives me a huge output, but not the one I need.

library(lsmeans)
summary(lsmeans(mod.lme, pairwise~Treatment*date), infer=TRUE)
Solution

  • Try this:

    Obtain the LS means:

    library("lsmeans")
mod.lsm <- lsmeans::lsmeans(mod.lme, ~ Treatment * date)

    (Calling lsmeans::lsmeans prevents it from using the same function in the lmerTest package, if it happens to be loaded.)

    List the LS means:

    mod.lsm

    Do pairwise comparisons, separately for each date:

    pairs(mod.lsm, by = "date")

    (The default result shows $t$ tests and adjusted $P$ values using the Tukey method.)