изменить контрасты члена взаимодействия с двоеточием в lm ()

Question

Можно ли изменить contrasts условий взаимодействия, которые были указаны в lm с использованием двоеточия :?

В приведенном ниже примере ссылочной категорией по умолчанию является последний из шести терминов, сгенерированных как gear:vs (то есть gear5:vs1). Вместо этого я хотел бы использовать первую из шести в качестве ссылки (т.е. gear3:vs0).

mtcars.1 <- mtcars %>%
  mutate(gear = as.factor(gear)) %>%
  mutate(vs = as.factor(vs))

lm(data=mtcars.1, mpg ~ gear:vs) %>%
  tidy
#> # A tibble: 6 x 5
#>   term        estimate std.error statistic      p.value
#>   <chr>          <dbl>     <dbl>     <dbl>        <dbl>
#> 1 (Intercept)    30.4       4.13      7.36 0.0000000824
#> 2 gear3:vs0     -15.4       4.30     -3.57 0.00143     
#> 3 gear4:vs0      -9.40      5.06     -1.86 0.0747      
#> 4 gear5:vs0     -11.3       4.62     -2.44 0.0218      
#> 5 gear3:vs1     -10.1       4.77     -2.11 0.0447      
#> 6 gear4:vs1      -5.16      4.33     -1.19 0.245

Задание контрастов для gear и vs по отдельности, похоже, не оказывает влияния:

lm(data=mtcars.1, mpg ~ gear:vs, 
             contrasts = list(gear = contr.treatment(n=3,base=3),
                              vs = contr.treatment(n=2,base=2))) %>%
  tidy
#> # A tibble: 6 x 5
#>   term        estimate std.error statistic      p.value
#>   <chr>          <dbl>     <dbl>     <dbl>        <dbl>
#> 1 (Intercept)    30.4       4.13      7.36 0.0000000824
#> 2 gear3:vs0     -15.4       4.30     -3.57 0.00143     
#> 3 gear4:vs0      -9.40      5.06     -1.86 0.0747      
#> 4 gear5:vs0     -11.3       4.62     -2.44 0.0218      
#> 5 gear3:vs1     -10.1       4.77     -2.11 0.0447      
#> 6 gear4:vs1      -5.16      4.33     -1.19 0.245

И я не уверен, как указать контраст для gear:vs напрямую:

lm(data=mtcars.1, mpg ~ gear:vs,
             contrasts = list("gear:vs" = contr.treatment(n=6,base=6))) %>%
  tidy
#> Warning in model.matrix.default(mt, mf, contrasts): variable 'gear:vs' is
#> absent, its contrast will be ignored
#> # A tibble: 6 x 5
#>   term        estimate std.error statistic      p.value
#>   <chr>          <dbl>     <dbl>     <dbl>        <dbl>
#> 1 (Intercept)    30.4       4.13      7.36 0.0000000824
#> 2 gear3:vs0     -15.4       4.30     -3.57 0.00143     
#> 3 gear4:vs0      -9.40      5.06     -1.86 0.0747      
#> 4 gear5:vs0     -11.3       4.62     -2.44 0.0218      
#> 5 gear3:vs1     -10.1       4.77     -2.11 0.0447      
#> 6 gear4:vs1      -5.16      4.33     -1.19 0.245

^{Создано в 2019-01-21 с помощью представительного пакета (v0.2.1)}

Answer 1

Одним из способов решения этой проблемы является предварительное вычисление члена взаимодействия до регрессии.

Чтобы продемонстрировать, мы можем создать столбец фактора GV в mtcars с теми же уровнями, что и в вашем lm выход.Он генерирует те же значения:

mtcars %>% 
  mutate(GV = interaction(factor(gear), factor(vs)), 
         GV = factor(GV, levels = c("5.1", "3.0", "4.0", "5.0", "3.1", "4.1"))) %>% 
  lm(mpg ~ GV, .) %>% 
  tidy() 

# A tibble: 6 x 5
  term        estimate std.error statistic      p.value
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>
1 (Intercept)    30.4       4.13      7.36 0.0000000824
2 GV3.0         -15.4       4.30     -3.57 0.00143     
3 GV4.0          -9.4       5.06     -1.86 0.0747      
4 GV5.0         -11.3       4.62     -2.44 0.0218      
5 GV3.1         -10.1       4.77     -2.11 0.0447      
6 GV4.1          -5.16      4.33     -1.19 0.245 

Теперь мы опускаем второй член mutate, поэтому уровни составляют 3,0, 4,0, 5,0, 3,1, 4,1, 5,1.

mtcars %>% 
  mutate(GV = interaction(factor(gear), factor(vs))) %>% 
  lm(mpg ~ GV, .) %>% 
  tidy()

  # A tibble: 6 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    15.1       1.19     12.6  1.38e-12
2 GV4.0           5.95      3.16      1.88 7.07e- 2
3 GV5.0           4.08      2.39      1.71 9.96e- 2
4 GV3.1           5.28      2.67      1.98 5.83e- 2
5 GV4.1          10.2       1.77      5.76 4.61e- 6
6 GV5.1          15.4       4.30      3.57 1.43e- 3

Используйте interaction(factor(gear), factor(vs), lex.order = TRUE), чтобы получить уровни 3.0, 3.1, 4.0, 4.1, 5.0, 5.1.

mtcars %>% 
  mutate(GV = interaction(factor(gear), factor(vs), lex.order = TRUE)) %>%
  lm(mpg ~ GV, .) %>% 
  tidy()

# A tibble: 6 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    15.0       1.19     12.6  1.38e-12
2 GV3.1           5.28      2.67      1.98 5.83e- 2
3 GV4.0           5.95      3.16      1.88 7.07e- 2
4 GV4.1          10.2       1.77      5.76 4.61e- 6
5 GV5.0           4.07      2.39      1.71 9.96e- 2
6 GV5.1          15.3       4.30      3.57 1.43e- 3