Surowa czy ortogonalna regresja wielomianowa?

22

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

Dlaczego nie mogę po prostu wykonać „normalnej” regresji, aby uzyskać współczynniki dlaβjay=ja=05βjaxja (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?

l7ll7
źródło
1
@Sycorax Wiem, że polyma to coś wspólnego z wielomianami ortogonalnymi, a ja (x ^ 2) nie (choć nie znam szczegółów) - ale dlaczego autorzy ISLR zalecają metodę, która nie działa ? Wydaje się to bardzo mylące, jeśli oba polecenia wydają się robić to samo, ale tylko jedno jest w porządku.
l7ll7,
1
@gung Przejrzałem dokumentację polyi spędziłem już trochę czasu z tym problemem, ale nie mogę zrozumieć, dlaczego poli (x, 2) i x + I (x ^ 2) mają znaczenie? Czy mógłbyś proszę oświecić mnie tutaj w komentarzach, jeśli pytanie jest nie na temat?
l7ll7,
1
@gung Zrobiłem całkowicie ponownie edytować moje pytanie. Ten wybór surowy / ortogonalny wprawia mnie jeszcze bardziej w zakłopotanie - wcześniej myślałem, że to tylko drobna Rtechnika, której nie rozumiałem, ale teraz wydaje się, że jest to problem z pełną statystyką, który utrudnia mi kodowanie regresji, która nie powinna być tak trudne do zakodowania.
l7ll7
2
@gung To faktycznie pomieszało mnie bardziej niż pomogło. Wcześniej myślałem, że powinienem po prostu iść z wielomianami ortogonalnymi, ponieważ wydawało się to właściwą drogą, ale w tej odpowiedzi używane są surowe wielomiany. O dziwo, wszyscy w sieci krzyczą „RTFM”, ale tak naprawdę nie ma jasnej odpowiedzi, kiedy z czego korzystać. (Twój link również nie daje odpowiedzi na to pytanie, tylko przykład, gdy orth. Pol. Pójdzie źle)
l7ll7
2
O ile nie pracujesz w jakiejś dziedzinie fizycznej lub inżynierskiej, która stwierdza, że ​​odpowiedź będzie kwintycznym wielomianem, prawie na pewno właściwym podejściem nie jest regresja wielomianowa. Zainwestuj swoje stopnie swobody w splajn lub coś, co byłoby znacznie bardziej elastyczne i stabilne niż dopasowanie wielomianowe.
whuber

Odpowiedzi:

10

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.

użytkownik5957401
źródło
3
-1. Większe standardowe błędy, które widać na współczynnikach niższego rzędu, to czerwony śledź. Współczynniki niższego rzędu w twoich dwóch modelach szacują zupełnie różne rzeczy, więc porównywanie ich standardowych błędów nie ma sensu. Współczynnik najwyższego rzędu jest jedynym, który ocenia to samo w obu modelach, a zobaczysz, że statystyka t jest identyczna, niezależnie od tego, czy wielomiany są ortogonalne, czy nie. Twoje dwa modele są statystycznie równoważne pod względem dopasowanych wartości, R ^ 2 itd., Różnią się one głównie tylko interpretacją współczynników
Jake Westfall
@JakeWestfall, nie sądzę, że się z tobą zgadzam. Po pierwsze, uruchomienie kodu generuje wartości, które są różne dla wszystkich rzędów wielomianów, nie dla wszystkich oprócz jednego - w zasadzie bierze wielomian i robi na nim PCA. Po drugie i, co ważniejsze, statystyki t są zasadniczo różne - uruchomienie kodu w mojej odpowiedzi potwierdzi, że - funkcjonalnie rozwiązujemy problem wielokoliniowości. Masz rację, że dopasowane wartości, R ^ 2, testy F itp. Się nie zmieniają. To w rzeczywistości jest powodem ortogonalizacji - nie zmienia to nic poza naszą zdolnością do wykrywania warunków wielomianowych .
user5957401
1
Re: pierwszy punkt, przepraszam, chciałem odnieść się do statystyki t terminu najwyższego rzędu, a nie jego współczynnika. Ten predyktor jest skalowany + przesuwany między modelami, więc tak zmienia się cewka, ale testuje ten sam merytoryczny efekt, jak pokazuje t
Jake Westfall,
Re: drugi punkt, powód „statystyki t są zasadniczo różne” dla terminów niższego rzędu, znowu, ponieważ szacują zupełnie różne rzeczy w dwóch modelach. Rozważmy efekt liniowy: w raw.modnim szacuje nachylenie krzywej przy x = 0, w orthogonal.modnim szacuje krańcowe nachylenie (tj. Identyczne z tym, lm(y ~ poly(x,1))gdzie pominięto terminy wyższego rzędu). Nie ma powodu, aby oszacowania tych całkowicie różnych oszacowań miały porównywalne błędy standardowe. Łatwo można skonstruować kontrofilm, w którym statystyki raw.modsą znacznie wyższe
Jake Westfall,
@JakeWestfall. Nadal uważam, że brakuje ci wielokoliniowości. Wydaje się jednak, że rozmawiamy obok siebie i być może istnieje rozwiązanie. Mówisz, że możesz łatwo zbudować kontrprzykład, proszę. Myślę, że zobaczenie dgp, o którym myślisz, bardzo by mnie wyjaśniło. W tej chwili jedyne rzeczy, które udało mi się wymyślić, które mogą się zachowywać tak, jak to opisujesz, to poważne błędne określenie modelu.
user5957401,
8

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

0.40.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

Haitao Du
źródło
6
Wygląda na to, że używasz pływaków o pojedynczej precyzji i cytujesz je, aby zwiększyć czterokrotnie precyzję! Jak to się stało? Z wyjątkiem GPU, prawie wszystkie obliczenia statystyczne wykorzystują co najmniej podwójną precyzję. Np. RNa wyjściu print(0.4, digits=20)jest 0.40000000000000002.
whuber
6

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?

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

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.

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

Utworzono 25.10.2019 przez pakiet reprezentx (v0.3.0)

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

Noah
źródło
1
Doskonały wkład. Nie mogę odtworzyć wyników marginalnych (funkcja marginesu wyświetla błąd dotyczący poli, gdy próbuję uruchomić pierwszy blok kodu - nie znam pakietu marginesu) - ale są one dokładnie takie, jak się spodziewam. Jako mała sugestia - należy również uwzględnić wyniki analizy marginesów w modelu surowym. Twój argument jest podważony (nieznacznie) przez zmianę wartości p ze streszczenia na funkcje marginesu (nie mniej zmieniając nasze wnioski!) - co wydaje się być spowodowane przez użycie normalnej zamiast rozkładu. Twój punkt regularyzacji jest doskonały.
user5957401
1
Dziękuję za miłe słowa. Myślę, że trzeba zaliczyć stats::w wywołaniu poly()w lm()dla marginsrozpoznać go (co jest głupie). Chciałem skoncentrować mój argument na szacunkach punktowych i standardowych błędach, i wiem, że przedstawiono wiele obcych i rozpraszających informacji, ale mam nadzieję, że tekst ilustruje moje argumenty.
Noah
To nie to. Dokładnie skopiowałem twój kod, a ty go używasz stats::poly(). Błąd mówi 'degree' must be less than number of unique points- co niewiele mi pomaga. Niemniej jednak margin()tworzy kopie sprawdzalnych stwierdzeń, więc nie jest to ważne.
user5957401,
4

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.

10dR2adj R2

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.

Peter Leopold
źródło
1

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.

Shoeburg
źródło
Komentarze nigdy nie powinny być używane jako odpowiedzi bez względu na numery reputacji.
Michael R. Chernick
Jeśli chodzi o twój ostatni punkt, dotyczy to również wielomianów nieortogonalnych. Współczynnik t-test jest równy testowi F porównując model ze współczynnikiem w nim i model bez dla wszystkich współczynników w regresji (wzięty pojedynczo).
Noah,