我有一个BACI试验,涉及三个站点在2007年和2011年之前和之后的测量。我在设置对比时遇到了问题,以便我能够看到每个站点的治疗效果(与对照组相比)的年份效应是否不同。
我尝试过扩展vignette中的示例,但没有成功。我也尝试过lsmeans包。
下面是示例数据框和代码,以及我想测试的明确对比。
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's 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")
ls <- lsmeans(mod1, specs = ~Year*Treatment|Site); pairs(tt, interaction = FALSE)
提供了正确的结构,但会测试每个年份和治疗方法的组合。 - Walden