Определение инфиксного оператора для использования в формуле
Я пытаюсь создать более экономную версию этого решения, которое влечет за собой указание RHS формулы в форме d1 + d1:d2.
Учитывая, что * в контексте формулы является надежным условием полного взаимодействия (т.е. d1 * d2 дает d1 + d2 + d1:d2), мой подход состоял в том, чтобы попытаться определить альтернативный оператор, скажем %+:% используя подход infix, к которому я привык в других приложениях, a la:
"%+:%" <- function(d1,d2) d1 + d2 + d1:d2
Однако это предсказуемо терпит неудачу, потому что я не был осторожен в оценке; давайте представим пример для иллюстрации моего прогресса:
set.seed(1029)
v1 <- runif(1000)
v2 <- runif(1000)
y <- .8*(v1 < .3) + .2 * (v2 > .25 & v2 < .8) -
.4 * (v2 > .8) + .1 * (v1 > .3 & v2 > .8)
В этом примере, надеюсь, ясно, почему простое изложение двух терминов может быть нежелательным:
y ~ cut(v2, breaks = c(0, .25, .8, 1)) +
cut(v2, breaks = c(0, .25, .8, 1)):I(v1 < .3)
Одно решение, близкое к моему желаемому результату, состоит в том, чтобы определить всю формулу как функцию:
plus.times <- function(outvar, d1, d2){
as.formula(paste0(quote(outvar), "~", quote(d1),
"+", quote(d1), ":", quote(d2)))
}
Это дает ожидаемые коэффициенты при передаче в lm, но с именами, которые сложнее интерпретировать напрямую (особенно в реальных данных, где мы позаботимся дать описательные имена d1 и d2, в отличие от этого общий пример):
out1 <- lm(y ~ cut(v2, breaks = c(0, .25, .8, 1)) +
cut(v2, breaks = c(0, .25, .8, 1)):I(v1 < .3))
out2 <- lm(plus.times(y, cut(v2, breaks = c(0, .25, .8, 1)), I(v1 < .3)))
any(out1$coefficients != out2$coefficients)
# [1] FALSE
names(out2$coefficients)
# [1] "(Intercept)" "d1(0.25,0.8]" "d1(0.8,1]" "d1(0,0.25]:d2TRUE"
# [5] "d1(0.25,0.8]:d2TRUE" "d1(0.8,1]:d2TRUE"
Так что это меньше оптимального.
Есть ли способ определить настройку кода так, чтобы оператор инфикса, о котором я говорил выше, работает так, как ожидалось? Как насчет изменения формы plus.times, чтобы переменные не были переименованы?
Я выкалывал (?formula, ?"~", ?":", getAnywhere(formula.default), этот ответ и т.д.), но не видно, как именно R интерпретирует *, когда он встречается в формуле, чтобы я мог внести желаемые незначительные корректировки.
Ответы
Ответ 1
В этом случае вам не нужно определять новый оператор: в формуле d1/d2 расширяется до d1 + d1:d2. Другими словами d1/d2 указывает, что d2 вложен в d1. Продолжая ваш пример:
out3 <- lm(y ~ cut(v2,breaks=c(0,.25,.8,1))/I(v1 < .3))
all.equal(coef(out1), coef(out3))
# [1] TRUE
Дополнительные комментарии
Факторы могут быть пересечены или вложены. Два фактора пересекаются, если можно наблюдать каждую комбинацию уровней двух факторов, например. пола и лечения, температуры и рН и т.д. Фактор вложен в другой, если каждый уровень этого фактора можно наблюдать только на одном из уровней другого фактора, например. города и страны, сотрудника и магазина и т.д.
Эти отношения отражены в параметризации модели. Для скрещенных множителей мы используем d1*d2 или d1 + d2 + d1:d2, чтобы дать основной эффект каждого фактора плюс взаимодействие. Для вложенных факторов мы используем d1/d2 или d1 + d1:d2, чтобы дать отдельную подмодуль формы 1 + d2 для каждого уровня d1.
Идея гнездования не ограничивается факторами, например, мы можем использовать sex/x для установки отдельной линейной регрессии на x для мужчин и женщин.
В формуле %in% эквивалентна :, но ее можно использовать для выделения вложенной или иерархической структуры данных/модели. Например, a + b %in% a совпадает с a + a:b, но чтение его как "плюс b в пределах" дает лучшее описание устанавливаемой модели. Тем не менее использование / имеет преимущество упрощения формулы модели в то же время, что и подчеркивание структуры.
