Как указать конкретные контрасты для повторных измерений ANOVA с использованием автомобиля?

Я пытаюсь выполнить повторные измерения Anova в R с последующими некоторыми конкретными контрастами в этом наборе данных. Я думаю, что правильный подход будет использовать Anova()из пакета автомобилей.

Позволяет проиллюстрировать мой вопрос на примере взят из ?Anovaиспользуя OBrienKaiserданные (Примечание: Я опущен гендерный фактор из примера):
Мы имеем конструкцию с одним между субъектами фактором, лечением (3 уровня: контроль, A, B), и 2 повторными -измерения (в рамках предметов) факторы, фаза (3 уровня: предварительный тест, пост-тест, наблюдение) и час (5 уровней: от 1 до 5).

Стандартная таблица ANOVA задается как (в отличие от примера (Anova) я переключился на тип 3 Суммы квадратов, это то, чего хочет мое поле):

require(car)
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)
av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)
summary(av.ok, multivariate=FALSE)

Теперь представьте, что взаимодействие наивысшего порядка было бы значительным (что не так), и мы хотели бы изучить его далее со следующими контрастами:
есть ли разница между часами 1 и 2 по сравнению с часами 3 (контраст 1) и между часами 1 и 2 по сравнению с 4 и 5 часами (контраст 2) в условиях лечения (A и B вместе)?
Другими словами, как мне указать эти контрасты:

1. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) против ((treatment %in% c("A", "B")) & (hour %in% 3))
2. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) против ((treatment %in% c("A", "B")) & (hour %in% 4:5))

Моя идея состояла бы в том, чтобы запустить другую ANOVA, пропуская ненужное условие лечения (контроль):

mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")
av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)
summary(av2, multivariate=FALSE)

Тем не менее, я до сих пор не знаю, как настроить соответствующую контрастную матрицу внутри объекта, сравнивая часы 1 и 2 с 3 и 1 и 2 с 4 и 5. И я не уверен, что если вы пропустите ненужную группу лечения, это действительно хорошая идея, так как она меняет общий термин ошибки.

Прежде чем идти, Anova()я тоже думал о том, чтобы пойти lme. Однако существуют небольшие различия в значениях F и p между ANOVA из учебника и тем, что возвращается из- anove(lme) за возможных отрицательных отклонений в стандартном ANOVA (которые не допускаются вlme ). Относительно того, кто-то указал мне, glsчто позволяет подогнать повторные измерения ANOVA, однако, это не имеет контрастного аргумента.

Чтобы уточнить: я хочу F или t-тест (используя суммы квадратов типа III), который отвечает, являются ли желаемые контрасты значительными или нет.

Обновить:

Я уже задавал очень похожий вопрос по R-help, ответа не было .

Подобные вопросы были заданы на R-help некоторое время назад. Тем не менее, ответы также не решили проблему.

Обновление (2015):

Поскольку этот вопрос все еще порождает некоторую активность, определение тезисов и, в основном, всех других контрастов теперь может быть сделано относительно легко с помощью afexпакета в сочетании с lsmeansпакетом, как описано в виньетке afex .

Этот метод обычно считается «старомодным», поэтому, хотя это и возможно, синтаксис сложен, и я подозреваю, что меньше людей знают, как манипулировать командами anova, чтобы получить то, что вы хотите. Более распространенным методом является использование glhtмодели на основе вероятности из nlmeили lme4. (Я, конечно, могу быть ошибочным при других ответах.)

Тем не менее, если бы мне нужно было сделать это, я бы не стал беспокоиться о командах anova; Я просто подгонял бы эквивалентную модель, используя lm, выбирал правильный термин ошибки для этого контраста и сам вычислял F-тест (или, эквивалентно, t-тест, поскольку есть только 1 df). Это требует, чтобы все было сбалансировано и имело сферичность, но если у вас его нет, вам, вероятно, все равно следует использовать модель, основанную на вероятности. Возможно, вы сможете несколько исправить несферичность, используя поправки Гринхауса-Гейзера или Хюйна-Фельдта, которые (я полагаю) используют ту же F-статистику, но изменяют df члена ошибки.

Если вы действительно хотите использовать car, вы можете найти полезные виньетки heplot ; они описывают, как определяются матрицы в carпакете.

Используя метод Каракала (для контрастов 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))

Получить р-значения.

> 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, вы можете сделать следующее. К сожалению, это будет долго, и я не знаю, как это сделать удобнее. Тем не менее, я думаю, что результаты верны, поскольку они проверены по Максвеллу и Делани (см. Ниже).

Вы хотите сравнить группы первых внутри фактора hourв схеме SPF-p.qr (запись Кирка (1995): схема сплит-график -факториал 1 между фактором treatmentс p группами, сначала внутри фактора hourс q группами, затем внутри фактора prePostFupс р группы). Следующее предполагает одинаковые размеры 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, таким образом, переключается на более простую схему SPF-pq, которая содержит только treatmentи hourкак IV.

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групп можно проверить так же, как в одностороннем ANOVA между субъектами, просто подставив термин ошибки и соответствующие степени свободы. Самый простой способ получить 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)

$\alpha$

Anova()car$\hat{\epsilon}$