这种方法通常被认为是“过时的”，所以虽然它可能是可能的，但语法很困难，我怀疑很少有人知道如何操纵 anova 命令来获得你想要的东西。更常见的方法是使用glht来自nlme或的基于似然的模型lme4。（不过，我当然欢迎被其他答案证明是错误的。）

也就是说，如果我需要这样做，我不会打扰 anova 命令。我只需使用 拟合等效模型lm，为这种对比选择正确的误差项，然后自己计算 F 检验（或等效地，t 检验，因为只有 1 个 df）。这要求一切都平衡并具有球形度，但如果您没有，那么您可能应该使用基于可能性的模型。您可以使用 Greenhouse-Geiser 或 Huynh-Feldt 校正在一定程度上校正非球形度，这些校正（我相信）使用相同的 F 统计量，但修改了误差项的 df。

如果你真的想使用car，你可能会发现heplot vignettes 很有帮助；它们描述了car包中的矩阵是如何定义的。

使用 caracal 的方法（对于对比 1&2 - 3 和 1&2 - 4&5），我得到

      psiHat      tStat          F         pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777  0.3719508 5.445988e-01

这就是我获得相同 p 值的方式：

将数据重新整形为长格式并运行lm以获取所有 SS 术语。

library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable), 
                   levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)

为小时项制作一个备用对比矩阵。

foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0

检查我的对比是否提供与默认对比相同的 SS（并且与完整模型中的相同）。

anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))

获取我想要的两个对比的 SS 和 df。

anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))

获取 p 值。

> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999

可选择调整球形度。

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)

如果您想要/必须使用与相应 ANOVA 中的合并误差项的对比，您可以执行以下操作。不幸的是，这会很长，我不知道如何更方便地做到这一点。尽管如此，我认为结果是正确的，因为它们已针对 Maxwell & Delaney 进行了验证（见下文）。

您想hour在 SPF-p.qr 设计中比较您的第一个内因子组（来自Kirk (1995) 的符号：Split-Plot-Factorial design 1 between factor treatmentwith p groups, first in factor hourwith q groups, second within factor prePostFupwith组）。以下假设相同大小的treatment组和球形。

Nj    <- 10                                             # number of subjects per group
P     <- 3                                              # number of treatment groups
Q     <- 5                                              # number of hour groups
R     <- 3                                              # number of PrePostFup groups
id    <- factor(rep(1:(P*Nj), times=Q*R))                                  # subject
treat <- factor(rep(LETTERS[1:P], times=Q*R*Nj), labels=c("CG", "A", "B")) # treatment
hour  <- factor(rep(rep(1:Q, each=P*Nj), times=R))                         # hour
ppf   <- factor(rep(1:R, each=P*Q*Nj), labels=c("pre", "post", "fup"))     # prePostFup
DV    <- round(rnorm(Nj*P*Q*R, 15, 2), 2)               # some data with no effects
dfPQR <- data.frame(id, treat, hour, ppf, DV)           # data frame long format

summary(aov(DV ~ treat*hour*ppf + Error(id/(hour*ppf)), data=dfPQR)) # SPF-p.qr ANOVA

首先请注意，在对 进行hour平均之后， 的主要效果是相同的prePostFup，因此切换到仅包含treatment和hour作为 IV 的更简单的 SPF-pq 设计。

dfPQ <- aggregate(DV ~ id + treat + hour, FUN=mean, data=dfPQR)  # average over ppf
# SPF-p.q ANOVA, note effect for hour is the same as before
summary(aov(DV ~ treat*hour + Error(id/hour), data=dfPQ))

现在请注意，在 SPF-pq ANOVA 中，hour针对交互作用测试的效果id:hour，即，这种交互作用提供了测试的误差项。hour现在可以通过简单地替换误差项和相应的自由度来测试组的对比，就像在单向受试者之间的方差分析中一样。获得这种交互的 SS 和 df 的简单方法是用 拟合模型lm()。

(anRes <- anova(lm(DV ~ treat*hour*id, data=dfPQ)))
SSE    <- anRes["hour:id", "Sum Sq"]     # SS interaction hour:id -> will be error SS
dfSSE  <- anRes["hour:id", "Df"]         # corresponding df

但是，让我们也在这里手动计算所有内容。

# substitute DV with its difference to cell / person / treatment group means
Mjk   <- ave(dfPQ$DV,           dfPQ$treat, dfPQ$hour, FUN=mean)  # cell means
Mi    <- ave(dfPQ$DV, dfPQ$id,                         FUN=mean)  # person means
Mj    <- ave(dfPQ$DV,           dfPQ$treat,            FUN=mean)  # treatment means
dfPQ$IDxIV <- dfPQ$DV - Mi - Mjk + Mj                             # interaction hour:id
(SSE  <- sum(dfPQ$IDxIV^2))               # SS interaction hour:id -> will be error SS
dfSSE <- (Nj*P - P) * (Q-1)               # corresponding df
(MSE  <- SSE / dfSSE)                     # mean square

现在我们有了正确的误差项，我们可以为计划的比较建立通常的检验统计量：其中是对比向量，是它的长度，是对比度估计，是交互作用的均方（合适的误差项）。$t = \frac{\hat{\psi} - 0}{||c|| \sqrt{MS_{E}}}$$c$$||c||$$\hat{\psi} = \sum\limits_{k=1}^{q} c_{k} M_{.k}$$MS_{E}$hour:id

Mj     <- tapply(dfPQ$DV, dfPQ$hour, FUN=mean)  # group means for hour
Nj     <- table(dfPQ$hour)                      # cell sizes for hour (here the same)
cntr   <- rbind(c(1, 1, -2,  0, 0),
                c(1, 1, -1, -1, 0))             # matrix of contrast vectors
psiHat <- cntr   %*% Mj                         # estimates psi-hat
lenSq  <- cntr^2 %*% (1/Nj)                     # squared lengths of contrast vectors
tStat  <- psiHat / sqrt(lenSq*MSE)              # t-statistics
pVal   <- 2*(1-pt(abs(tStat), dfSSE))           # p-values
data.frame(psiHat, tStat, pVal)

对于多重比较，您必须考虑校正方法，例如 Bonferroni。$\alpha$

Maxwell & Delaney (2004) 示例的相应计算见第 3 页。599f 可以在这里找到。请注意，M&D 计算 F 值，要查看结果是否相同，您必须对 t 统计量的值进行平方。该代码还包括使用Anova()from完成的分析car，以及对内因子主效应的$\hat{\epsilon}$