Odzyskiwanie surowych współczynników i wariancji z ortogonalnej regresji wielomianowej

14

Wydaje się, że jeśli mam model regresji, taki jak , mogę albo dopasować surowy wielomian i uzyskać niewiarygodne wyniki, albo dopasować ortogonalny wielomian i uzyskać współczynniki które nie mają bezpośredniej fizycznej interpretacji (np. nie mogę ich użyć do znalezienia lokalizacji ekstremy w oryginalnej skali). Wydaje się, że powinienem być w stanie mieć to, co najlepsze z obu światów i być w stanie przekształcić dopasowane współczynniki ortogonalne i ich wariancje z powrotem do surowej skali. Ukończyłem kurs dyplomowy z zastosowanej regresji liniowej (używając Kutnera, 5ed) i przejrzałem rozdział dotyczący regresji wielomianowej w Draper (3ed, o którym wspomina Kutner), ale nie znalazłem dyskusji na temat tego, jak to zrobić. Tekst pomocy dlayiβ0+β1xi+β2xi2+β3xi3)poly()funkcja w R nie. Nie znalazłem też niczego w wyszukiwaniu w sieci, w tym tutaj. Rekonstruuje surowe współczynniki (i uzyskuje ich wariancje) ze współczynników dopasowanych do ortogonalnego wielomianu ...

  1. niemożliwe do zrobienia i tracę czas.
  2. być może możliwe, ale nie wiadomo jak w ogólnym przypadku.
  3. możliwe, ale nie omówione, ponieważ „kto by chciał?”
  4. możliwe, ale nie omówione, ponieważ „to oczywiste”.

Jeśli odpowiedź brzmi 3 lub 4, byłbym bardzo wdzięczny, gdyby ktoś miał cierpliwość do wyjaśnienia, jak to zrobić lub wskazał źródło, które to robi. Jeśli jest to 1 lub 2, nadal byłbym ciekawy, co to jest przeszkoda. Dziękuję bardzo za przeczytanie tego i z góry przepraszam, jeśli przeoczyłem coś oczywistego.

f1r3br4nd
źródło
1
Nie rozumiem twoich uwag. x, x 2 i x 3 nie są ortogonalne. Stąd są one skorelowane, a parametry regresji mogą być niestabilne, ale nie jest tak, że są one niewiarygodne. Konwersja na wielomiany ortognonalne może być bardziej niezawodna. Ale co sprawia, że ​​współczynnik pierwotnych potęg x jest bardziej interpretowalny niż współczynniki ortogonalnych wielomianów? Jeśli x jest jedyną zmienną jak w modelu y = a + bx, wówczas ∆y = yi-yi-1 = b∆x i b można interpretować jako zmianę y na jednostkę zmiany x. Ale przy zaangażowanych mocach taka interpretacja zostaje utracona. 23
Michael R. Chernick,
Dla uproszczenia zastosowałem model z tylko x jako zmienną, ale w rzeczywistości porównuję krzywe między grupami leczenia. Tak więc, w zależności od tego, które terminy są znaczące i ich wielkości, mogę je interpretować - na przykład ogólne przesunięcie w górę / w dół lub większe / mniejsze początkowe nachylenie. Ponadto, jak mówi moje pytanie, naturalnym porównaniem krzywych jest lokalizacja maksimów / minimów, co jest łatwiejsze do interpretacji, jeśli jest w oryginalnej skali. Twój głos jest na wybór 3, rozumiem?
f1r3br4nd
Nie, nie wiem, czy jest to jeszcze możliwe. Właśnie zrozumiałem, dlaczego chcesz to zrobić.
Michael R. Chernick,
4
Dobrze pamiętać, że model pasuje wielomianów ortogonalnych będzie miał dokładnie taką samą Fit (czyli tyle samo , te same wartości są zamontowane, itd.), Jak model pasuje do wielomianu surowych warunkach. Jeśli więc chcesz powiązać to z pierwotnymi danymi, możesz spojrzeć na współczynniki dla surowych warunków, ale użyć ortogonalnych wielomianów, aby wnioskować o poszczególnych terminach w sposób, który „uwzględnia” zależność między nimi . R2
Makro
1
Jak się okazuje, splajny sześcienne i splajny B same w sobie należą do klasy i są najlepsze z dwóch światów.
Carl

Odpowiedzi:

6

Tak, to możliwe.

Niech być nie ciągłe części wielomianów ortogonalnych obliczane z x i . (Każdy jest wektorem kolumnowym.) Regresowanie ich względem x i musi zapewniać idealne dopasowanie. Możesz to zrobić za pomocą oprogramowania, nawet jeśli nie dokumentuje ono swoich procedur obliczania ortogonalnych wielomianów. Regresja z j daje współczynniki γ i j, dla którychz1,z2,z3)xixizjγij

zij=γj0+xiγj1+xi2γj2+xi3γj3.

Wynikiem jest macierz Γ, która po prawidłowym pomnożeniu przekształca macierz projektową X = ( 1 ; x ; x 2 ; x 3 ) w Z = ( 1 ; z 1 ; z 2 ; z 3 ) = X Γ .4×4ΓX=(1;x;x2;x3)

(1)Z=(1;z1;z2;z3)=XΓ.

Po dopasowaniu modelu

E(Y)=Zβ

β^(1)

Y^=Zβ^=(XΓ)β^=X(Γβ^).

Γβ^x

Poniższy Rkod ilustruje te procedury i testuje je przy użyciu danych syntetycznych.

n <- 10        # Number of observations
d <- 3         # Degree
#
# Synthesize a regressor, its powers, and orthogonal polynomials thereof.
#
x <- rnorm(n)
x.p <- outer(x, 0:d, `^`); colnames(x.p) <- c("Intercept", paste0("x.", 1:d))
z <- poly(x, d)
#
# Compute the orthogonal polynomials in terms of the powers via OLS.
#
xform <- lm(cbind(1, z) ~ x.p-1)
gamma <- coef(xform)
#
# Verify the transformation: all components should be tiny, certainly
# infinitesimal compared to 1.
#
if (!all.equal(as.vector(1 + crossprod(x.p %*% gamma - cbind(1,z)) - 1), 
    rep(0, (d+1)^2)))
  warning("Transformation is inaccurate.")
#
# Fit the model with orthogonal polynomials.
#
y <- x + rnorm(n)
fit <- lm(y ~ z)
#summary(fit)
#
# As a check, fit the model with raw powers.
#
fit.p <- lm(y ~ .-1, data.frame(x.p))
#summary(fit.p)
#
# Compare the results.
#
(rbind(Computed=as.vector(gamma %*% coef(fit)), Fit=coef(fit.p)))

if (!all.equal(as.vector(gamma %*% coef(fit)), as.vector(coef(fit.p))))
  warning("Results were not the same.")
Whuber
źródło
Γ
110161
Dwa lata później ... @whuber, czy można to również rozszerzyć na 95% CI współczynników?
user2602640,
@ user2602640 Tak. Musisz wyodrębnić macierz wariancji-kowariancji współczynników (użyć vcovw R), aby przekonwertować wariancje obliczone jako jedna podstawa na wariancje w nowej podstawie, a następnie ręcznie obliczyć CI w zwykły sposób.
whuber
@whuber Śledziłem twój komentarz w połowie, a potem straciłem cię całkowicie ... czy jest szansa, że ​​zlitujesz się nad matematykiem, który podważa matematykę, i napiszesz go w kodzie?
user2602640,