Poprzedni odcinek [Odcinek 10] zakończyłem na stwierdzeniu, że porównywanie wielkości R2 w dwóch modelach (np. pełnym, ze wszystkimi predyktorami, i ograniczony, ze zmniejszoną liczbą predyktorów) daje nam istotne informacje na temat zmiennych uwzględnionych w modelu, jeśli zmienne te są ze sobą nieskorelowane. W przypadku jednak, gdy istnieje większa lub mniejsza korelacja między poszczególnymi predyktorami (sytuację taką nazywamy współliniowością), zmiana R2 z modelu na model przestaje w sposób jednoznaczny pokazywać nam ważność wprowadzonych do modelu zmiennych. Współliniowość jest niestety dość powszechnym problemem w badaniach psychologicznych i – wstyd powiedzieć – wielu autorów wciąż nie radzi sobie zbyt dobrze z właściwą i pełną interpretacją modeli regresji, w których ten efekt występuje. Wypada więc, byśmy przyjrzeli się tej kwestii nieco dokładniej.
Podobnie jak w poprzednim odcinku, rozważmy całkowicie fikcyjny przykład, w którym zmienną poziom lęku będziemy przewidywać na podstawie wzrostu, szybkości reakcji w jakimś zadaniu motorycznym oraz pojemności pamięci bezpośredniej, przy czym wzrost, szybkość reakcji i pamięć są względem siebie ortogonalne, tj. korelacje między nimi wynoszą zero. Równocześnie, każda z tych zmiennych niech koreluje z lękiem na poziomie r=0.4.
Szacując modele regresji oddzielnie dla każdej z tych zmiennych, otrzymujemy dla wzrostu model z R2=0.16, a standaryzowany współczynnik regresji dla tej zmiennej (tzw. beta) wynosi 0.4. W przypadku szybkości reakcji, R2=0.16, a beta wynosi 0.4, a w przypadku pamięci… R2 wynosi 0.16, a beta 0.4. Jak widać, w każdym przypadku beta jest równa korelacji między predyktorem a zmienną przewidywaną (poziom lęku), a R2, czyli współczynnik determinacji, jest po prostu wartością tej korelacji podniesioną do kwadratu. I tak właśnie być powinno w przypadku jednozmiennowego modelu regresji.
Szacując następnie model regresji uwzględniający wszystkie trzy zmienne, otrzymujemy R2=0.48 (przy R2adj=0.475, a R2pred=0.466), a standaryzowane współczynniki regresji dla tych zmiennych wynoszą odpowiednio 0.4 dla wzrostu, 0.4 dla szybkości reakcji i 0.4 dla pamięci. Na pierwszy rzut oka widać więc, że przy ortogonalnych predyktorach, całkowita zmienność wyjaśniana przez model regresji jest sumą zmienności wyjaśnianej przez poszczególne zmienne, a bety w jedno i wielozmiennowych modelach są identyczne i odzwierciedlają korelacje predyktorów ze zmienną przewidywaną. Łatwo dzięki temu możemy określić, która ze zmiennych w modelu jest najsilniej powiązana ze zmienną, którą chcemy przewidywać (lub – jak w naszym przypadku – stwierdzić, że wszystkie zmienne są równie ważne).
Dokonajmy jednak modyfikacji naszego przykładu – oto wzrost, szybkość reakcji i pamięć wciąż korelują z lękiem na poziomie 0.4 każda, ale dodatkowo niech pamięć koreluje z wzrostem i szybkością reakcji na poziomie 0.6, a wzrost i szybkość reakcji korelują ze sobą minimalnie, na poziomie 0.1. Oszacowanie oddzielnych modeli regresji dla każdej z tych zmiennych dam nam efekty takie same jak poprzednio (każdorazowo R2=0.16, beta wynosząca 0.4), ale już model zawierający wszystkie trzy predyktory równocześnie okazuje się być diametralnie odmienny.
Dla takiego modelu, całkowita wyjaśniona zmienność, tj. R2 wynosi 0.29. Bety dla wzrostu i szybkości reakcji są identyczne i wynoszą 0.42, zaś beta dla pamięci jest ujemna i wynosi -0.11. Współliniowość predyktorów obniżyła więc proporcję wariancji wyjaśnionej i „zamieszała” we współczynnikach regresji. Ten pierwszy fakt jest dość łatwy do wytłumaczenia – skoro predyktory są ze sobą skorelowane, to znaczy, że część ich zmienności to zmienność wspólna. A więc w R2 z modeli jednozmiennowych, w wartości 0.16, tylko pewna część tych szesnastu setnych to unikalna zmienność związana z pojedynczym predyktorem. Jak duża? Najczęściej określamy to, wyliczając delta R2 [Odcinek 10] dla modelu, w którym interesujący nas predyktor dodawany jest jako ostatni. Nazywamy to czasami miarą użyteczności zmiennej, i, co ciekawe, jest ona miarą odpowiadającą współczynnikowi eta2 wyliczanemu dla poszczególnych czynników w analizie wariancji (czyli: ΔR2 dla każdej ze zmiennych włączanych w ostatniej kolejności do modelu jest miarą analogiczną do wartości eta2 dla każdej z tych zmiennych, jaką otrzymamy na podstawie analizy modelu pełnego… ufff).
Dla wzrostu więc będzie to delta R2 z porównania modelu ze wzrostem, szybkością reakcji i pamięcią z modelem zawierającym tylko szybkość reakcji i pamięć, i wynosi ona ΔR2=0.09. Analogicznie wyliczona ΔR2 dla szybkości reakcji także wynosi 0.09, a ΔR2 wyliczona w ten sposób dla pamięci wynosi… 0. Suma wskaźników użyteczności poszczególnych zmiennych, czyli suma delt wynosi 0.18. Pozostała część R2 (czyli 0.29 minus 0.18, a więc 0.11) można uznać za wariancję wspólną dla naszych trzech predyktorów.
Z przykładu tego widać też, że, po uwzględnieniu wzrostu i szybkości reakcji, trzecia zmienna, czyli pamięć, która jest umiarkowanie (r=0.6) skorelowana z poprzedniczkami, nie wnosi już żadnych unikalnych informacji do naszego modelu, czyli nie poprawia przewidywań poziomu lęku, jakie możemy czynić na podstawie wartości dwóch pierwszych zmiennych. W naszym modelu użyteczność tego predyktora więc jest zerowa, i jest to jak najbardziej poprawna interpretacja, choć możemy mieć poczucie, że jest nieco krzywdząca dla naszej zmiennej, która wszak – identycznie jak pozostałe – koreluje z lękiem na poziomie 0.4.
Co więcej – mimo że korelacja z poziomem lęku jest dodatnia (czyli, hipotetycznie, wraz ze wzrostem pojemności pamięci zwiększa się poziom lęku – przypominam, że przykład jest z gatunku absurdalnych), to w pełnym modelu beta dla tej zmiennej ma wartość ujemną, co oznacza, że wraz ze wzrostem pojemności pamięci poziom lęku maleje, przy założeniu oczywiście, że poziom pozostałych zmiennych (wzrost i szybkość reakcji) pozostaje ten sam. Ta interpretacja też jest poprawna…
Jak więc się dzieje, że – zależnie od podejścia – wzrost pojemności pamięci będzie się wiązał z wzrostem albo ze spadkiem poziomu lęku? Odpowiedź na to pytanie kryje się w przedostatnim zdaniu poprzedniego akapitu, a konkretnie we fragmencie: „wraz ze wzrostem pojemności pamięci poziom lęku maleje, przy założeniu oczywiście, że poziom pozostałych zmiennych (…) pozostaje ten sam”. Problem jest zaś taki, że jest to sytuacja bardzo hipotetyczna – skoro pamięć jest skorelowana z pozostałymi predyktorami, to zmiana wartości tej zmiennej może pociągać za sobą zmianę wartości pozostałych predyktorów, a więc cytowane założenie nie daje się utrzymać!
Wnioski z powyższego przykładu są zarówno pocieszające, jak i zasmucające. Te weselsze są związane z tym, że choć współliniowość zmiennych wpływa na oszacowanie parametrów regresji, czyniąc bety trudnymi lub wręcz niemożliwymi do sensownej interpretacji, to ogólne dopasowanie modelu wyrażone poprzez całkowite R2 (czy jego skorygowane wersje) nie jest w znaczący sposób tą współliniowością obciążone i rzetelnie pokazuje nam, jak dużą część zmienności wyników jesteśmy w stanie powiązać ze zbiorem wybranych przez nas predyktorów.
Mniej wesoły wniosek dotyczy zaś faktu, że w modelu takim oszacowanie znaczenia poszczególnych zmiennych robi się bardzo trudne – można oczywiście postępować pragmatycznie i patrzeć po prostu na ich użyteczność, ale w ten sposób możemy wyeliminować zmienne, które są merytorycznie (np. w związku z teorią psychologiczną) zmiennymi ważnymi, tyle, że – pechowo – inne, „silniejsze” zmienne eliminują je z ostatecznego modelu. Może być też tak, że jakaś zmienna okaże się statystycznie mało „użyteczna” w przewidywaniu interesującej nas cechy, ale w praktyce to właśnie wartość tej zmiennej podlega największym zmianom u osób badanych, i dopiero wtórnie „pociąga za sobą” zmianę wartości pozostałych zmiennych (np. mając model ze skorelowanymi predyktorami takimi jak wiek i wzrost, może okazać się, że wzrost, jako lepszy predyktor, wyparł nam z modelu wiek – choć dobrze wiemy, że zmiany we wzroście wynikają ze zmian wieku, a nie odwrotnie, i na poziomie interpretacji wyników trudno byłoby tą zmienną zupełnie zignorować…).
Pewnym rozwiązaniem, które nieco bardziej precyzyjnie niż użyteczność pozwala ocenić, na ile znaczące są poszczególne predyktory w przewidywaniu interesującej nas zmiennej, są tzw. miary relatywnej ważności dla współczynników regresji. Wadą użyteczności (czyli ΔR2 dla zmiennych wprowadzanych do modelu na ostatniej pozycji) jest m.in. to, że wartości tej miary dla poszczególnych zmiennych nie sumują się i nie mają bezpośredniego przełożenia na R2 całego modelu (w naszym modelu sumowanie wartości ΔR2 daje nam w efekcie wartość 0.18, czyli znacznie mniej, niż 0.29 dla modelu pełnego).
Jednym z celów stosowania miar relatywnej ważności jest więc poradzenie sobie z tym problemem i pokazanie „proporcjonalnego wkładu” każdego z predyktorów do całkowitego R2. Miary te, choć nie są one pomysłem szczególnie nowym, są jednak wciąż mało popularne, a ich największy rozwój przypada głównie na ostatnie dwudziestolecie. Jest to prawdopodobnie związane z tym, że oparte są one na dość złożonych i wymagających sporych zasobów algorytmach obliczeniowych.
Najlepiej ugruntowaną w literaturze, a przez to najczęściej polecaną miarą tego typu jest współczynnik lmg, nazwany tak od nazwisk autorów (Lindeman, Merenda i Gold, 1980). W dużym uogólnieniu, współczynnik ten, dla każdej ze zmiennych, jest uśrednioną miarą jej kolejnych wartości ΔR2 uzyskiwanych w modelach regresji z różnymi konfiguracjami predyktorów. W przypadku naszego przykładu, lmg dla wzrostu wynosi 0.11, dla szybkości reakcji także 0.11, a dla pamięci 0.07. Jak widać, potwierdza się wniosek o tym, że w naszym, składającym się z trzech predyktorów modelu pamięć odgrywa rolę najmniejszą, jej wkład (związany z częścią zmienności, którą dzieli z innymi predyktorami) został jednak „doceniony”. Zsumowanie wartości lmg daje nam zaś wartość 0.29, czy R2 naszego modelu pełnego .
Nowszą alternatywę stanowi miara PMVD opisana przez Feldmana (2005), która – zachowując ogólne cechy miar relatywnej ważności – miała na celu wyeliminowanie tej właściwości miary lmg, która sprawiała, że nawet predyktory o beta równym 0 (ale znacząco skorelowane z pozostałymi predyktorami) zostawały „obdzielone” sporym kawałkiem wyjaśnianej przez model zmienności. W naszym przypadku PMVD dla wzrostu wynosi 0.14, dla szybkości reakcji podobnie, zaś dla pamięci 0.01. Jak widać, wnioski są tu podobne jak na podstawie analizy z zastosowaniem ΔR2, ale – po pierwsze – rola pamięci, choć minimalna, nie jest wyrugowana całkowicie; pod drugie zaś, w odróżnieniu do ΔR2 i podobnie jak przypadku lmg, suma wartości PMVD wynosi 0.29, czyli daje nam wartość całkowitego R2 naszego modelu regresji. Co ciekawe, algorytm służący do wyliczania PMVD został w USA opatentowany i z tego względu jego rozpowszechnianie w postaci darmowych skryptów do pakietów statystycznych czy darmowego oprogramowania jest na terenie Stanów Zjednoczonych niezgodne z prawem, co – na chwilę obecną – znacząco ogranicza jego popularność.
Jedną z najnowszych propozycji jest zaś parametr CAR (Zuber i Strimmer, 2011), który autorzy sugerują interpretować jako „korelacje brzegowe skorygowane o korelacje między predyktorami” (zainteresowanych, co to dokładnie znaczy, odsyłam do oryginalnego artykułu). W naszym przypadku, CAR dla wzrostu wynosi 0.13, dla szybkości reakcji podobnie, zaś dla pamięci 0.03. Wartości te ponownie sumują nam się od 0.29 i lokują się gdzieś między wartościami lmg i PMVD, choć są zdecydowanie bliższe tych drugich.
Podsumowując – trudno definitywnie wskazać, jaka metoda oceny „ważności” predyktorów w równaniu regresji jest najlepsza w sytuacji, gdy predyktory te są ze sobą skorelowane. Dużo zależy od pytań, jakie zadajemy (czy chodzi nam np. o wkład w wariancję wyjaśnioną przez model, czy o bezpośrednią siłę związku ze zmienną przewidywaną czy może o to, jak dużą zmianę w zmiennej przewidywanej pociąga za sobą zmiana wartości naszego predyktora). Problem jest złożony i najlepiej podejść do niego z różnych perspektyw. Po tym krótkim tekście powinno być jednak jasne, że najbardziej popularne metody, czyli porównywanie wartości beta (metoda całkowicie chybiona!) czy porównywanie użyteczności zmiennej w modelu (ΔR2) to jedynie wierzchołek góry lodowej i drobna część możliwości, jakie w tym zakresie oferuje nam współczesna statystyka.
dr Piotr Zieliński
Wojskowy Instytut Medycyny Lotniczej
Wszystkie opisane w tekście miary relatywnej ważności predyktorów zostały obliczone za pomocą pakietu relaimpo (Groemping, 2006) w środowisku statystycznym R.
Literatura cytowana:
Feldman, B. (2005). Relative importance and value. Manuscript version 1.1.
Groemping, U. (2006). Relative importance for linear regression in R: the package relaimpo. Journal of Statistical Software, 17, 1-27.
Lindeman, R.H, Merenda, P.F., Gold, R.Z. (1980). Introduction to bivariate and multivariate analysis. Scott, Foresman, Glenview, IL.
Zuber, V., Strimmer, K. (2011). High-dimensional regression and variable selection using CAR scores. Statistical Applications in Genetics and Molecular Biology, 10(1).