Почему lme и aov возвращают разные результаты для повторных измерений ANOVA в R?

24

Я пытаюсь перейти от использования ezпакета lmeдля повторных измерений ANOVA (как я надеюсь, я смогу использовать пользовательские контрасты с lme).

Следуя совету из этого поста в блоге, я смог настроить одну и ту же модель, используя обе функции aov(как это делается ezпри запросе) и lme. Однако, в то время как в примере , приведенном в этой должности в F -значения делать совершенно согласен между aovи lme(я проверил это, и они делают), это не относится к моим данным. Хотя значения F похожи, они не одинаковы.

aovвозвращает значение f 1,3399, lmeвозвращает 1,3664. Я готов принять aovрезультат как «правильный», поскольку это также то, что возвращает SPSS (и это то, что имеет значение для моего поля / руководителя).

Вопросов:

  1. Было бы здорово, если бы кто-то мог объяснить, почему существует такая разница и как я могу использовать ее lmeдля получения достоверных результатов. (Я также хотел бы использовать lmerвместо lmeэтого типа материала, если он дает «правильный» результат. Однако я до сих пор не использовал его.)

  2. После решения этой проблемы я хотел бы провести анализ контраста. Особенно меня будет интересовать контраст объединения первых двух уровней фактора (т. Е. c("MP", "MT")) И сравнение его с третьим уровнем фактора (т. Е. "AC"). Кроме того, тестирование третьего по сравнению с четвертым уровнем фактора (то есть, по "AC"сравнению с "DA").

Данные:

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

И код:

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))
Хенрик
источник
Похоже, вы только что ответили на часть о контрастах в своем ответе здесь ; если нет, пожалуйста, отредактируйте этот вопрос, чтобы мы знали, какие трудности остаются.
Аарон - Восстановить Монику
2
@ Аарон, если есть различия в lmeрезультатах из стандартного учебника ANOVA (предоставленного aovи который мне нужен), это не вариант для меня. В своей статье я хочу сообщить о ANOVA, а не о чем-то вроде ANOVA. Интересно, что Venables & Ripley (2002, стр. 285) показывают, что оба подхода приводят к одинаковым оценкам. Но различия в значениях F оставляют у меня плохое предчувствие. Кроме того, Anova()(from car) возвращает только значения Chi² для lmeобъектов. Поэтому для меня мой первый вопрос еще не ответил.
Хенрик,
Я понимаю (но не разделяю) вашу настороженность lme; но для контрастов, glhtработает на lmподгонки тоже, не только lmeподходит. (Кроме того, lmeрезультаты являются стандартными результатами учебника тоже.)
Аарон - Восстановить Монику
К сожалению, вы не можете указать lmдля анализа повторных измерений. Только aovможет иметь дело с повторными мерами, но вернет объект класса, aovlistкоторый, к сожалению, не обрабатывается glht.
Хенрик
3
lmиспользует остаточную ошибку в качестве члена ошибки для всех эффектов; когда есть эффекты, которые должны использовать другой термин ошибки, aovнеобходимо (или вместо этого, используя результаты lmдля вычисления F-статистики вручную). В вашем примере термин ошибки для factor- это id:factorвзаимодействие, которое является термином остаточной ошибки в аддитивной модели. Сравните свои результаты с anova(lm(value~factor+id)).
Аарон - Восстановить Монику

Ответы:

28

Они отличаются, потому что модель Ime заставляет компонент дисперсии idбыть больше нуля. Глядя на необработанную таблицу anova для всех терминов, мы видим, что среднеквадратичная ошибка для id меньше, чем для остаточных значений.

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131   

Когда мы вычисляем компоненты дисперсии, это означает, что дисперсия из-за id будет отрицательной. Моя память об ожидаемых средних квадратах шаткая, но вычисление что-то вроде

(0.15052-0.16131)/3 = -0.003597.

Это звучит странно, но может случиться. Это означает, что средние значения для каждого идентификатора ближе друг к другу, чем вы ожидаете, учитывая количество остаточных изменений в модели.

Напротив, использование Ime заставляет эту дисперсию быть больше нуля.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

Это сообщает о стандартных отклонениях, возводя в квадрат, чтобы получить выход дисперсии 9.553e-10для идентификатора дисперсии и 0.1586164для остаточной дисперсии.

Теперь вы должны знать, что использование aovдля повторных измерений уместно, только если вы считаете, что корреляция между всеми парами повторных измерений идентична; это называется сложной симметрией. (Технически, сферичности требуется , но этого достаточно для теперь.) Одной из причин для использования lmeболее aov, что он может обрабатывать различные виды корреляционных структур.

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

Использование lmeсо случайным эффектом эквивалентно подгонке модели составной симметрии, где эта корреляция должна быть неотрицательной; мы можем подобрать модель, в которой корреляция может быть отрицательной, используя gls:

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

Этот стол ANOVA согласуется со столом по aovформе и lmформе.

Ок, и что? Что ж, если вы считаете, что отклонение idи корреляция между наблюдениями должны быть неотрицательными, то lmeподгонка на самом деле является более подходящей, чем использование подгонки, aovили lmее оценка остаточной дисперсии немного лучше. Тем не менее, если вы считаете , что корреляция между наблюдениями может быть отрицательным, aovили lmили glsлучше.

Вы также можете быть заинтересованы в дальнейшем изучении корреляционной структуры; чтобы посмотреть на общую структуру корреляции, вы бы сделали что-то вроде

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Здесь я ограничиваю вывод только структурой корреляции. Значения от 1 до 4 представляют четыре уровня factor; мы видим, что фактор 1 и фактор 4 имеют довольно сильную отрицательную корреляцию:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

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

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965
Аарон - Восстановить Монику
источник
2
Фактически можно использовать составную симметрию с, lmeчтобы получить те же результаты, что и с aov(и, следовательно, включить lmeвсе ANOVA), а именно, использовать аргумент корреляции в вызове lme:anova(lme(value ~ factor, data = tau.base, random = ~1|id, correlation = corCompSymm(form = ~1|id)))
Хенрик,
1
Хорошая находка. Но нет ли дополнительного параметра в этом соответствии? Он имеет три параметра дисперсии; дисперсия для id, остаточная дисперсия и корреляция, в то время как gls имеет только остаточную дисперсию и корреляцию.
Аарон - Восстановить Монику
1
Ваш аргумент звучит правдоподобно, однако результаты не согласны. Все таблицы ANOVA ( aov, lmeбез соединения симметрии, и lmeс соединением симметрии) имеют точно такое же количество ДФСА.
Хенрик
1
Вы должны будете убедить меня в том, что эти три параметра действительно являются чрезмерной параметризацией первых двух. Вы выяснили, как они связаны?
Аарон - Восстановить Монику
1
Я доверяю выводу anova.lme(). Из вашего ответа я понял, что связь между ANOVA и смешанными моделями заключается в их корреляционной структуре. Затем я прочитал, что введение симметричной структуры корреляции компоновки приводит к равенству между двумя подходами. Поэтому я навязал это. Я понятия не имею, если это съедает другой дф. Вывод однако не согласуется с этой интерпретацией.
Хенрик
2

aov()соответствует модели с lm()использованием наименьших квадратов, lmeсоответствует максимальной вероятности. Эта разница в том, как оцениваются параметры линейной модели, вероятно, объясняет (очень маленькую) разницу в ваших f-значениях.

На практике (например, для проверки гипотез) эти оценки одинаковы, поэтому я не понимаю, как одно можно считать «более достоверным», чем другое. Они приходят из разных модельных парадигм.

Для контрастов вам нужно установить контрастную матрицу для ваших факторов. Venebles и Ripley показывают, как это сделать, на стр. 143, стр. 146 и стр. 293-294 4-го издания.

Крис
источник
Хм, но почему тогда иногда есть различия, а иногда результаты абсолютно равны? Далее, кажется, что невозможно использовать lmeили lmerдля расчета ANOVA (строго говоря), поскольку он использует метод, который похож, но не идентичен. Так нет ли способа вычисления контрастов для ANOVA с повторным измерением в R?
Хенрик
Если система вашего моделирования действительно линейная, то наименьшие квадраты и ML должны давать одинаковую f-статистику. Только при наличии другой структуры данных эти два метода будут давать разные результаты. Пинейро и Бейтс рассказывают об этом в своей книге о моделях смешанных эффектов. Кроме того, они, вероятно, не «точно» равны, если бы вы пошли достаточно далеко в цифрах, я уверен, что вы найдете некоторую разницу. Но для всех практических целей они одинаковы.
Крис