Полиномиальная регрессия: почему вам нужно определить полиномиальные термины вне lm (), чтобы избежать особенностей?

Если вы попытаетесь запустить полиномиальную регрессию, где x^2 определяется в функции lm(), полиномиальный член будет отброшен из-за особенностей. Однако, если мы определяем полиномиальный член вне lm(), модель подходит правильно.

Похоже, что это должно работать одинаково в обоих направлениях. Почему нам нужно определить полиномиальный член вне функции lm()? 

x <- round(rnorm(100, mean = 0, sd = 10))
y <- round(x*2.5 + rnorm(100))

# Trying to define x^2 in the model, x^2 is dropped
model_wrong <- lm(y ~ x + x^2)

# Define x^2 as its own object
x2 <- x^2
model_right <- lm(y ~ x + x2)
1 2
regression linear-regression

1 ответ:

lm не знает, где начинается и заканчивается термин в Формуле, если вы не скажете ему, обычно обертывая его в функцию. Для произвольных вычислений вы можете обернуть их в I(...), который говорит функции использовать его как есть: 

set.seed(47)
x <- round(rnorm(100, mean = 0, sd = 10))
y <- round(x*2.5 + rnorm(100))

lm(y ~ x + I(x^2))
#> 
#> Call:
#> lm(formula = y ~ x + I(x^2))
#> 
#> Coefficients:
#> (Intercept)            x       I(x^2)  
#>   2.563e-01    2.488e+00   -3.660e-06

Действительно, Вы можете обернуть x^2 в почти любой вызов функции, который вернет вычисленный вектор, который может быть использован в матрице модели. В некоторых случаях cbind может быть очень удобно, хотя c, identity, или даже {...} будет работать. I является специально построенным, хотя.

В качестве альтернативы, вы можете использовать функцию poly, чтобы сделать оба члена для вас, что очень полезно для полиномов более высокой степени. По умолчанию он генерирует ортогональные полиномы, которые заставят коэффициенты выглядеть иначе: 
lm(y ~ poly(x, 2))
#> 
#> Call:
#> lm(formula = y ~ poly(x, 2))
#> 
#> Coefficients:
#> (Intercept)  poly(x, 2)1  poly(x, 2)2  
#>    1.500000   243.485357    -0.004319

Даже если они будут оценивать одно и то же: 

new <- data.frame(x = seq(-1, 1, .5))

predict(lm(y ~ x + I(x^2)), new)
#>          1          2          3          4          5 
#> -2.2317175 -0.9876930  0.2563297  1.5003505  2.7443695

predict(lm(y ~ poly(x, 2)), new)
#>          1          2          3          4          5 
#> -2.2317175 -0.9876930  0.2563297  1.5003505  2.7443695

Если вы действительно хотите, чтобы ваши коэффициенты были одинаковыми, добавьте raw = TRUE: 

lm(y ~ poly(x, 2, raw = TRUE))
#> 
#> Call:
#> lm(formula = y ~ poly(x, 2, raw = TRUE))
#> 
#> Coefficients:
#>             (Intercept)  poly(x, 2, raw = TRUE)1  poly(x, 2, raw = TRUE)2  
#>               2.563e-01                2.488e+00               -3.660e-06
4