Contrast adjustment for tripartite interaction using multcomp

Question

I have a BACI experiment with measures before (2007) and after (2011) on three sites. I'm having trouble adjusting the contrasts, so I can see if the effect of the year affects the cf control on each of the three sites.

I tried to extend this example in vignette without success. Also tried the lsmeans package.

The df example and the code below, as well as the explicit contrasts I would like to check.

library(nlme) library(multcomp) require(lsmeans) # I have these versions: # R version 3.2.4 Revised (2016-03-16 r70336) # Platform: x86_64-w64-mingw32/x64 (64-bit) # Running under: Windows 7 x64 (build 7601) Service Pack 1 # nlme_3.1-126 # multcomp_1.4-4 exdf <- expand.grid(Site = c("Z", "Y", "X"), Plot = 1:30, Year = c("2007", "2011"), Treatment = c("Control", "A", "B", "A+B")) exdf$Site <- factor(exdf$Site) exdf$Plot <- factor(exdf$Plot) exdf$Year <- factor(exdf$Year) exdf$Treatment <- factor(exdf$Treatment) exdf$Treatment <- relevel(exdf$Treatment, ref = "Control") exdf$Response <- rnorm(mean = 1150, sd = 130, n = nrow(exdf)) mod1 <- lme(data = exdf, fixed = Response ~ Year * Treatment * Site, random = ~1|Plot) ## would like to test (obviously that not actual code) ## How do the treatments A, B, A+B differ in their interaction with ## year, compared to the control treatment (ie do the slopes differ) SiteX:Year:TreatmentControl - SiteX:Year:TreatmentA == 0 SiteX:Year:TreatmentControl - SiteX:Year:TreatmentB == 0 SiteX:Year:TreatmentControl - SiteX:Year:TreatmentA+B == 0 SiteY:Year:TreatmentControl - SiteY:Year:TreatmentA == 0 SiteY:Year:TreatmentControl - SiteY:Year:TreatmentB == 0 SiteY:Year:TreatmentControl - SiteY:Year:TreatmentA+B == 0 SiteZ:Year:TreatmentControl - SiteZ:Year:TreatmentA == 0 SiteZ:Year:TreatmentControl - SiteZ:Year:TreatmentB == 0 SiteZ:Year:TreatmentControl - SiteZ:Year:TreatmentA+B == 0 ## Also potentially - does the interaction between site and year ## differ across sites? SiteX:Year:TreatmentA - SiteY:Year:TreatmentA == 0 SiteZ:Year:TreatmentA - SiteY:Year:TreatmentA == 0 SiteZ:Year:TreatmentA - SiteX:Year:TreatmentA == 0 SiteX:Year:TreatmentB - SiteY:Year:TreatmentB == 0 SiteZ:Year:TreatmentB - SiteY:Year:TreatmentB == 0 SiteZ:Year:TreatmentB - SiteX:Year:TreatmentB == 0 SiteX:Year:TreatmentA+B - SiteY:Year:TreatmentA+B == 0 SiteZ:Year:TreatmentA+B - SiteY:Year:TreatmentA+B == 0 SiteZ:Year:TreatmentA+B - SiteX:Year:TreatmentA+B == 0 ## Started doing the example from #https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf # Got a bit stuck tmp <- expand.grid(Year = levels(exdf$Year), Treatment = levels(exdf$Treatment), Site = levels(exdf$Site)) X <- model.matrix(~ Year * Treatment * Site, data = tmp) glht(mod1, linfct = X) ## Gives all the interactions not in the form I want.. tt <- lsmeans(mod1, specs = ~Year:Treatment|Site) comps <- pairs(tt, interaction = TRUE) summary(comps, adjust = "holm") 

Answer 1

Build it step by step. First, you want to get a difference of two years for each combination of other factors. This is easy to do in lsmeans using

 tt = lsmeans(mod1, specs = ~ Year | Treatment:Site) dd = pairs(tt, reverse = TRUE) # (2011) - (2008) for each T*S comb summary(dd, by = NULL) 

Now you can compare these differences as you wish:

 pairs(dd, by = "Site") pairs(dd, by = "Treatment") 

... or compare them slightly:

 pairs(lsmeans(dd, "Treatment")) pairs(lsmeans(dd, "Site"))