Regressão polinomial bruta ou ortogonal?

22

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 5yx,x2,,x5

Por que não posso simplesmente fazer uma regressão "normal" para obter os coeficientes de y = 5 i = 0 β i x iβiy=i=05βixi (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?

l7ll7
fonte
1
@ Sycorax Eu sei que polytem algo a ver com polinômios ortogonais e eu (x ^ 2) não (embora eu não conheça os detalhes) - mas ainda assim, por que os autores da ISLR recomendariam um método que não funciona? ? Parece muito enganador se ambos os comandos parecem fazer o mesmo, mas apenas um realmente está ok.
L7ll7
1
@gung Eu olhei para a documentação polye já gastei um tempo com esse problema, mas não consigo descobrir por que poli (x, 2) e x + I (x ^ 2) fazem a diferença? Poderia, por favor, me esclarecer aqui nos comentários, se a pergunta for offtopic?
L7ll7
1
@gung Fiz uma reedição completa da minha pergunta. Essa opção bruta / ortogonal está me confundindo ainda mais - antes eu pensava que isso era apenas um detalhe Rtécnico menor , que eu não entendia, mas agora parece ser um problema estatístico completo que me impede de codificar uma regressão que não deveria ser difícil de codificar.
L7ll7
2
@ gung Isso realmente me confundiu mais do que ajudou. Anteriormente, eu pensava que deveria apenas usar polinômios ortogonais, porque esse parecia ser o caminho certo, mas nessa resposta polinômios brutos são usados. Surpreendentemente, todos na rede estão gritando "RTFM", mas na verdade não há uma resposta clara quando usar o quê. (O seu link também não dá uma resposta para isso, apenas um exemplo, quando Orth pol dar errado..)
l7ll7
2
A menos que você esteja trabalhando em algum domínio físico ou de engenharia que afirme que a resposta será um polinômio quintico, quase certamente a abordagem correta é não fazer regressão polinomial em primeiro lugar. Invista seus graus de liberdade em um spline ou em algo que seria muito mais flexível e estável do que o ajuste polinomial.
whuber

Respostas:

10

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.

user5957401
fonte
3
-1. Os erros padrão maiores que você vê nos coeficientes de ordem inferior são um arenque vermelho. Os coeficientes de ordem inferior em seus dois modelos estão estimando coisas completamente diferentes; portanto, comparar seus erros padrão não faz sentido. O coeficiente de ordem mais alta é o único que estima a mesma coisa nos dois modelos, e você verá que a estatística t é idêntica, independentemente de os polinômios serem ortogonais ou não. Seus dois modelos são estatisticamente equivalentes em termos de valores ajustados, R ^ 2 etc., eles diferem principalmente apenas na interpretação dos coeficientes
Jake Westfall
@JakeWestfall, acho que não concordo com você. Primeiro, executar o código produz valores diferentes para todas as ordens polinomiais, não todas, exceto uma - em essência, ele pega o polinômio e faz PCA nele. Em segundo lugar, e mais importante, as estatísticas t são substancialmente diferentes - executar o código na minha resposta confirmará isso - funcionalmente, estamos resolvendo o problema da multicolinearidade. Você está certo que os valores ajustados, R ^ 2, testes F etc. não mudam. Na verdade, esse é um motivo para ortogonalizar - nada muda, exceto nossa capacidade de detectar termos polinomiais .
User5957401 25/10
1
Re: o primeiro ponto, desculpe, eu pretendia me referir ao estatuto t do termo de maior ordem, não ao seu coeficiente. Esse preditor é escalado + deslocado entre os modelos, então sim, o coef muda, mas testa o mesmo efeito substantivo, como mostra t
Jake Westfall
Re: o segundo ponto, a razão "as estatísticas t são substancialmente diferentes" para os termos de ordem inferior é, novamente, porque eles estão estimando coisas completamente diferentes nos dois modelos. Considere o efeito linear: raw.modcalcula a inclinação da curva em x = 0, orthogonal.modcalcula a inclinação marginal (ou seja, idêntica a lm(y ~ poly(x,1))onde os termos de ordem superior são omitidos). Não há razão para que as estimativas dessas estimativas completamente diferentes tenham erros padrão comparáveis. Pode-se facilmente construir um contra-exemplo onde raw.modhá estatísticas t muito mais altas
Jake Westfall
@JakeWestfall. Eu ainda acho que você está perdendo a multicolinearidade. No entanto, parece que estamos falando um do outro, e talvez haja uma solução. Você diz que pode facilmente construir um contra-exemplo, por favor. Eu acho que ver o dgp que você tem em mente esclareceria muito para mim. No momento, as únicas coisas que pude apresentar que podem se comportar conforme você descreve envolvem uma especificação incorreta do modelo.
User5957401 29/10
8

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 aqui0,40000000596046447753906250.40.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

Haitao Du
fonte
6
Você parece estar usando flutuadores de precisão única e citando-os para quadruplicar a precisão! Como isso aconteceu? Exceto pelas GPUs, quase todo o cálculo estatístico usa pelo menos precisão dupla. Por exemplo, na Rsaída de print(0.4, digits=20)é 0.40000000000000002.
whuber
6

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?

R2XYX=0X=0X

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.

YXYX

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.0010.0030.0050.9270.9270.0200.0050.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,

Noé
fonte
1
Excelente contribuição. Não consigo replicar seus resultados marginais (a função de margem exibe um erro sobre poli quando tento executar seu primeiro bloco de código - não estou familiarizado com o pacote de margens) - mas eles são exatamente o que eu espero. Como uma pequena sugestão - você deve incluir também a saída da análise de margem no modelo bruto. Seu argumento é minado (levemente) pela mudança nos valores de p das funções sumário e margem (alterando nossas conclusões não menos!) - que parece ser causada pelo uso de uma distribuição normal em vez de na distribuição. Seu ponto de regularização é excelente.
user5957401
1
Obrigado pelas amáveis ​​palavras. Eu acho que você tem que incluir o stats::na chamada para poly()a lm()para marginsa reconhecê-lo (o que é estúpido). Eu queria focar meu argumento nas estimativas pontuais e nos erros padrão, e sei que há muitas informações estranhas e perturbadoras apresentadas, mas espero que o texto ilustre meus pontos.
Noah
Não é isso. Copiei seu código exatamente e você usa stats::poly(). O erro diz 'degree' must be less than number of unique points- o que não me ajuda muito. No entanto, margin()está fazendo backup de declarações prováveis, portanto não é importante.
user5957401
4

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.

10dR2adj R2

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.

Peter Leopold
fonte
1

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.

shoeburg
fonte
Os comentários nunca devem ser usados ​​como respostas, independentemente do seu número de reputação.
Michael R. Chernick
Em relação ao seu último ponto, isso também se aplica a polinômios não ortogonais. O teste t do coeficiente é igual ao teste F comparando um modelo com o coeficiente nele e um modelo sem todos os coeficientes em regressão (um de cada vez).
Noah