Surowa czy ortogonalna regresja wielomianowa?

Chcę regresować zmienną na . Czy powinienem to zrobić przy użyciu surowych czy ortogonalnych wielomianów? Spojrzałem na pytanie na stronie, które się nimi zajmują, ale tak naprawdę nie rozumiem, jaka jest różnica między ich używaniem. $y$$x,x^2,\ldots,x^5$

Dlaczego nie mogę po prostu wykonać „normalnej” regresji, aby uzyskać współczynniki dla$\beta_i$$y=\sum_{i=0}^5 \beta_i x^i$ (wraz z wartościami p i wszystkimi innymi fajnymi rzeczami) i zamiast tego musisz się martwić, czy używasz wielomianów surowych czy ortogonalnych? Ten wybór wydaje mi się poza zakresem tego, co chcę zrobić.

W książce statystyk, którą obecnie czytam (ISLR Tibshirani i in.), Te rzeczy nie zostały wspomniane. W rzeczywistości byli w pewien sposób lekceważeni.
Powodem jest, AFAIK, że w lm()funkcji w R, używając y ~ poly(x, 2)kwot do używania ortogonalnych wielomianów i używając y ~ x + I(x^2)ilości do używania surowych. Ale na str. 116 autorzy twierdzą, że korzystamy z pierwszej opcji, ponieważ ta druga jest „uciążliwa”, co nie pozostawia żadnych wskazówek, że te polecenia faktycznie mają zupełnie inne rzeczy (i w konsekwencji mają różne wyniki).
(trzecie pytanie) Dlaczego autorzy ISLR tak mylą swoich czytelników?

Uważam, że odpowiedź nie dotyczy stabilności numerycznej (choć odgrywa to pewną rolę), a bardziej zmniejszenia korelacji.

W gruncie rzeczy - kwestia sprowadza się do tego, że kiedy cofamy się w stosunku do wielomianów wysokiego rzędu, zmienne towarzyszące, w których się cofamy, stają się wysoce skorelowane. Przykładowy kod poniżej:

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

To jest niezwykle ważne. W miarę, jak zmienne towarzyszące stają się bardziej skorelowane, nasza zdolność do określania, które są ważne (i jaka jest ich wielkość) gwałtownie maleje. Jest to zwykle określane jako problem wielokoliniowości. Na granicy, gdybyśmy mieli dwie zmienne, które były w pełni skorelowane, kiedy regresujemy je przeciwko czemuś, niemożliwe jest ich rozróżnienie - można to potraktować jako skrajną wersję problemu, ale problem ten wpływa na nasze szacunki dla mniejszy stopień korelacji. Zatem w prawdziwym sensie - nawet jeśli niestabilność numeryczna nie stanowiła problemu - korelacja z wielomianów wyższego rzędu wyrządza ogromne szkody naszym procedurom wnioskowania. Przejawi się to jako większe standardowe błędy (a tym samym mniejsze statystyki t), które w przeciwnym razie byście zobaczyli (patrz przykładowa regresja poniżej).

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)

Jeśli uruchomisz ten kod, interpretacja będzie trudna, ponieważ wszystkie współczynniki się zmieniają, więc trudno jest porównać. Patrząc jednak na statystyki T, widzimy, że zdolność do wyznaczania współczynników była DUŻO większa w przypadku wielomianów ortogonalnych. Dla 3 odpowiednich współczynników otrzymałem t-statystyki (560 21449) dla modelu ortogonalnego i tylko (28, -38, 121) dla surowego modelu wielomianowego. To ogromna różnica dla prostego modelu, który miał tylko kilka względnie niskich terminów wielomianowych, które miały znaczenie.

Nie oznacza to, że przychodzi to bez kosztów. Należy pamiętać o dwóch podstawowych kosztach. 1) tracimy pewną interpretowalność w przypadku wielomianów ortogonalnych. Możemy zrozumieć, co x**3oznacza współczynnik na , ale interpretacja współczynnika na x**3-3x(trzeci poli-pustelnik - niekoniecznie to, czego użyjesz) może być znacznie trudniejsza. Po drugie - gdy mówimy, że są to wielomiany, które są ortogonalne - mamy na myśli, że są one ortogonalne w odniesieniu do pewnej miary odległości. Wybranie miary odległości odpowiedniej dla danej sytuacji może być trudne. Jednak powiedziawszy to, uważam, że polyfunkcja ma na celu takie dobranie, aby była ortogonalna w odniesieniu do kowariancji - co jest przydatne w regresjach liniowych.

Dlaczego nie mogę po prostu wykonać „normalnej” regresji, aby uzyskać współczynniki?

$0.4$$0.4000000059604644775390625$

Użycie surowego wielomianu spowoduje problem, ponieważ będziemy mieli ogromną liczbę. Oto mały dowód: porównujemy liczbę warunków macierzy z surowym i ortogonalnym wielomianem.

> 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

Możesz również sprawdzić moją odpowiedź tutaj na przykład.

Dlaczego istnieją duże współczynniki dla wielomianu wyższego rzędu

Wydaje mi się, że kilka z tych odpowiedzi całkowicie mija się z celem. Odpowiedź Haitao rozwiązuje problemy obliczeniowe związane z dopasowaniem surowych wielomianów, ale jasne jest, że OP pyta o różnice statystyczne między tymi dwoma podejściami. Oznacza to, że gdybyśmy mieli idealny komputer, który mógłby dokładnie reprezentować wszystkie wartości, dlaczego wolelibyśmy jedno podejście od drugiego?

$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

^{Utworzono 25.10.2019 przez pakiet reprezentx (v0.3.0)}

Efekt krańcowy Petal.Widthprzy 0 z dopasowania ortogonalnego i jego błąd standardowy są dokładnie równe efektom z surowego dopasowania wielomianowego. Korzystanie z wielomianów ortogonalnych nie poprawia dokładności oszacowań tej samej wielkości między dwoma modelami.

$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

^{Utworzono 25.10.2019 przez pakiet reprezentx (v0.3.0)}

$0.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.927$. Z ortogonalnego modelu wielomianu, ale nie z surowego modelu wielomianu, wiemy, że większość wariancji wyjaśnionej w wyniku jest spowodowana składnikiem liniowym, przy czym bardzo niewiele pochodzi od wyrażenia kwadratowego, a jeszcze mniej od wyrażenia sześciennego. Surowe wartości wielomianowe nie opowiadają tej historii.

Teraz, niezależnie od tego, czy chcesz uzyskać tę korzyść interpretacyjną w porównaniu z korzyścią interpetacyjną polegającą na faktycznym zrozumieniu współczynników modelu, powinieneś użyć wielomianów ortogonalnych. Jeśli wolisz spojrzeć na współczynniki i dokładnie wiedzieć, co one oznaczają (chociaż wątpię, że jeden zwykle tak robi), powinieneś użyć surowych wielomianów. Jeśli cię to nie obchodzi (tzn. Chcesz kontrolować tylko pomieszanie lub generować prognozowane wartości), to naprawdę nie ma znaczenia; obie formy zawierają te same informacje w odniesieniu do tych celów. Argumentowałbym również, że ortogonalne wielomiany powinny być preferowane w regularyzacji (np. Lasso), ponieważ usunięcie terminów wyższego rzędu nie wpływa na współczynniki terminów niższego rzędu, co nie jest prawdą w przypadku wielomianów surowych,

Potwierdzam doskonałą odpowiedź od @ user5957401 i dodaję komentarze na temat interpolacji, ekstrapolacji i raportowania.

Nawet w dziedzinie stabilnych wartości parametrów współczynniki / parametry modelowane przez wielomiany ortogonalne będą miały znacznie mniejsze błędy standardowe niż współczynniki / parametry modelowane przez parametry surowe. Zasadniczo ortogonalne wielomiany są wolnym zestawem deskryptorów zerow kowariancji. To PCA za darmo!

Jedyną potencjalną wadą jest konieczność wyjaśnienia tego komuś, kto nie rozumie zalet deskryptorów zerowej kowariancji. Współczynniki nie są natychmiast interpretowalne w kontekście efektów pierwszego rzędu (podobny do prędkości) lub drugiego rzędu (podobny do przyspieszenia). Może to być dość potworne w otoczeniu biznesowym.

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

Byłbym więc „rzędami wielkości” bardziej pewny, zgłaszając model ortogonalny niż surowy. W praktyce interpolowałbym z dowolnym modelem, ale ekstrapolowałbym tylko z modelem ortogonalnym.

Chciałbym właśnie to skomentować, aby o tym wspomnieć, ale nie mam wystarczającej liczby przedstawicieli, więc spróbuję rozwinąć się w odpowiedź. Być może zainteresuje Cię to, że w części laboratoryjnej 7.8.1 w „Wprowadzenie do uczenia statystycznego” (James i in., 2017, poprawione 8. drukowanie), omawiają pewne różnice między używaniem ortogonalnych wielomianów lub nie, a mianowicie raw=TRUElub raw=FALSEw poly()funkcji. Na przykład szacunki współczynników zmienią się, ale dopasowane wartości nie:

# 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    

Książka omawia również, w jaki sposób przy stosowaniu wielomianów ortogonalnych, wartości p uzyskane za pomocą anova()zagnieżdżonego testu F (w celu zbadania, w jakim stopniu wielomian może być uzasadniony) są takie same, jak te uzyskane przy zastosowaniu standardowego testu t, wyprowadzone przez summary(fit). To pokazuje, że statystyka F jest równa kwadratowi statystyki t w niektórych sytuacjach.