Regressão polinomial bruta ou ortogonal?

Eu quero regredir uma variável em . Devo fazer isso usando polinômios brutos ou ortogonais? Eu olhei para a pergunta no site que lida com isso, mas eu realmente não entendo qual é a diferença entre usá-lo. x , x 2 , … , x $y$$x,x^2,\ldots,x^5$

Por que não posso simplesmente fazer uma regressão "normal" para obter os coeficientes de y = ∑ 5 i = 0 β i x $\beta_i$$y=\sum_{i=0}^5 \beta_i x^i$ (junto com os valores-p e todas as outras coisas legais) e, em vez disso precisa se preocupar se está usando polinômios brutos ou ortogonais? Essa escolha me parece estar fora do escopo do que eu quero fazer.

No livro de estatísticas que estou lendo atualmente (ISLR de Tibshirani et al), essas coisas não foram mencionadas. Na verdade, eles foram subestimados de certa forma.
O motivo é, AFAIK, que na lm()função em R, use y ~ poly(x, 2)quantidades para usar polinômios ortogonais e use y ~ x + I(x^2)quantidades para usar brutos. Mas na p. 116 os autores dizem que usamos a primeira opção porque esta é "complicada", o que não deixa indicação de que esses comandos realmente atinjam coisas completamente diferentes (e tenham resultados diferentes como conseqüência).
(terceira pergunta) Por que os autores da ISLR confundiram seus leitores assim?

Acredito que a resposta seja menos sobre estabilidade numérica (embora isso tenha um papel) e mais sobre como reduzir a correlação.

Em essência - a questão se resume ao fato de que, quando regredimos contra um monte de polinômios de alta ordem, as covariáveis ​​contra as quais estamos regredindo se tornam altamente correlacionadas. Exemplo de código abaixo:

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

Isso é tremendamente importante. À medida que as covariáveis ​​se tornam mais correlacionadas, nossa capacidade de determinar quais são importantes (e qual o tamanho de seus efeitos) diminui rapidamente. Isso geralmente é chamado de problema da multicolinearidade. No limite, se tivéssemos duas variáveis ​​totalmente correlacionadas, quando as regredimos contra algo, é impossível distinguir entre as duas - você pode pensar nisso como uma versão extrema do problema, mas esse problema afeta nossas estimativas para graus menores de correlação também. Assim, em um sentido real - mesmo que a instabilidade numérica não tenha sido um problema - a correlação de polinômios de ordem superior causa um tremendo dano às nossas rotinas de inferência. Isso se manifestará como erros padrão maiores (e, portanto, estatísticas t menores) que você veria de outra forma (veja o exemplo de regressão abaixo).

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)

Se você executar esse código, a interpretação é um pouco difícil, porque todos os coeficientes mudam e, portanto, as coisas são difíceis de comparar. Porém, olhando para as estatísticas T, podemos ver que a capacidade de determinar os coeficientes era MUITO maior com os polinômios ortogonais. Para os três coeficientes relevantes, obtive estatísticas t de (560,21,449) para o modelo ortogonal e apenas (28, -38,121) para o modelo polinomial bruto. Essa é uma grande diferença para um modelo simples, com apenas alguns termos polinomiais de ordem relativamente baixa que são importantes.

Isso não quer dizer que isso ocorra sem custos. Há dois custos principais a serem considerados. 1) perdemos alguma interpretabilidade com polinômios ortogonais. Podemos entender o que x**3significa o coeficiente em , mas interpretar o coeficiente em x**3-3x(o terceiro poli hermita - não necessariamente o que você usará) pode ser muito mais difícil. Segundo - quando dizemos que esses polinômios são ortogonais - queremos dizer que são ortogonais em relação a alguma medida de distância. Escolher uma medida de distância relevante para sua situação pode ser difícil. No entanto, tendo dito isso, acredito que a polyfunção seja projetada para escolher de forma que seja ortogonal em relação à covariância - o que é útil para regressões lineares.

Por que não posso simplesmente fazer uma regressão "normal" para obter os coeficientes?

Porque não é numericamente estável. Lembre-se de que o computador está usando um número fixo de bits para representar um número flutuante. Verifique IEEE754 para obter detalhes, você pode surpreender que mesmo o número simples , o computador precise armazená-lo como . Você pode tentar outros números aqui$0.4$$0.4000000059604644775390625$

Usar polinômio bruto causará problemas, pois teremos um número enorme. Aqui está uma pequena prova: estamos comparando o número da condição da matriz com o polinômio bruto e ortogonal.

> 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

Você também pode verificar minha resposta aqui para um exemplo.

Por que existem grandes coeficientes para polinômios de ordem superior

Sinto que várias dessas respostas estão completamente erradas. A resposta de Haitao aborda os problemas computacionais com o ajuste de polinômios brutos, mas é claro que o OP está perguntando sobre as diferenças estatísticas entre as duas abordagens. Ou seja, se tivéssemos um computador perfeito que pudesse representar todos os valores exatamente, por que preferiríamos uma abordagem à outra?

$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

^{Criado em 2019-10-25 pelo pacote reprex (v0.3.0)}

O efeito marginal de Petal.Widthem 0 do ajuste ortogonal e seu erro padrão são exatamente iguais aos do ajuste polinomial bruto. O uso de polinômios ortogonais não melhora a precisão das estimativas da mesma quantidade entre os dois modelos.

$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

^{Criado em 2019-10-25 pelo pacote reprex (v0.3.0)}

$0.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.927$. Do modelo polinomial ortogonal, mas não do modelo polinomial bruto, sabemos que a maior parte da variância explicada no resultado se deve ao termo linear, com muito pouco advindo do termo quadrado e menos ainda do termo cúbico. Os valores polinomiais brutos não contam essa história.

Agora, se você deseja esse benefício interpretativo sobre o benefício interpetacional de ser capaz de realmente entender os coeficientes do modelo, use polinômios ortogonais. Se você preferir olhar para os coeficientes e saber exatamente o que eles significam (embora eu duvide que um seja tipicamente), use os polinômios brutos. Se você não se importa (ou seja, você deseja apenas controlar a confusão ou gerar valores previstos), isso realmente não importa; os dois formulários carregam as mesmas informações com relação a esses objetivos. Eu também argumentaria que polinômios ortogonais devem ser preferidos na regularização (por exemplo, laço), porque a remoção de termos de ordem superior não afeta os coeficientes dos termos de ordem inferior, o que não é verdade em polinômios brutos,

Confirmo a excelente resposta de @ user5957401 e adiciono comentários sobre interpolação, extrapolação e geração de relatórios.

Mesmo no domínio dos valores estáveis ​​dos parâmetros, os coeficientes / parâmetros modelados pelos polinômios ortogonais terão erros padrão substancialmente menores que os coeficientes / parâmetros modelados pelos parâmetros brutos. Essencialmente, os polinômios ortogonais são um conjunto livre de descritores de covariância zero. Isso é PCA de graça!

A única desvantagem potencial é ter que explicar isso para alguém que não entende a virtude dos descritores de covariância zero. Os coeficientes não são imediatamente interpretáveis ​​no contexto dos efeitos de primeira ordem (tipo velocidade) ou de segunda ordem (tipo aceleração). Isso pode ser bastante prejudicial em um ambiente de negócios.

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

Portanto, eu seria "ordens de grandeza" mais confiante em relatar o modelo ortogonal do que o modelo bruto. Na prática, eu interpolava com qualquer um dos modelos, mas extrapolava apenas com o modelo ortogonal.

Eu teria comentado apenas para mencionar isso, mas não tenho representante suficiente, então tentarei expandir para uma resposta. Você pode estar interessado em ver que, na Seção 7.8.1 do laboratório, em "Introdução à aprendizagem estatística" (James et. Al., 2017, 8ª impressão corrigida), eles discutem algumas diferenças entre o uso de polinômios ortogonais ou não, que está usando o método raw=TRUEou raw=FALSEna poly()função. Por exemplo, as estimativas do coeficiente serão alteradas, mas os valores ajustados não:

# 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    

O livro também discute como quando polinômios ortogonais são usados, os valores p obtidos com o anova()teste F aninhado (para explorar qual grau de polinômio pode ser garantido) são os mesmos que os obtidos ao usar o teste t padrão, produzido por summary(fit). Isso ilustra que a estatística F é igual ao quadrado da estatística t em determinadas situações.