3 Фактор вложенных ANOVA в R

Question

Я пытаюсь воспроизвести 3-факторный вложенный ANOVA anlaysis в документе: Underwood, AJ (1993) Механика пространственно-тиражированных программ выборки для обнаружения воздействия на окружающую среду в переменном мире.

данные для примера (из таблицы 3, Underwood, 1993), может быть получено:

dat <- 
structure(list(B = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L 
), .Label = c("A", "B"), class = "factor"), C = structure(c(2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("C", "I"), class = "factor"), 
    Times = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
    4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), 
    Locations = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 
    1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
    3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 
    2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 
    1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
    3L), X = c(59L, 51L, 45L, 46L, 40L, 32L, 39L, 32L, 25L, 51L, 
    44L, 37L, 55L, 47L, 41L, 31L, 38L, 45L, 41L, 47L, 55L, 43L, 
    36L, 29L, 23L, 30L, 37L, 57L, 50L, 43L, 36L, 44L, 51L, 39L, 
    29L, 23L, 38L, 44L, 52L, 31L, 38L, 45L, 42L, 35L, 28L, 52L, 
    44L, 37L, 51L, 43L, 37L, 38L, 31L, 24L, 60L, 52L, 46L, 30L, 
    37L, 44L, 41L, 34L, 27L, 53L, 46L, 39L, 40L, 34L, 26L, 21L, 
    27L, 35L), Times.unique = structure(c(5L, 5L, 5L, 5L, 5L, 
    5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 
    7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
    4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("A_1", "A_2", "A_3", 
    "A_4", "B_1", "B_2", "B_3", "B_4"), class = "factor")), .Names = c("B", 
"C", "Times", "Locations", "Y", "Times.unique"), row.names = c(NA, 
-72L), class = "data.frame") 

dat 

Кадр данных DAT имеет 4 фактора:

B - имеет два уровня "А" и «B» (до v после)

Времена - 8 уровней, 4 в пределах до «B» и 4 внутри «A», закодированные как 1: 4 в каждом. обратите внимание, что переменная Times.unique это то же самое, но с уникальным кодом для каждого времени (до и после)

Locations - имеет три уровня, все измеряется каждый раз, как до, так и после того, как

C - имеет два уровня (C) и (I). примечание: два местоположения являются контрольными, а один - ударом.

Хотя я понимаю, как анализировать такой дизайн с использованием смешанных моделей (lmer), я хотел бы точно воспроизвести его пример, чтобы я мог запускать некоторые симуляции, чтобы сравнить его метод.

В частности, я пытаюсь воспроизвести значения SS, представленные в таблице 4 в столбце «a». Он подходит дизайн, который имеет SS Д.Ф. значения для следующих условиях:

B -> SS = 66,13, DF = 1

Times (B) -> SS = 280,64, ДФ = 6

Местоположения -> SS = 283.86, DF = 2

В х МЕСТОПОЛОЖЕНИЕ -> SS = 29,26, DF = 2

Времена (В) х Locations-> СС = 575,45, DF = 12

Остаточная -> SS = 2420,00, df = 48

Итого -> СС = 6208,34, DF = 71

Я полагаю, термин Таймс (В) представляет Времена вложенными в пределах до/после обработки "B". В этом примере он игнорирует, что местоположения относятся к методам контроля и воздействия и вообще не учитывают фактор C.

Я пробовал все возможные комбинации. Я могу придумать, чтобы воспроизвести эту вложенную anova, используя как уникальную кодировку Times, так и Times, закодированную как 1: 4 внутри B (до и после). Я попытался использовать аргументы% in /,/и Error(), а также Anova из автомобиля, чтобы изменить тип вычисляемой СС. Примеры% в% и/вложенных припадки включают:

aov(Y~B+Locations+Times%in%B+B:Locations+Times%in%B:Locations, data=dat) 
aov(Y~B+Locations+B/Times+B:Locations+B/Times:Locations, data=dat) 

Я, кажется, не в состоянии повторить значения SS Андервуда точно, в частности, для двух членов взаимодействия. Друг позволяет мне подгонять модель в statistix, где значения SS могут быть воспроизведены точно, поэтому можно получить вышеуказанные значения SS для этой модели.

Может ли кто-нибудь помочь мне подогнать эту модель в R? Я хочу встроить его в более крупную симуляцию и действительно нужно иметь возможность запускать модель в R, чтобы точно определить значения Underwood 1993 SS?

Answer 1

Ваша проблема в том, что dat$Locations является целым числом, когда оно должно быть фактором (три уникальных местоположения). Один намек, что ваш ANOVA линия думает Места занимает только 1 Д.Ф., в то время как Underwood отдает 2.

Просто добавьте строку:

dat$Locations = factor(dat$Locations) 

И тогда ваша строка кода воспроизводит Underwood результаты отлично:

aov(Y~B+Locations+B/Times+B:Locations+B/Times:Locations, data=dat) 
#Call: 
# aov(formula = Y ~ B + Locations + B/Times + B:Locations + B/Times:Locations, 
# data = dat) 
# 
#Terms: 
#      B Locations B:Times B:Locations B:Locations:Times 
#Sum of Squares 66.1250 2836.8611 280.6389  29.2500   575.4444 
#Deg. of Freedom   1   2   6   2    12 
#    Residuals 
#Sum of Squares 2420.0000 
#Deg. of Freedom  48