Split-Plot ANOVA: сравнение моделей с тестами в R

Как я могу протестировать эффекты в ANOVA с разделенным графиком, используя подходящие сравнения моделей для использования с аргументами Xи в R? Я знаком с Далгаардом и (2007) [1]. К сожалению, это только чистит дизайн Split-Plot. Делать это в полностью рандомизированной схеме с двумя факторами внутри субъекта:Manova.mlm()?anova.mlm

N  <- 20  # 20 subjects total
P  <- 3   # levels within-factor 1
Q  <- 3   # levels within-factor 2
DV <- matrix(rnorm(N* P*Q), ncol=P*Q)           # random data in wide format
id <- expand.grid(IVw1=gl(P, 1), IVw2=gl(Q, 1)) # intra-subjects layout of data matrix

library(car)        # for Anova()
fitA <- lm(DV ~ 1)  # between-subjects design: here no between factor
resA <- Anova(fitA, idata=id, idesign=~IVw1*IVw2)
summary(resA, multivariate=FALSE, univariate=TRUE)  # all tests ...

Следующие сравнения моделей приводят к тем же результатам. Ограниченная модель не включает рассматриваемый эффект, но все другие эффекты того же порядка или ниже, полная модель добавляет рассматриваемый эффект.

anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitA, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                      X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Дизайн Split-Splot с одним фактором внутри и одним фактором:

idB  <- subset(id, IVw2==1, select="IVw1")          # use only first within factor
IVb  <- gl(2, 10, labels=c("A", "B"))               # between-subjects factor
fitB <- lm(DV[ , 1:P] ~ IVb)                        # between-subjects design
resB <- Anova(fitB, idata=idB, idesign=~IVw1)
summary(resB, multivariate=FALSE, univariate=TRUE)  # all tests ...

Это anova()команды для копирования тестов, но я не знаю, почему они работают. Почему тесты следующих сравнений моделей дают одинаковые результаты?

anova(fitB, idata=idB, X=~1, test="Spherical") # IVw1, IVw1:IVb
anova(fitB, idata=idB, M=~1, test="Spherical") # IVb

Два фактора внутри субъекта и один фактор между субъектами:

fitC <- lm(DV ~ IVb)  # between-subjects design
resC <- Anova(fitC, idata=id, idesign=~IVw1*IVw2)
summary(resC, multivariate=FALSE, univariate=TRUE)  # all tests ...

Как мне воспроизвести результаты, приведенные выше, с соответствующими модельными сравнениями для использования с аргументами Xи ? Какова логика этих сравнений моделей?Manova.mlm()

РЕДАКТИРОВАТЬ: Suncoolsu указал, что для всех практических целей данные из этих проектов должны быть проанализированы с использованием смешанных моделей. Тем не менее, я все еще хотел бы понять, как повторить результаты summary(Anova())с anova.mlm(..., X=?, M=?).

[1]: Далгаард, П. 2007. Новые функции для многомерного анализа. Р Новости, 7 (2), 2-7.

XИ в Mосновном определяют две модели , которые вы хотите сравнить, но только с точки зрения последствий внутригрупповых субъекта; затем он показывает результаты для взаимодействия всех эффектов между субъектами (включая перехват) с эффектами внутри объекта, которые изменились между Xи M.

Ваши примеры fitBлегче понять, если мы добавим значения по умолчанию для Xи M:

anova(fitB, idata=idB, M=~1, X=~0, test="Spherical") # IVb
anova(fitB, idata=idB, M=diag(3), X=~1, test="Spherical") # IVw1, IVw1:IVb

Первая модель - это изменение не изнутри предметных эффектов (все имеют одинаковое среднее значение) на разные средние значения для каждого, поэтому мы добавили idслучайный эффект, который является правильным решением для проверки общего перехвата и общего между предметным эффектом. на.

Вторая модель рекламирует id:IVw1взаимодействие, которое является правильной проверкой IVw1и IVw1:IVbусловиями против. Поскольку существует только один внутрисубъектный эффект (с тремя уровнями), значение по умолчанию diag(3)во второй модели будет учитывать его; это было бы эквивалентно бежать

anova(fitB, idata=idB, M=~IVw1, X=~1, test="Spherical") # IVw1, IVw1:IVb

Для вас fitC, я считаю, что эти команды воссоздают Anovaрезюме.

anova(fitC, idata=id, M=~1, X=~0, test="Spherical") #IVb
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitC, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                  X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Теперь, как вы обнаружили, эти команды действительно хитры. К счастью, нет больше причин использовать их больше. Если вы хотите принять сферичность, вы должны просто использовать aovили, для еще более простого синтаксиса, просто использовать lmи вычислять правильные F-тесты самостоятельно. Если вы не желаете принимать сферичность, использование lmeдействительно является подходом, поскольку вы получаете гораздо большую гибкость, чем с поправками GG и HF.

Например, вот aovи lmкод для вашего fitA. Вы должны сначала иметь данные в длинном формате; Вот один из способов сделать это:

library(reshape)
d0 <- data.frame(id=1:nrow(DV), DV)
d0$IVb <- IVb
d0 <- melt(d0, id.vars=c(1,11), measure.vars=2:10)
id0 <- id
id0$variable <- factor(levels(d0$variable), levels=levels(d0$variable))
d <- merge(d0, id0)
d$id <- factor(d$id)

А вот lm andкод AOV`:

anova(lm(value ~ IVw1*IVw2*id, data=d))
summary(aov(value ~ IVw1*IVw2 + Error(id/(IVw1*IVw2)), data=d))

Сплит-сюжеты возникли в сельском хозяйстве, отсюда и название. Но они часто встречаются, и я бы сказал - рабочая лошадка большинства клинических испытаний. Основной график обрабатывается с уровнем одного фактора, в то время как уровни некоторых других факторов могут варьироваться в зависимости от подзаговоров. Конструкция возникает в результате ограничения полной рандомизации. Например: поле может быть разделено на четыре участка. Может быть возможно высаживать различные сорта на участках, но для всего поля можно использовать только один тип полива. Не различие между расколами и блоками, Блоки - это особенности экспериментальных блоков, которые мы можем использовать в экспериментальном проекте, потому что мы знаем, что они есть. Разделение, с другой стороны, накладывает ограничение на то, какие присвоения факторов возможны. Они предъявляют требования к дизайну, которые препятствуют полной рандомизации.

Они широко используются в клинических испытаниях, где, когда один фактор легко изменить, а другому требуется гораздо больше времени для изменения. Если экспериментатор должен последовательно выполнить все прогоны для каждого уровня коэффициента, который трудно изменить, получается разделенный график, а коэффициент, который трудно изменить, представляет весь коэффициент графика.

Вот пример: в полевых испытаниях на сельском хозяйстве цель заключалась в определении воздействия двух сортов сельскохозяйственных культур и четырех различных методов орошения. Было доступно восемь полей, но только один тип орошения может быть применен к каждому полю. Поля могут быть разделены на две части с различным разнообразием в каждой части. Весь фактор участка - это орошение, которое должно быть случайным образом распределено по полям. В каждом поле присваивается сорт.

Вот как вы делаете это в R:

install.packages("faraway")
data(irrigation)
summary(irrigation)

library(lme4)

R> (lmer(yield ~ irrigation * variety + (1|field), data = irrigation))
Linear mixed model fit by REML 
Formula: yield ~ irrigation * variety + (1 | field) 
   Data: irrigation 
  AIC  BIC logLik deviance REMLdev
 65.4 73.1  -22.7     68.6    45.4
Random effects:
 Groups   Name        Variance Std.Dev.
 field    (Intercept) 16.20    4.02    
 Residual              2.11    1.45    
Number of obs: 16, groups: field, 8

Fixed effects:
                       Estimate Std. Error t value
(Intercept)               38.50       3.02   12.73
irrigationi2               1.20       4.28    0.28
irrigationi3               0.70       4.28    0.16
irrigationi4               3.50       4.28    0.82
varietyv2                  0.60       1.45    0.41
irrigationi2:varietyv2    -0.40       2.05   -0.19
irrigationi3:varietyv2    -0.20       2.05   -0.10
irrigationi4:varietyv2     1.20       2.05    0.58

Correlation of Fixed Effects:
            (Intr) irrgt2 irrgt3 irrgt4 vrtyv2 irr2:2 irr3:2
irrigation2 -0.707                                          
irrigation3 -0.707  0.500                                   
irrigation4 -0.707  0.500  0.500                            
varietyv2   -0.240  0.170  0.170  0.170                     
irrgtn2:vr2  0.170 -0.240 -0.120 -0.120 -0.707              
irrgtn3:vr2  0.170 -0.120 -0.240 -0.120 -0.707  0.500       
irrgtn4:vr2  0.170 -0.120 -0.120 -0.240 -0.707  0.500  0.500

По сути, эта модель говорит о том, что ирригация и разнообразие являются фиксированными эффектами, а разнообразие вложено в ирригацию. Поля являются случайными эффектами и наглядно будут что-то вроде

I_1 | I_2 | I_3 | I_4

V_1 V_2 | V_1 V_2 | V_1 V_2 | V_1 V_2

Но это был особый вариант с фиксированным эффектом целого сюжета и эффектом сюжета. Могут быть варианты, в которых один или несколько случайны. Там могут быть более сложные проекты, такие как сплит-сплит. В принципе, вы можете пойти дикий и сумасшедший. Но, учитывая базовую структуру и распределение (то есть фиксированное или случайное, вложенное или скрещенное, ..), ясно понятное, lmer-Ninjaу моделирования не возникнет проблем. Может быть, интерпретация будет беспорядок.

Что касается сравнений, скажем, у вас есть lmer1и lmer2:

anova(lmer1, lmer2)

даст вам соответствующий тест, основанный на статистике теста chi-sq со степенями свободы, равными разности параметров.

cf: Faraway, J., Расширение линейных моделей с помощью R.

Казелла Г., Статистический дизайн