Как правильно настроить контрасты в R

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

Для наглядности взгляните на этот пример, и x, и y являются факторами с тремя уровнями.

x <- y <- gl(3, 2)
# [1] 1 1 2 2 3 3
# Levels: 1 2 3

Первый подход создает контрастную матрицу для квадратичного полинома, т. е. с линейным (.L) и квадратичным трендом (.Q). 3 означает: Создайте 3 - 1th полином.

contrasts(x) <- contr.poly(3)
# [1] 1 1 2 2 3 3
# attr(,"contrasts")
#              .L         .Q
# 1 -7.071068e-01  0.4082483
# 2 -7.850462e-17 -0.8164966
# 3  7.071068e-01  0.4082483
# Levels: 1 2 3

Напротив, второй подход дает полином первого порядка (т. е. только линейный тренд). Это связано с аргументом how.many = 1. Следовательно, создается только 1 контраста.

contrasts(y, how.many = 1) <- contr.poly(3)
# [1] 1 1 2 2 3 3
# attr(,"contrasts")
#              .L
# 1 -7.071068e-01
# 2 -7.850462e-17
# 3  7.071068e-01
# Levels: 1 2 3

Если вас интересует только линейный тренд, второй вариант кажется вам более подходящим.

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

Использование полного ("nlevels - 1") набора контрастов создает ортогональный набор контрастов, которые исследуют полный набор (независимых) конфигураций отклика. Сокращение только до одного предотвращает соответствие модели одной конфигурации (в данном случае квадратичной составляющей, которой фактически обладают наши данные.

Чтобы увидеть, как это работает, используйте встроенный набор данных mtcars и проверьте (смешанное) отношение передач к галлонам. Мы предположим, что чем больше передач, тем лучше (по крайней мере, до определенного момента).

df = mtcars # copy the dataset
df$gear = as.ordered(df$gear) # make an ordered factor

Упорядоченные факторы по умолчанию полиномиальные контрасты, но мы установим их здесь, чтобы они были явными:

contrasts(df$gear) <- contr.poly(nlevels(df$gear))

Затем мы можем смоделировать отношения.

m1 = lm(mpg ~ gear, data = df);
summary.lm(m1)
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  20.6733     0.9284  22.267  < 2e-16 ***
# gear.L        3.7288     1.7191   2.169  0.03842 *  
# gear.Q       -4.7275     1.4888  -3.175  0.00353 ** 
# 
# Multiple R-squared:  0.4292,  Adjusted R-squared:  0.3898 
# F-statistic:  10.9 on 2 and 29 DF,  p-value: 0.0002948

Обратите внимание, что у нас есть F(2,29) = 10,9 для общей модели и p = 0,038 для нашего линейного эффекта с расчетным дополнительным расходом 3,7 мили на галлон на передачу.

Теперь давайте запросим только линейный контраст и запустим «тот же самый» анализ.

contrasts(df$gear, how.many = 1) <- contr.poly(nlevels(df$gear))
m1 = lm(mpg ~ gear, data = df)
summary.lm(m1)
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   21.317      1.034  20.612   <2e-16 ***
# gear.L         5.548      1.850   2.999   0.0054 ** 
# Multiple R-squared:  0.2307,  Adjusted R-squared:  0.205 
# F-statistic: 8.995 on 1 and 30 DF,  p-value: 0.005401

Линейный эффект передачи теперь больше (5,5 миль на галлон) и p ‹‹ 0,05 - Победа? За исключением того, что общее соответствие модели теперь значительно хуже: учитываемая дисперсия теперь составляет всего 23% (было 43%)! Почему ясно, если мы построим отношения:

plot(mpg ~ gear, data = df) # view the relationship

Итак, если вас интересует линейный тренд, но вы также ожидаете (или не знаете) дополнительных уровней сложности, вам также следует протестировать эти более высокие полиномы. Квадратичный (или, вообще, тренд до уровня-1).

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

Если кто-то хочет проверить гипотезу о том, что 4 передачи лучше, чем 3, он может ответить этот вопрос :-)