Сырая или ортогональная полиномиальная регрессия?

Я хочу регрессировать переменную $y$ на $x,x^2,\ldots,x^5$ . Должен ли я сделать это, используя сырые или ортогональные полиномы? Я посмотрел на вопрос на сайте, который касается этих вопросов, но я не совсем понимаю, в чем разница между их использованием.

Почему я не могу просто сделать «нормальный» регрессии , чтобы получить коэффициенты $\beta_i$ в $y=\sum_{i=0}^5 \beta_i x^i$ (вместе с р-значения и все другие хороший материал) , а вместо этого придется беспокоиться ли использование необработанные или ортогональные полиномы? Мне кажется, что этот выбор выходит за рамки того, что я хочу сделать.

В статистической книге, которую я сейчас читаю (ISLR Tibshirani et al), эти вещи не упоминались. На самом деле, они были преуменьшены в пути.
Причина в том, AFAIK, что в lm()функции в R использование y ~ poly(x, 2)сумм равно использованию ортогональных многочленов, а использование y ~ x + I(x^2)сумм - использованию необработанных. Но на стр. 116 авторы говорят, что мы используем первый вариант, потому что последний «громоздок», что не оставляет указаний на то, что эти команды на самом деле совершенно разные вещи (и, как следствие, имеют разные выходы).
(третий вопрос) Почему авторы ISLR так путают своих читателей?

Я считаю, что ответ не столько о числовой стабильности (хотя это играет роль), сколько о снижении корреляции.

По сути, проблема сводится к тому, что когда мы регрессируем против группы многочленов высокого порядка, ковариаты, против которых мы регрессируем, становятся сильно коррелированными. Пример кода ниже:

x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
cor(raw.poly)
cor(orthogonal.poly)

Это чрезвычайно важно. По мере того, как ковариаты становятся более коррелированными, наша способность определять, какие из них важны (и каков их размер), быстро разрушается. Обычно это называют проблемой мультиколлинеарности. На пределе, если у нас было две переменные, которые были полностью коррелированы, когда мы регрессировали их против чего-либо, невозможно различить эти две - вы можете думать об этом как об экстремальной версии проблемы, но эта проблема влияет на наши оценки для меньшие степени корреляции также. Таким образом, в реальном смысле - даже если численная нестабильность не была проблемой - корреляция из полиномов более высокого порядка наносит огромный ущерб нашим процедурам вывода. Это будет проявляться в виде более крупных стандартных ошибок (и, следовательно, меньших значений t-статистики), которые вы могли бы увидеть (см. Пример регрессии ниже).

y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
summary(raw.mod)
summary(orthogonal.mod)

Если вы запустите этот код, интерпретация будет трудной задачей, потому что все коэффициенты меняются, и поэтому трудно сравнивать. Глядя на Т-статистику, мы видим, что способность определять коэффициенты была НАМНОГО больше с ортогональными полиномами. Для 3 соответствующих коэффициентов я получил t-stats (560,21,449) для ортогональной модели и только (28, -38,121) для необработанной полиномиальной модели. Это огромная разница для простой модели с несколькими полиномиальными членами относительно низкого порядка, которые имели значение.

Это не означает, что это происходит без затрат. Необходимо учитывать две основные затраты. 1) мы теряем некоторую интерпретируемость с ортогональными полиномами. Мы могли бы понять, что x**3означает коэффициент на , но интерпретировать коэффициент на x**3-3x(третье поле Эрмита - не обязательно то, что вы будете использовать) может быть гораздо сложнее. Второе - когда мы говорим, что эти многочлены ортогональны, мы имеем в виду, что они ортогональны относительно некоторой меры расстояния. Выбор меры расстояния, соответствующей вашей ситуации, может быть затруднен. Тем не менее, сказав это, я считаю, что polyфункция разработана таким образом, чтобы она была ортогональной по отношению к ковариации - что полезно для линейных регрессий.

Почему я не могу просто сделать "нормальную" регрессию, чтобы получить коэффициенты?

$0.4$$0.4000000059604644775390625$

Использование сырого полинома вызовет проблемы, потому что у нас будет огромное количество. Вот небольшое доказательство: мы сравниваем число условий матрицы с необработанным и ортогональным многочленом.

> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13

Вы также можете проверить мой ответ здесь для примера.

Почему существуют большие коэффициенты для полинома высшего порядка?

Я чувствую, что некоторые из этих ответов полностью упускают суть. Ответ Хайтао решает вычислительные проблемы с подгонкой необработанных полиномов, но ясно, что OP задает вопрос о статистических различиях между двумя подходами. То есть, если бы у нас был идеальный компьютер, который мог бы точно представлять все значения, почему бы мы предпочли один подход другому?

$R^2$$X$$Y$$X=0$$X=0$$X$

data("iris")

#Raw:
fit.raw <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
summary(fit.raw)

#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        1.1034     0.1304   8.464 2.50e-14 ***
#> Petal.Width        1.1527     0.5836   1.975  0.05013 .  
#> I(Petal.Width^2)   1.7100     0.5487   3.116  0.00221 ** 
#> I(Petal.Width^3)  -0.5788     0.1408  -4.110 6.57e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3898 on 146 degrees of freedom
#> Multiple R-squared:  0.9522, Adjusted R-squared:  0.9512 
#> F-statistic: 969.9 on 3 and 146 DF,  p-value: < 2.2e-16

#Orthogonal
fit.orth <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), data = iris)

#Marginal effect of X at X=0 from orthogonal model
library(margins)
summary(margins(fit.orth, variables = "Petal.Width", 
                at = data.frame(Petal.Width = 0)))
#> Warning in check_values(data, at): A 'at' value for 'Petal.Width' is
#> outside observed data range (0.1,2.5)!
#>       factor Petal.Width    AME     SE      z      p  lower  upper
#>  Petal.Width      0.0000 1.1527 0.5836 1.9752 0.0482 0.0089 2.2965

^{Создано в 2019-10-25 пакетом представлением (v0.3.0)}

Предельный эффект Petal.Widthпри 0 от ортогонального соответствия и его стандартная ошибка в точности равны таковым из необработанного полиномиального соответствия. Использование ортогональных полиномов не повышает точность оценок одной и той же величины между двумя моделями.

$Y$$X$$Y$$X$

library(jtools)
data("iris")

fit.raw3 <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
fit.raw1 <- lm(Petal.Length ~ Petal.Width, data = iris)

round(summ(fit.raw3, part.corr = T)$coef, 3)
#>                    Est.  S.E. t val.     p partial.r part.r
#> (Intercept)       1.103 0.130  8.464 0.000        NA     NA
#> Petal.Width       1.153 0.584  1.975 0.050     0.161  0.036
#> I(Petal.Width^2)  1.710 0.549  3.116 0.002     0.250  0.056
#> I(Petal.Width^3) -0.579 0.141 -4.110 0.000    -0.322 -0.074

round(summ(fit.raw1, part.corr = T)$coef, 3)
#>              Est.  S.E. t val. p partial.r part.r
#> (Intercept) 1.084 0.073 14.850 0        NA     NA
#> Petal.Width 2.230 0.051 43.387 0     0.963  0.963

fit.orth3 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), 
               data = iris)
fit.orth1 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3)[,1], 
               data = iris)

round(summ(fit.orth3, part.corr = T)$coef, 3)
#>                                Est.  S.E.  t val. p partial.r part.r
#> (Intercept)                   3.758 0.032 118.071 0        NA     NA
#> stats::poly(Petal.Width, 3)1 20.748 0.390  53.225 0     0.975  0.963
#> stats::poly(Petal.Width, 3)2 -3.015 0.390  -7.735 0    -0.539 -0.140
#> stats::poly(Petal.Width, 3)3 -1.602 0.390  -4.110 0    -0.322 -0.074

round(summ(fit.orth1, part.corr = T)$coef, 3)
#>                                    Est.  S.E. t val. p partial.r part.r
#> (Intercept)                       3.758 0.039 96.247 0        NA     NA
#> stats::poly(Petal.Width, 3)[, 1] 20.748 0.478 43.387 0     0.963  0.963

^{Создано в 2019-10-25 пакетом представлением (v0.3.0)}

$0.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.927$, Из ортогональной полиномиальной модели, но не исходной полиномиальной модели, мы знаем, что большая часть дисперсии, объясненной в результате, обусловлена ​​линейным членом, с очень небольшим числом, полученным от квадратного члена и еще меньше от кубического члена. Необработанные полиномиальные значения не рассказывают эту историю.

Теперь, если вы хотите получить эту интерпретационную выгоду по сравнению с интерпетационной выгодой от фактической способности понимать коэффициенты модели, тогда вам следует использовать ортогональные полиномы. Если вы предпочитаете смотреть на коэффициенты и точно знать, что они означают (хотя я сомневаюсь, что один из них обычно это делает), то вам следует использовать необработанные полиномы. Если вам все равно (то есть вы хотите контролировать только смешанные или генерировать прогнозные значения), то это действительно не имеет значения; Обе формы несут одинаковую информацию об этих целях. Я также утверждал бы, что ортогональные многочлены должны быть предпочтительнее при регуляризации (например, лассо), потому что удаление членов более высокого порядка не влияет на коэффициенты членов более низкого порядка, что не верно для необработанных многочленов,

Я подтверждаю отличный ответ от @ user5957401 и добавляю комментарии по интерполяции, экстраполяции и отчетности.

Даже в области стабильных значений параметров коэффициенты / параметры, моделируемые ортогональными полиномами, будут иметь существенно меньшие стандартные ошибки, чем коэффициенты / параметры, моделируемые необработанными параметрами. По существу, ортогональные полиномы являются свободным набором дескрипторов нулевой ковариации. Это PCA бесплатно!

Единственный потенциальный недостаток - это объяснение тому, кто не понимает достоинства дескрипторов нулевой ковариации. Коэффициенты не могут быть немедленно интерпретированы в контексте эффектов первого порядка (подобного скорости) или второго порядка (подобного ускорению). Это может быть довольно обидно в деловой обстановке.

$10^{-d}$$\sim$$R^2$$adj ~R^2$

Так что я бы был на «порядок» более уверенно сообщать об ортогональной модели, чем необработанную. На практике я буду интерполировать любую модель, но я буду экстраполировать только ортогональную.

Я бы просто прокомментировал, чтобы упомянуть об этом, но мне не хватает представителя, поэтому я постараюсь расширить ответ. Возможно, вам будет интересно увидеть, что в Разделе 7.8.1 лабораторной работы «Введение в статистическое обучение» (James et. Al., 2017, исправлено 8-е издание) они обсуждают некоторые различия между использованием ортогональных полиномов или нет, которое использует raw=TRUEили raw=FALSEв poly()функции. Например, оценки коэффициента будут меняться, но установленные значения не будут:

# using the Wage dataset in the ISLR library
fit1 <- lm(wage ~ poly(age, 4, raw=FALSE), data=Wage)
fit2 <- lm(wage ~ poly(age, 4, raw=TRUE), data=Wage)
print(coef(fit1)) # coefficient estimates differ
print(coef(fit2))
all.equal(predict(fit1), predict(fit2)) #returns TRUE    

В книге также обсуждается, как при использовании ортогональных многочленов p-значения, полученные с помощью anova()вложенного F-критерия (чтобы выяснить, какой степени многочлен может быть оправдан), совпадают с теми, которые получены при использовании стандартного t-критерия, выведенного с помощью summary(fit). Это показывает, что F-статистика в некоторых ситуациях равна квадрату t-статистики.