R
我正在对(lme4
包)中的线性混合效应模型执行事后测试。我正在使用multcomp
包(glht()
函数)来执行事后测试。
我的实验设计是重复测量,具有随机块效应。模型指定为:
mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)
我没有在此处附加我的数据,而是使用包中调用的warpbreaks
数据multcomp
。
data <- warpbreaks
warpbreaks$rand <- NA
我添加了一个额外的随机变量来模仿我的“块”效果:
warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)
这模仿了我的模型:
mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)
我知道“ Additional Multcomp Examples - 2 Way Anova”中的例子。这个例子引导你比较wool
.
如果我想反其道而行之——比较 的 水平之wool
内的水平tension
怎么办?(在我的情况下,这将是在时间水平(三 - 六月、七月、八月)内比较治疗水平(二 - 0、1)。
我已经提出了以下代码来执行此操作,但它似乎不起作用(请参阅下面的错误消息)。
首先,从示例中(有wool
和tension
交换的地方):
tmp <- expand.grid(wool = unique(warpbreaks$wool), tension = unique(warpbreaks$tension))
X <- model.matrix(~ tension * wool, data = tmp)
glht(mod, linfct = X)
Tukey <- contrMat(table(warpbreaks$wool), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(warpbreaks$tension)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[2], rownames(K2), sep = ":")
从这里到底部,我自己的代码:
K3 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[3], rownames(K3), sep = ":")
K <- rbind(K1, K2, K3)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))
> summary(glht(mod, linfct = K %*% X))
Error in summary(glht(mod, linfct = K %*% X)) :
error in evaluating the argument 'object' in selecting a method for function 'summary': Error in K %*% X : non-conformable arguments