Как сделать SS-ANOVA типа III в R с контрастными кодами?

26

Пожалуйста, предоставьте R код, который позволяет проводить ANOVA между субъектами с -3, -1, 1, 3 контрастами. Я понимаю, что есть спор относительно соответствующего типа суммы квадратов (SS) для такого анализа. Однако тип SS по умолчанию, используемый в SAS и SPSS (тип III), считается стандартом в моей области. Поэтому я хотел бы, чтобы результаты этого анализа идеально соответствовали тем, что генерируются этими статистическими программами. Чтобы быть принятым ответ должен напрямую вызывать aov (), но за другие ответы может быть проголосовано (особенно, если они просты для понимания / использования).

sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))

Изменить: Обратите внимание, что запрашиваемый мной контраст не является простым линейным или полиномиальным контрастом, а является контрастом, полученным из теоретического предсказания, то есть типа контрастов, обсуждаемых Розенталем и Росновом.

russellpierce
источник
5
Я понимаю, что вам нужна сумма типа III, но эта статья ( stats.ox.ac.uk/pub/MASS3/Exegeses.pdf ) хороша для чтения. Это иллюстрирует некоторые интересные моменты.
Suncoolsu
Что касается вашего вопроса, вас может заинтересовать следующее обсуждение: stats.stackexchange.com/questions/60362/… Выбор между ANOVA типа I, II и III не так прост, как кажется.
phx
Ваш вопрос был назван полезным, поскольку он вызвал несколько заученных ответов, но я также отмечаю, что вы согласились с респондентом, который в основном сказал, что предпосылка вопроса неверна. Я надеюсь, что суммирую позицию StaGuy, сказав, что определенные контрасты были по определению «тип I», а обсуждение других типов стало актуальным только при оценке статистики частичной регрессии, предположительно, наиболее важной, когда «машина управляет» автоматизированными методами.
DWin
@DWin: Я не уверен, что полностью следую за тобой. Можно законно использовать другие типы СС, не позволяя «машине управлять» (по крайней мере, как я понимаю эту фразу). Я мог бы быть немного ржавым здесь, но если память служит, другие типы могут быть релевантными, когда не используется частичная регрессия. Например, SS типа III не разделяет основные эффекты взаимодействия. Различие между типами имеет значение именно там, потому что Тип III не является частичным, тогда как Тип I делает. Проблема, как указано, включала только один контраст, и поэтому различие между типами СС было / остается спорным.
Расселпирс
Насколько я понимаю, обоснование, данное SAS для выбора SSS типа III (и, похоже, именно поэтому люди считают, что тип III предпочтительнее), заключается в том, что он лучше поддерживает процесс обратного и прямого выбора.
DWin

Ответы:

22

Сумма квадратов типа III для ANOVA легко доступна через Anova()функцию из автомобильного пакета.

Контраст кодирование может быть сделано несколькими способами, используя C(), в contr.*семью (как указано @nico), или непосредственно contrasts()функции / аргумент. Это подробно описано в §6.2 (стр. 144-151) « Современной прикладной статистики с S» (Springer, 2002, 4-е изд.). Обратите внимание, что aov()это просто функция-обертка для lm()функции. Интересно, когда кто-то хочет контролировать погрешность модели (как в рамках внутрисубъектного дизайна), но в противном случае они оба дают одинаковые результаты (и независимо от того, как вы подходите к вашей модели, вы все равно можете выводить ANOVA или LM- как резюме с summary.aovили summary.lm).

У меня нет SPSS для сравнения двух выходов, но что-то вроде

> library(car)
> sample.data <- data.frame(IV=factor(rep(1:4,each=20)),
                            DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))
> Anova(lm1 <- lm(DV ~ IV, data=sample.data, 
                  contrasts=list(IV=contr.poly)), type="III")
Anova Table (Type III tests)

Response: DV
            Sum Sq Df F value    Pr(>F)    
(Intercept)  18.08  1  21.815  1.27e-05 ***
IV          567.05  3 228.046 < 2.2e-16 ***
Residuals    62.99 76                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

стоит попробовать в первую очередь.

О факторном кодировании в R против SAS: R рассматривает базовый или опорный уровень как первый уровень в лексикографическом порядке, тогда как SAS рассматривает последний. Таким образом, чтобы получить сопоставимые результаты, либо вы должны использовать contr.SAS()или relevel()ваш фактор R.

хл
источник
1
Я не думаю, что в этом ответе используется контраст -3, -1,1,3, который я указал, и, похоже, он не дает 1-го теста контраста.
Расселпирс
@drknexus Да, ты прав. Писал слишком быстро. Что-то вроде Anova(lm(DV ~ C(IV, c(-3,-1,1,3),1), data=sample.data), type="III")должно быть лучше. Пожалуйста, дайте мне знать, если это хорошо с вами.
ЧЛ
Благодарность! Это выглядит хорошо, я проверю это на SPSS завтра и вернусь к вам.
Russellpierce
1
Кстати, взгляните на пакет ez ( cran.r-project.org/web/packages/ez/index.html ) для упаковки кода Anova ...
Тал Галили,
2
@drknexus: Если бы только была страница с запросом о новых функциях и отправкой сообщений для ez ... github.com/mike-lawrence/ez/issues :)
Майк Лоуренс,
11

Это может выглядеть как самореклама (и я полагаю, что это так). Но я разработал пакет lsmeans для R (доступный на CRAN), который предназначен именно для такой ситуации. Вот как это работает для вашего примера:

> sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))
> sample.aov <- aov(DV ~ factor(IV), data = sample.data)

> library("lsmeans")
> (sample.lsm <- lsmeans(sample.aov, "IV"))
 IV    lsmean        SE df   lower.CL  upper.CL
  1 -3.009669 0.2237448 76 -3.4552957 -2.564043
  2 -3.046072 0.2237448 76 -3.4916980 -2.600445
  3  1.147080 0.2237448 76  0.7014539  1.592707
  4  3.049153 0.2237448 76  2.6035264  3.494779

> contrast(sample.lsm, list(mycon = c(-3,-1,1,3)))
 contrast estimate       SE df t.ratio p.value
 mycon    22.36962 1.000617 76  22.356  <.0001

Вы можете указать дополнительные контрасты в списке, если хотите. В этом примере вы получите те же результаты со встроенным линейным полиномиальным контрастом:

> con <- contrast(sample.lsm, "poly")
> con
 contrast   estimate        SE df t.ratio p.value
 linear    22.369618 1.0006172 76  22.356  <.0001
 quadratic  1.938475 0.4474896 76   4.332  <.0001
 cubic     -6.520633 1.0006172 76  -6.517  <.0001

Чтобы подтвердить это, обратите внимание, что "poly"спецификация направляет его на вызов poly.lsmc, что приводит к следующим результатам:

> poly.lsmc(1:4)
  linear quadratic cubic
1     -3         1    -1
2     -1        -1     3
3      1        -1    -3
4      3         1     1

Если вы хотите выполнить совместную проверку нескольких контрастов, используйте testфункцию с joint = TRUE. Например,

> test(con, joint = TRUE)

Это даст тест типа III. В отличие от этого car::Anova(), он будет делать это правильно, независимо от контрастного кодирования, используемого на этапе подбора модели. Это связано с тем, что тестируемые линейные функции указываются напрямую, а не неявно, посредством сокращения модели. Дополнительная особенность заключается в том, что обнаруживается случай, когда проверяемые контрасты являются линейно зависимыми, и вырабатываются правильные статистические данные и степени свободы.

RVL
источник
10

Вы можете посмотреть этот пост в блоге:

Получение того же результата ANOVA приводит к R, как и в SPSS - трудности с суммами квадратов типа II и типа III

( Спойлер: добавьте options(contrasts=c("contr.sum", "contr.poly"))в начале вашего скрипта)

Nico
источник
Обратите внимание на дополнительную информацию о совместном тестировании в ответе rvl.
rvl
7

Когда вы делаете контрасты, вы делаете конкретную, заявленную линейную комбинацию значений ячеек в контексте соответствующего термина ошибки. Как таковое, понятие «Тип СС» не имеет смысла с контрастами. Каждый контраст по сути является первым эффектом, использующим SS типа I. «Тип СС» имеет отношение к тому, что частично выделено или учитывается другими терминами. Для контрастов ничто не частично или не учитывается. Контраст стоит сам по себе.

StatGuy
источник
Вы абсолютно правы.
Расселпирс
3

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

Более естественный и менее подверженный ошибкам способ получить общие контрасты и возможность описать то, что вы сделали, обеспечивается функцией rmsпакета R. contrast.rmsКонтрасты могут быть очень сложными, но для пользователя они очень просты, потому что они указаны с точки зрения различий в прогностических значениях. Тесты и одновременные контрасты поддерживаются. Это обрабатывает эффекты нелинейной регрессии, эффекты нелинейного взаимодействия, частичные контрасты, все виды вещей.

Фрэнк Харрелл
источник
Это все хорошо и хорошо для вас, как признанного авторитетного человека. У других нет влияния, чтобы не соглашаться с рецензентами. Поскольку интерпретация статистики различается, вы бы просили новых людей придерживаться принципа и неоправданно дорого. Я говорю это как кто-то, кто умер мою долю раз на вершине этого (и подобных) холмов. За изменение IMO в этой области отвечают привратники, то есть редакторы и рецензенты.
Расселпирс
Люди, которые действительно хорошо разбираются в данных, имеют широкий выбор рабочих мест и могут иметь возможность работать в областях, где их навыки и мнения уважаются.
Фрэнк Харрелл
1
... и это то, что я делаю сейчас. Но люди, которые приходят к этому вопросу, часто не принадлежат к этому классу. Так же, как я не был 7 лет назад. Я только защищаю немного сочувствия новичку и его обстоятельствам.
Расселпирс
2

Попробуйте команду Anova в автомобильной библиотеке. Используйте аргумент type = "III", так как по умолчанию используется тип II. Например:

library(car)
mod <- lm(conformity ~ fcategory*partner.status, data=Moore, contrasts=list(fcategory=contr.sum, partner.status=contr.sum))
Anova(mod, type="III")
jebyrnes
источник
3
Я знаю, что Мур находится в автомобильной библиотеке, но когда предоставляются образцы данных, спрашивающему легче понять ваш ответ, если вы используете образцы данных.
russellpierce
0

Также для саморекламы я написал функцию именно для этого: https://github.com/samuelfranssens/type3anova

Установите следующим образом:

library(devtools)
install_github(samuelfranssens/type3anova)
library(type3anova)

sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))

type3anova(lm(DV ~ IV, data = sample.data))

Вам также нужно будет установить carпакет.

sam_f
источник
Как бы вы применили это к контрастной части вопроса?
Расселпирс
1
Извините, я не прочитал вопрос должным образом. Моя функция только упростит проведение III типа Anova. Как и StatGuy выше, я не вижу, где SS вступает в игру при тестировании конкретных контрастов.
sam_f