Aov() в R: какая разница bw Error (id) и Error (id/timevar) спецификация?
В чем разница между спецификациями формул aov(depvar~timevar+Error(id))
и aov(depvar~timevar+Error(id/timevar))
? Эти два варианта дают несколько разные результаты.
Тот же вопрос был задан здесь: https://stats.stackexchange.com/info/60108/how-to-write-the-error-term-in-repeated-measures-anova-in-r
Тем не менее, я хотел бы повторить это с более подходящим примером.
Вот пример, который я создал:
var=rep(NA,180)
id=rep(1:20,each=180/20)
group=rep(rep(1:2,each=9),180/(9*2))
time1=rep(rep(1:3,each=3),180/(3*3))
time2=rep(c(8,15,20),180/3)
var[group==1&time1==1&time2==8]=runif(10,105,115)
var[group==2&time1==1&time2==8]=runif(10,105,115)
var[group==1&time1==1&time2==15]=runif(10,95,105)
var[group==2&time1==1&time2==15]=runif(10,95,105)
var[group==1&time1==1&time2==20]=runif(10,85,95)
var[group==2&time1==1&time2==20]=runif(10,85,95)
var[group==1&time1==2&time2==8]=runif(10,95,105)
var[group==2&time1==2&time2==8]=runif(10,95,105)
var[group==1&time1==2&time2==15]=runif(10,85,95)
var[group==2&time1==2&time2==15]=runif(10,75,85)
var[group==1&time1==2&time2==20]=runif(10,75,85)
var[group==2&time1==2&time2==20]=runif(10,65,75)
var[group==1&time1==3&time2==8]=runif(10,95,105)
var[group==2&time1==3&time2==8]=runif(10,95,105)
var[group==1&time1==3&time2==15]=runif(10,85,95)
var[group==2&time1==3&time2==15]=runif(10,75,85)
var[group==1&time1==3&time2==20]=runif(10,75,85)
var[group==2&time1==3&time2==20]=runif(10,65,75)
df=data.frame(id,var,group,time1,time2)
df$id=factor(df$id)
df$group=factor(df$group)
df$time1=factor(df$time1)
df$time2=factor(df$time2)
Выполнение aov() при этом получает несколько разные результаты в зависимости от спецификации term():
Только один раз:
> summary(aov(var~time1+Error(id),data=df))
Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 19 958.4 50.44
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
time1 2 7538 3769 30.41 6.72e-12 ***
Residuals 158 19584 124
> summary(aov(var~time1+Error(id/time1),data=df))
Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 19 958.4 50.44
Error: id:time1
Df Sum Sq Mean Sq F value Pr(>F)
time1 2 7538 3769 211.5 <2e-16 ***
Residuals 38 677 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 120 18907 157.6
Или для обоих сроков (не вводите вывод здесь ради пространства, вы можете проверить его самостоятельно):
summary(aov(var~group*time1*time2+Error(id/(group*time1*time2)),data=df))
summary(aov(var~group*time1*time2+Error(id),data=df))
Почему это происходит? Какой вариант правильный?
Ответы
Ответ 1
Вот сообщение в блоге, которое поможет сломать то, что каждый означает в разделе "Случайные эффекты в классическом ANOVA".
Из блога здесь приведено краткое описание того, что означает "деление" в терминах Error
.
aov(Y ~ Error(A), data=d) # Lone random effect
aov(Y ~ B + Error(A/B), data=d) # A random, B fixed, B nested within A
aov(Y ~ (B*X) + Error(A/(B*X)), data=d) # B and X interact within levels of A
Итак, из вашего вопроса,
aov(depvar~timevar+Error(id/timevar))
означает, что у вас есть случайный эффект от id
, но затем исправить timevar
с timevar
, вложенным в уровни id
в сравнении с
aov(depvar~timevar+Error(id))
который просто принимает id
как случайные эффекты без ограничения на другие переменные.
Источник: http://conjugateprior.org/2013/01/formulae-in-r-anova/
Это также может оказаться полезным, а это некоторый код, проходящий анализ дисперсии, содержащий некоторые рекомендации по изучению ANOVA.
Ответ 2
Разница между aov(depvar~timevar+Error(id))
и aov(depvar~timevar+Error(id/timevar))
заключается в том, включите ли вы timevar
в качестве случайного эффекта.
Обратите внимание, что существует более одного способа включить переменную как случайный эффект. Вы также можете использовать aov(depvar~timevar+Error(id*timevar))
или aov(depvar~timevar+Error(id + timevar))
. Каждое из них означает нечто совершенно иное, но это может сбивать с толку, потому что они часто дают вам аналогичные результаты при применении к одному набору данных из-за ограничений самих данных.
Слэш /
, используемый в aov()
, обозначает nesting. Когда вы используете /
, R автоматически расширяет его до основного эффекта нижней переменной плюс взаимодействие между дном и вершиной. Например, A/B
автоматически расширяется до A + A:B
. Это похоже на то, как A*B
автоматически расширяется до A + B + A:B
, но при вложенности переменная в гнезде никогда не появляется вне своего гнезда (т.е. Сам по себе не может быть основного эффекта B
).
Вы можете увидеть это расширение в вашем выходе:
> summary(aov(var~time1+Error(id / time1)))
Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 52.24 52.24
Error: id:time1
Df Sum Sq Mean Sq
time1 1 4291 4291
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
time1 1 1239 1238.7 10.19 0.00167 **
Residuals 176 21399 121.6
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Термины Error
обозначают случайные эффекты. Обратите внимание, что вы получаете один для основного эффекта id
, потому что он является базой гнезда, а другой для взаимодействия между id
и time1
, потому что time1
вложен в id
(вы также получаете Error
для Within
, который является основным остаточным членом для модели, т.е. случайным эффектом самих отдельных наблюдений).
Итак, какой правильный подход для ваших данных?
Это зависит от 1) того, как ваши данные на самом деле структурированы, и 2) какую модель вы собираетесь запускать. Примечание. Нет окончательного теста, который вы можете выполнить для данных, чтобы определить структуру или правильную модель; это мыслящее упражнение, а не вычислительное.
В приведенных примерах моделей вы получаете результат var
, а затем то, что, как представляется, группирует переменные group
и id
, а затем два временных переменных time1
и time2
. Каждый id
находится только в 1 группе, а не в обеих группах, что указывает на то, что id вложен внутри группы.
> table(group, id)
id
group 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0
2 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9
Я предполагаю, что id
относится к одному участнику, а 9 измерений на time1
и time2
- это тесты внутри предмета для каждого участника (т.е. каждый участник был измерен 9 раз на var
, так что это повторный расчет мер).
Чтобы сделать это конкретным, скажем, var
- оценка по задаче решения задачи, а time1
и time2
- это минуты, участники могут изучить проблему и количество времени, которое они дадут выполните задачу, соответственно. Поскольку time1
и time2
пересекаются, каждый участник завершает задачу 9 раз при каждой комбинации обстоятельств.
> table(time1, time2)
time2
time1 8 15 20
1 20 20 20
2 20 20 20
3 20 20 20
> table(time1, time2, id)
, , id = 1
time2
time1 8 15 20
1 1 1 1
2 1 1 1
3 1 1 1
, , id = 2
time2
time1 8 15 20
1 1 1 1
2 1 1 1
3 1 1 1
(output truncated)
Участники тестируются группами, причем половина участников в группе 1 и другая половина в группе 2. Возможно, исследование проводилось в классных комнатах, а группа 1 - это один класс, а группа 2 - второй класс. Вероятно, групповая идентичность на самом деле не представляет собой переменную, представляющую интерес, но мы не должны оставлять ее вне модели, потому что может быть какая-то неприятная дисперсия, вызванная различиями между группами. Например, возможно, в первом классе было лучшее освещение, дающее всем членам группы 1 и лучший шанс на то, чтобы хорошо забить головоломки, чем члены группы 2.
Оценки, идентификатор и группа должны быть случайными, а time1 и time2 должны быть фиксированными эффектами (обратите внимание, что это может отличаться для тех же данных, если у вас были разные мысли в модели, например, вы можете рассмотреть группу как фиксированную в зависимости от вашего исследовательского вопроса).
Учитывая эту модель, это будет самая полная версия модели, используя aov()
:
aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df)
Здесь вывод:
> summary(aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df))
Error: group
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 771.7 771.7
Error: group:id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 18 243.8 13.55
Error: group:id:time1
Df Sum Sq Mean Sq F value Pr(>F)
time1 2 7141 3571 181.6 <2e-16 ***
Residuals 38 747 20
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: group:id:time2
Df Sum Sq Mean Sq F value Pr(>F)
time2 2 16353 8176 434.6 <2e-16 ***
Residuals 38 715 19
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: group:id:time1:time2
Df Sum Sq Mean Sq F value Pr(>F)
time1:time2 4 214.5 53.63 5.131 0.00103 **
Residuals 76 794.3 10.45
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In aov(var ~ time1 * time2 + Error(group/id/(time1 * time2)), data = df) :
Error() model is singular
(Наряду со ссылками выше, вот несколько дополнительных рекомендаций по случайным и фиксированным эффектам)