0

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)
K_R
  • 29
  • 1
  • 2
  • 6
  • In the `lsmeans` call, specify `pairwise ~ Treatment | date`. If you look at `vignette("using-lsmeans")`, you can find out all sorts of things you can do. – Russ Lenth Mar 31 '17 at 23:28

1 Answers1

0

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.)

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • This is great in the sense that it does indeed split the results by date!:) >Unfortunately, they are not ranked, making the results difficult to interpret with my 6 treatments (15lines/date=120lines to interpret).
    Example of the result:
    date = t0:
    contrast estimate SE df t.ratio p.value
    TR_B-TR_A -0.221575173 0.06832637 30 -3.243 0.0313
    TR_C-TR_A -0.170258301 0.06832637 30 -2.492 0.1586
    TR_D-TR_A -0.142952834 0.06832637 30 -2.092 0.3182
    Would you happen to have a way to organize this with a ranking method? THANKS! :)     – K_R Apr 06 '17 at 00:32
  • Sorry for the long paragraph! I haven't figured out how to code "enter/carriage return" yet. Hope it's clear enough? – K_R Apr 06 '17 at 00:41
  • I can't really read your comment on my phone. But perhaps you want to try this: `cld(mod.lsm)`. It creates a compact letter display, which displays the LS means in numerical order for each "by" level, along with grouping letters. To use this you need to also have the **multcompView** package installed. – Russ Lenth Apr 06 '17 at 12:43
  • We are getting somewhere! The cld(mod.lsm) code works, but it returns ranks comparing the 120 lines of data to each other, rather than grouping them by "Treatment" within the dates. Nonetheless, I can go through and extract the information myself within the dates easier than before. Thanks again! – K_R Apr 07 '17 at 02:00