Intuicja stojąca za interakcjami produktu tensorowego w GAM (pakiet MGCV w R)

Uogólnione modele addytywne to takie, w których na przykład . funkcje są płynne i należy je oszacować. Zwykle przez karane splajny. MGCV jest pakietem w R, który to robi, a autor (Simon Wood) pisze książkę o swoim pakiecie z przykładami R. Ruppert i in. (2003) napisać znacznie bardziej dostępną książkę o prostszych wersjach tego samego.

$$y = \alpha + f_1(x_1) + f_2(x_2) + e_i$$

Moje pytanie dotyczy interakcji w tego rodzaju modelach. Co jeśli chcę zrobić coś takiego: jeśli byliśmy w krainie OLS (gdzie to tylko beta) , Nie miałbym problemu z interpretacją . Jeśli oceniamy za pomocą splajnowanych splajnów, nie mam również problemu z interpretacją w kontekście addytywnym.

$$y = \alpha + f_1(x_1) + f_2(x_2) + f_3(x_1\times x_2) + e_i$$ $f$$\hat{f}_3$

Ale pakiet MGCV w GAM ma te rzeczy zwane „wygładzaniem produktu tensorowego”. I google „produkt tensorowy”, a moje oczy natychmiast piorunują, próbując przeczytać wyjaśnienia, które znajduję. Albo nie jestem wystarczająco bystry, albo matematyka nie jest dobrze wyjaśniona, albo jedno i drugie.

Zamiast kodowania

normal = gam(y~s(x1)+s(x2)+s(x1*x2))

produkt tensorowy zrobiłby to samo (?)

what = gam(y~te(x1,x2))

kiedy robię

plot(what)

lub

vis.gam(what)

Dostaję naprawdę fajne wyjście. Ale nie mam pojęcia, co się dzieje w czarnej skrzynce te(), ani jak interpretować wyżej wspomnianego fajnego wyjścia. Pewnej nocy miałem koszmar, który prowadziłem na seminarium. Pokazałem wszystkim fajny wykres, zapytali mnie, co to znaczy, i nie wiedziałem. Potem odkryłem, że nie mam na sobie ubrań.

Czy ktokolwiek mógłby pomóc zarówno mnie, jak i potomstwu, dając odrobinę mechaniki i intuicji na temat tego, co dzieje się tutaj pod maską? Idealnie mówiąc trochę o różnicy między normalnym przypadkiem interakcji addytywnej a przypadkiem tensora? Punkty bonusowe za powiedzenie wszystkiego prostym językiem angielskim przed przejściem do matematyki.

Odpowiem na to w trzech krokach: po pierwsze, określmy dokładnie, co rozumiemy przez gładkość jednowymiarową. Następnie opiszemy gładką wielowymiarową (w szczególności gładką dwóch zmiennych). Na koniec postaram się jak najlepiej opisać gładki produkt tensorowy.

1) Płynna jednowymiarowa

Powiedzmy, że mamy pewne dane odpowiedzi które, jak przypuszczamy, są nieznaną funkcją zmiennej predykcyjnej plus pewien błąd . Model byłby:f $y$$f$$x$$ε$

$$y=f(x)+ε$$

Teraz, aby dopasować ten model, musimy zidentyfikować funkcjonalną formę . Sposób, w jaki to robimy, polega na identyfikowaniu funkcji bazowych, które są nakładane w celu reprezentowania funkcji w całości. Bardzo prostym przykładem jest regresja liniowa, w której funkcjami podstawowymi są po prostu i , punkt przecięcia. Stosując rozszerzenie podstawy, mamyf β 2 x β$f$$f$$β_2x$$β_1$

$$y=β_1+β_2x+ε$$

W postaci macierzy mielibyśmy:

$$Y=Xβ+ε$$

Gdzie jest wektorem kolumny n-na-1, jest matrycą modelu n-na-2, jest wektorem kolumny 2 na 1 o współczynnikach modelu, a jest wektorem kolumn n-na-1 błędów . ma dwie kolumny, ponieważ w naszym rozszerzeniu podstawowym istnieją dwa terminy: termin liniowy i punkt przecięcia.X β ε$Y$$X$$β$$ε$$X$

Ta sama zasada dotyczy rozszerzania bazy w MGCV, chociaż funkcje bazowe są znacznie bardziej wyrafinowane. W szczególności indywidualne funkcje bazowe nie muszą być definiowane w pełnej domenie zmiennej niezależnej . Dzieje się tak często w przypadku korzystania z baz opartych na węzłach (patrz „przykład oparty na węzłach”$x$). Model jest następnie reprezentowany jako suma funkcji podstawowych, z których każda jest oceniana przy każdej wartości zmiennej niezależnej. Jednakże, jak wspomniałem, niektóre z tych funkcji bazowych przyjmują wartość zero poza danym interwałem, a zatem nie przyczyniają się do rozszerzenia bazy poza ten interwał. Jako przykład weźmy pod uwagę sześcienną podstawę splajnu, w której każda funkcja bazowa jest symetryczna względem innej wartości (węzła) zmiennej niezależnej - innymi słowy, każda funkcja bazowa wygląda tak samo, ale jest tylko przesunięta wzdłuż osi zmiennej niezależnej (jest to nadmierne uproszczenie, ponieważ każda praktyczna podstawa będzie również obejmować przechwytywanie i termin liniowy, ale mam nadzieję, że masz pomysł).

Mówiąc wprost, podstawowe rozszerzenie wymiaru mogłoby wyglądać następująco:$i-2$

$$y=β_1+β_2x+β_3f_1(x)+β_4f_2(x)+...+β_if_{i-2} (x)+ε$$

gdzie każda funkcja jest być może funkcją sześcienną zmiennej niezależnej .$f$$x$

Równanie macierzowe można nadal wykorzystać do przedstawienia naszego modelu. Jedyna różnica polega na tym, że jest teraz macierzą n-by-i; oznacza to, że ma kolumnę dla każdego terminu w rozwinięciu podstawowym (w tym punkt przecięcia i termin liniowy). Ponieważ proces rozszerzania podstawy pozwolił nam przedstawić model w postaci równania macierzowego, możemy zastosować liniowe najmniejsze kwadraty, aby dopasować model i znaleźć współczynniki . X $Y=Xβ+ε$$X$$β$

Jest to przykład niezenalizowanej regresji, a jedną z głównych zalet MGCV jest oszacowanie gładkości za pomocą macierzy kar i parametru wygładzania. Innymi słowy, zamiast:

$$β=(X^TX)^{-1}X^TY$$

mamy:

$$β=(X^TX+λS)^{-1}X^TY$$

gdzie jest kwadratową macierzą kar -by- , a jest skalarnym parametrem wygładzania. Nie będę tu wchodził w specyfikację macierzy kar, ale wystarczy powiedzieć, że dla dowolnej podstawy rozwinięcie jakiejś niezależnej zmiennej i zdefiniowanie kwadratowej kary „wiggliness” (na przykład kara drugiej pochodnej), jedna można obliczyć matrycy kara .i i λ$S$$i$$i$$λ$$S$

MGCV może wykorzystywać różne metody szacowania optymalnego parametru wygładzania . Nie będę wchodził w ten temat, ponieważ moim celem tutaj było przedstawienie szerokiego przeglądu tego, jak zbudowana jest gładka jednowymiarowa, co, jak sądzę, zrobiłem.$λ$

2) Płynna wielowymiarowa

Powyższe objaśnienie można uogólnić na wiele wymiarów. Wróćmy do naszego modelu, który daje odpowiedź jako funkcję przewidywaczy i . Ograniczenie do dwóch niezależnych zmiennych zapobiegnie zaśmiecaniu wyjaśnienia notacją tajemną. Model jest wtedy:f x $y$$f$$x$$z$

$$y=f(x,z)+ε$$

Teraz intuicyjnie powinno być oczywiste, że będziemy reprezentować z rozszerzeniem podstawy (to jest superpozycją funkcji bazowych), tak jak to zrobiliśmy w przypadku jednowymiarowego przypadku powyżej. Powinno być również oczywiste, że co najmniej jedna, a prawie na pewno wiele innych z tych funkcji bazowych musi być funkcjami zarówno jak i (gdyby tak nie było, to domyślnie byłoby możliwe do rozdzielenia tak, że ). Wizualną ilustrację wielowymiarowej podstawy splajnu można znaleźć tutaj . Pełne dwuwymiarowe rozszerzenie wymiaru może wyglądać mniej więcej tak:f ( x ) x z f f ( x , z ) = f x ( x ) + f z ( z $f(x,z)$$f(x)$$x$$z$$f$$f(x,z)=f_x(x)+f_z(z)$$i-3$

$$y=β_1+β_2x+β_3z+β_4f_1(x,z)+...+β_if_{i-3} (x,z)+ε$$

Myślę, że jest całkiem jasne, że nadal możemy to przedstawić w postaci macierzy za pomocą:

$$Y=Xβ+ε$$

po prostu oceny poszczególnych funkcji bazowych na każdej unikalnej kombinacji i . Rozwiązaniem jest nadal:$x$$z$

$$β=(X^TX)^{-1}X^TY$$

Obliczenie macierzy kar za drugą pochodną jest bardzo takie samo jak w przypadku jednowymiarowym, z tym wyjątkiem, że zamiast całkowania drugiej pochodnej każdej funkcji bazowej w odniesieniu do jednej zmiennej, całkujemy sumę wszystkich drugich pochodnych (w tym częściowych) w odniesieniu do do wszystkich zmiennych niezależnych. Szczegóły powyższego nie są szczególnie ważne: chodzi o to, że nadal możemy konstruować macierz kar i użyć tej samej metody, aby uzyskać optymalną wartość parametru wygładzania , a biorąc pod uwagę ten parametr wygładzania, wektor współczynników jest nadal:$S$$λ$

$$β=(X^TX+λS)^{-1}X^TY$$

Ta dwuwymiarowa gładka ma teraz karę izotropową : oznacza to, że jedna wartość obowiązuje w obu kierunkach. Działa to dobrze, gdy zarówno jak i są w przybliżeniu w tej samej skali, na przykład w aplikacji przestrzennej. Ale co jeśli zastąpimy zmienną przestrzenną zmienną czasową ? Jednostki mogą być znacznie większe lub mniejsze niż jednostki , a to może zepsuć całkowanie naszych drugich pochodnych, ponieważ niektóre z tych pochodnych przyczynią się nieproporcjonalnie do ogólnej integracji (na przykład, jeśli mierzymy w nanosekundach i x z z t t x t x t x $λ$$x$$z$$z$$t$$t$$x$$t$$x$w latach świetlnych całka drugiej pochodnej w odniesieniu do może być znacznie większa niż całka drugiej pochodnej w odniesieniu do , a zatem „zawrotność” wzdłuż kierunku może pozostać w dużej mierze niezakategoryzowana). Slajd 15 „gładkiego zestawu narzędzi”, który podłączyłem, zawiera więcej szczegółów na ten temat.$t$$x$$x$

Warto zauważyć, że nie rozłożyliśmy funkcji bazowych na binarne podstawy i . Wynika z tego, że wygładzenia wielowymiarowe muszą być zbudowane z podstaw obsługujących wiele zmiennych. Wygładzenia produktu Tensor wspierają konstrukcję baz wielowymiarowych z bazowych jednowymiarowych, jak wyjaśnię poniżej.$x$$z$

3) Produkt Tensor wygładza

Wygładzanie produktu Tensor rozwiązuje problem modelowania odpowiedzi na interakcje wielu wejść z różnymi jednostkami. Załóżmy, że mamy odpowiedź która jest funkcją zmiennej przestrzennej i zmiennej czasowej . Nasz model to:f x $y$$f$$x$$t$

$$y=f(x,t)+ε$$

Co chcielibyśmy zrobić, to zbudować dwuwymiarowy podstawę dla zmiennych i . Będzie to o wiele łatwiejsze, jeśli będziemy mogli reprezentować jako:t $x$$t$$f$

$$f(x,t)=f_x(x)f_t(t)$$

W sensie algebraicznym / analitycznym niekoniecznie jest to możliwe. Pamiętajmy jednak, że dyskretyzujemy dziedziny i (wyobraźmy sobie dwuwymiarową „sieć” określoną przez położenie węzłów na osiach i ), tak że „prawdziwa” funkcja jest reprezentowana przez superpozycję funkcji podstawowych . Tak jak zakładaliśmy, że bardzo złożoną funkcję jednowymiarową można aproksymować prostą funkcją sześcienną w określonym przedziale jej dziedziny, możemy założyć, że funkcja nierozdzielalna może być aproksymowana przez iloczyn funkcji prostszych i$x$$t$$x$$t$$f$$f(x,t)$$f_x(x)$$f_t(t)$ w odstępach - pod warunkiem, że nasz wybór podstawowych wymiarów sprawia, że ​​odstępy te są wystarczająco małe!

Nasza ekspansja podstawy, biorąc pod uwagę wymiarową podstawę w oraz wymiarową podstawę w , wyglądałby wtedy następująco:$i$$x$$j$$t$

$$\begin{align} y = &β_{1} + β_{2}x + β_{3}f_{x1}(x)+β_{4}f_{x2}(x)+...+ \\ &β_{i}f_{x(i-3)}(x)+ β_{i+1}t + β_{i+2}tx + β_{i+3}tf_{x1}(x)+β_{i+4}tf_{x2}(x)+...+ \\ &β_{2i}tf_{x(i-3)}(x)+ β_{2i+1}f_{t1}(t) + β_{2i+2}f_{t1}(t)x + β_{2i+3}f_{t1}(t)f_{x1}(x)+β_{i+4}f_{t1}(t)f_{x2}(x){\small +...+} \\ &β_{2i}f_{t1}(t)f_{x(i-3)}(x)+\ldots+ \\ &β_{ij}f_{t(j-3)}(t)f_{x(i-3)}(x) + ε \end{align}$$

Co można interpretować jako iloczyn tensorowy. Wyobrazić, że oceniano każdej funkcji bazowych w i , tym samym konstruowania n-o-i wzór n-o-j macierzy i odpowiednio. Moglibyśmy obliczyć -by- tensorowy produkt z tych dwóch modeli matryce i reorganizacji do kolumn, tak, że każda kolumna reprezentuje unikalną kombinację . Przypomnijmy, że macierze modelu krańcowego miały odpowiednio kolumny i . Wartości te odpowiadają ich odpowiednim wymiarom bazowym. Nasza nowa podstawa z dwiema zmiennymi powinna zatem mieć wymiar$x$$t$$X$$T$$n^2$$ij$ $X \otimes T$$ij$$i$$j$$ij$, a zatem taka sama liczba kolumn w macierzy modelu.

UWAGA: Chciałbym zwrócić uwagę, że ponieważ wyraźnie skonstruowaliśmy funkcje bazowe iloczynu tensorowego, przyjmując iloczyny funkcji bazowych brzeżnych, podstawy iloczynu tensorowego można konstruować z baz brzeżnych dowolnego typu. Nie muszą obsługiwać więcej niż jednej zmiennej, w przeciwieństwie do wielowymiarowej gładkości omówionej powyżej.

W rzeczywistości proces ten powoduje ogólne rozszerzenie podstawy wymiaru ponieważ pełne zwielokrotnienie obejmuje pomnożenie każdej funkcji bazowej przez przecięcie x (więc odejmujemy ), a także mnożenie każdego funkcja podstawowa według t-przecięcia (więc odejmujemy ), ale musimy dodać przecięcie samo w sobie (więc dodajemy 1). Jest to znane jako stosowanie ograniczenia identyfikowalności.$ij-i-j+1$$t$$β_{x1}$$j$$x$$β_{t1}$$i$

Możemy to przedstawić jako:

$$y=β_1+β_2x+β_3t+β_4f_1(x,t)+β_5f_2(x,t)+...+β_{ij-i-j+1}f_{ij-i-j-2}(x,t)+ε$$

Gdzie każda z wielowymiarowych funkcji bazowych jest iloczynem pary marginalnych funkcji bazowych i . Znów jest całkiem jasne, skonstruując tę ​​podstawę, że nadal możemy to przedstawić za pomocą równania macierzowego:$f$$x$$t$

$$Y=Xβ+ε$$

Które (wciąż) ma rozwiązanie:

$$β=(X^TX)^{-1}X^TY$$

Gdzie macierz modelu ma kolumny. Jeśli chodzi o macierze kar i , są one konstruowane osobno dla każdej zmiennej niezależnej w następujący sposób:i j - i - j + 1 J x J $X$$ij-i-j+1$$J_x$$J_t$

$$J_x=β^T I_j \otimes S_x β$$

i,

$$J_t=β^T S_t \otimes I_i β$$

Pozwala to na ogólną karę anizotropową (różną w każdym kierunku) (uwaga: kary za drugą pochodną są sumowane przy każdym węźle na osi i odwrotnie). Parametry wygładzania i można teraz oszacować w podobny sposób, jak pojedynczy parametr wygładzania był dla wygładzeń jednowymiarowych i wielowymiarowych. W rezultacie ogólny kształt gładkości produktu tensorowego jest niezmienny w przypadku przeskalowywania jego zmiennych niezależnych.t λ x λ$x$$t$$λ_x$$λ_t$

Polecam przeczytanie wszystkich winiet na stronie MGCV, a także „ Uogólnione modele addytywne: i wprowadzenie za pomocą R. ” Niech żyje Simon Wood.