sobota, 12 września 2015

Parametryzacja w modelach liniowych - alternatywa dla B i Bety?

Update:
przykład pokazujący, że przeliczenie parametru B dla całego range nie zmniejsza wiarygodności tego parametru, a jedynie ułatwia jego interpretację (pomińmy bezsensowny intercept...):
Wykres 1: las w ha i standardowy parametr B (wzrostu wielkości lasu o 1ha daje o 0.0933 jarząbka więcej)

Wykres 2: to samo, ale parametr B przemnożony przez range (wzrostu wielkości lasu od najmniejszego do największego na danym terenie daje o 9.33 jarząbków więcej)
Wykres 3: las w kilometrach kwadratowych, standardowy parametr B (dla wzrostu wielkości lasu o 1km2 mamy o 9.33 jarząbków więcej)

Na czym polega wada parametru z modelu 2? Nie widzę tu wad innych niż dla "zwykłych" parametrów, a ma trochę zalet - patrz lista poniżej, ale przypomnę: B*range jest zawsze ten sam dla różnych jednostek predyktora (ha,km2,akr,mila morska) i jest wyrażony w jednostkach zmiennej zależnej, więc ułatwia interpretację (tu: w jarząbkach). Beta (policzona funkcją "lm.beta" z pakietu "QuantPsyc") wynosi: 0.9505 (niezależnie czy używamy ha czy km2) ale przyznam, że nie jest łatwa w interpretacji (poza tym wymaga normalności predyktora). A dziewięć jarząbków, to jest jakiś konkret :-)

I jeszcze cytat ze StackExchange, akurat znalazłem:
...the idea that standardizing independent variables makes it easier to compare the effects of one variable to another. This advantage is, in my opinion, somewhat illusory, since it depends on the range of data in your sample. Although it's a matter of some contention, I am generally against standardizing variables. Variables themselves are, in my view, easier to interpret than standard deviations of variables - we often have an intuitive sense about variables themselves.
 

Problem roboczy, niezbyt dobrze przemyślany:
parametry modeli liniowych są najczęściej prezentowane jako zmiana wartości zmiennej zależnej w reakcji na wzrost wartości danego predyktora o jedną jednostkę (standardowy output z Ra, SPSSa itp). Problem z interpretacją takiej tabeli polega na tym, że różne predyktory są wyrażane w różnych jednostkach (np. powierzchnia lasu i temperatura), co nie jest zbyt odkrywcze, oraz na tym, że różne predyktory mają różną "rozpiętość". 

Przykład: sprawdzamy wpływ powierzchni lasu i liczby strumieni na obecność jarząbka. Powierzchnię lasu mamy w hektarach, a liczba strumieni opisuje po prostu ile ich jest w danym lesie (wbrew pozorom te zmienne nie są mocno skorelowane). Z tym, że powierzchnia waha się od 5 do 150, a liczba strumieni od 0 do 5 (wartości zmyślone). Z modelu otrzymamy info jak rośnie szansa stwierdzenia gatunku wraz ze wzrostem pow lasu o 1 hektar i wzrostem liczby strumieni o 1 strumień. I tu moja wątpliwość: czy nie jest bardziej informatywne przedstawienie wzrostu szansy stwierdzenia gatunku wraz ze wzrostem powierzchni lasu o 145ha, zamiast o 1ha (a liczby strumieni o 5, zamiast 1)? Czyli dla pełnego zakresu zmienności danego predyktora, a nie o jedną jego jednostkę. A nawet jeśli takie przedstawienie parametrów nie jest substytutem tradycyjnego, to może warto je prezentować równolegle do niego?

Jakie są zalety takiej metody:
* możemy bezpośrednio porównać znaczenie predyktorów wyrażonych w różnych jednostkach i w tych samych jednostkach (być może nawet mówić o istotności różnic, na podstawie 95%CI - muszę o tym pomyśleć),
* mamy lepszą informację jak ważny jest faktycznie dany predyktor w rzeczywistości, bo pokazujemy jego potencjalny sumaryczny wpływ dla całej rozpiętości wielkości lasów, która faktycznie ma miejsce w danym krajobrazie.

Zrobiłem takie przeliczenie kiedyś w pracy o jarząbku właśnie, w publikacji z Łukaszem Kajtochem, bo wydawało mi się to rozsądne. Wygląda to tak (ostatnia kolumna, czyli B*range):
Nadal wydaje mi się to rozsądne ale bardzo możliwe, że popełniam błąd w rozumowaniu? Beta nadaje się do porównywania predyktorów ale same jej wartości są chyba mniej informatywne niż wartości B (ile nam mówi SD średniej powierzchni lasu?). A może istnieje jakaś oczywista alternatywa, której nie znam, lub nie kojarzę z tym problemem?
Będę wdzięczny za komentarz,

michał żmihorski


12 komentarzy:

  1. Używanie 'range' jest ryzykowne, czyni ten iloczyn podatny na wartości skrajne. Lepsze już byłoby używanie IQR.
    Ale aby oceniać siłę wpływu to raczej używa się % wyjaśnionej wariancji (w tabeli ANOVA) i jednostka dla X nie ma znaczenia.

    OdpowiedzUsuń
  2. Też wydaje mi się, że używanie zasięgu "range" nie jest najlepsze, właśnie ze względu na wartości skrajne, które niekoniecznie odzwierciedlają to co dzieje się wewnątrz zbioru obserwacji. Dlaczego nie postąpić "standardowo" i zestandaryzować wszystkich zmiennych ciągłych, i wyrazić ich jednostkach odchylenia standardowego? To generalnie bardzo dobrze pokazuje jak silny jest wpływ danej zmiennej niezależnej. Inna metoda, którą często stosuję jak mam zmienne ciągłe, to "hierarchical partitioning" - co pozwala określić ile % z wyjaśnionej wariancji można przypisać indywidualnemu wpływowi danej zmiennej, a ile wynika z łącznego wpływu wszystkich zmiennych (pakiet "hier.part" w R).

    OdpowiedzUsuń
  3. Przemek Chylarecki12 września 2015 23:43

    To nie jest "SD z powierzchni lasu" jak piszesz, ale SD z powierzchni lasu #w Twojej próbie#, czyli bierzesz pod uwagę "range" (jak chciałeś), czyli uwzględniasz specyfikę Twojej próby. Dlatego nie mam problemu/oporów z używniem beta. Mają one dokładnie te zalety, które wypunktowałeś jako poszukiwane

    OdpowiedzUsuń
    Odpowiedzi
    1. No tak, racja, SD jest z mojej próby. Ale czy widzisz jakieś konkretne wady B*range? I czy nie dostrzegasz, że niesie ona klarowną informację, bezpośrednio odpowiadając na pytania często stanowiące cel danej pracy?

      Usuń
  4. Myślę, że wartości skrajne nie mają tu znaczenia - jeśli uwzględniliśmy je w modelu i uznaliśmy, że nie są odstające, to konsekwentnie powinniśmy ich używać. No i używamy, bo parametr B również odnosi się - by użyć przykładu powyższego - do wzrostu jarząbka, gdy powierzchnia lasu rośnie z 149 do 150ha. Zakładamy zatem liniowy wzrost szans dla całej rozpiętości tej cechy, a mój pomysł w 100% bazuje na tym założeniu i jedynie inaczej prezentuje parametr B. Zobaczcie, że gdyby lasy wahały się od 1-330ha to zabieg który proponuję jest niczym innym tylko zmianą jednostki w której mierzymy lasy: jeśli zamiast hektarów używamy mili kwadratowej, to powierzchnia lasu wahałaby się od 0.003-1. Wtedy parametr B (przyrost jarząbka przy skoku lasu o 1 unit) byłby tym samym czym obecnie jest B*range.
    Ok, rozumiem że można patrzeć na wariancje wyjaśnianą, to jest dobre podejście. Standaryzacja też jest ok, ale czy takie parametry nie są mniej informatywne? Nie chodzi mi nawet o porównywanie znaczenia predyktorów, mnie głównie zachęca właśnie klarowność informacji, którą niesie parametr B*range: "między najmniejszym a największym lasem w danym krajobrazie szansa na stwierdzenie jarząbka rośnie 48-krotnie"
    Często właśnie o to pytamy rozpoczynając badania.

    OdpowiedzUsuń
    Odpowiedzi
    1. OK, rzeczywiście ma to sens. To co jest plusem tego B*range, może jednak być odbierane jako minus, bo w wielu przypadkach interesujący jest sam relatywny wpływ zmiennych niezależnych, na zależną. Ale jeśli podawać zarówno B oraz B*range ilość informacji jest faktycznie bardzo duża.

      Usuń
    2. Innymi słowy, słabość jaką widzę to fakt- że samo range może bardzo różnić się między predyktorami (nawet wyrażając je zestandaryzowanymi jednostkami). Weźmy hipotetyczny przykład np, range w przypadkow wielkości lasow może być 100, a w przypadku liczby dolin np. 50. Beta mogą wynosić 2 dla liczby lasów i 3 dla liczby dolin. B*range daje odpowiednio 200 i 150. Może to sugerować, że wielkość lasu jest ważniejsza od liczby dolin, co wcale nie musi być prawdą. Dlatego porównywanie ważności różnych predyktorów tym sposobem może być błędne, choć sama B*range niesie sporo informacji o próbie. Dobrze rozumuję?

      Usuń
    3. Nie wiem... Jeśli B*range dla lasów jest większe niż dla dolin, to de facto, na terenie badań lasy mają większe znaczenie niż doliny. A jakie znaczenie ma wysokość nad poziomem morza - może mieć ogromne, ale jeśli nasz teren badań jest płaski, to będzie bez znaczenia. Nasze wyniki są zawsze bardzo obciążone zmiennością predyktorów (jakoś nie spotykam zbyt wiele dyskusji na ten temat) - możemy przecież znaleźć lasy o tej samej wielkości i wtedy wielkość lasu będzie miała wpływ zerowy, a całą wariancję wyjaśnią doliny... Czy jesteśmy w stanie wystandaryzować te parametry tak, żeby powiedzieć, że generalnie (=niezależnie od danego terenu badań) np. lasy są ważniejsze niż liczba dolin (i czy w ogóle jest taka potrzeba?). Beta również jest zakotwiczona w terenie badań, a poza tym używa nieinterpretowalnych jednostek i zakłada normalność rozkładu predyktora.
      Jutro spróbuję o tym pogadać z dobrym statystykiem tu na SLU, to wrzucę coś może.

      Usuń
    4. No właśnie, stajemy przed dylematem, co jest ważniejsze - zmienność predyktora, czy sam predyktor ;-)

      Usuń
  5. Problem z range jest też taki, że dotyczy próby a nie populacji. Chcesz wnioskować o populacji a akurat range to bardzo kiepski (=obciążony) estymator.

    W przykładach: badasz zależność od powierzchni analizując 10 lasów. Czy za range weźmiesz min i max powierzchni lasów z tej dziesiątki, czy też min i max powierzchni lasów w Polsce (takie statystyki są dostępne) czy też min i max powierzchni lasów na świecie, czy też min i max lasów w których żyją jarząbki?

    A w ogóle czy minimalna wielkość lasu to nie jest 0? Czy też jest jakaś ustawowa wielkość minimalna o której wiedzą jarząbki.

    Min i max liczony na podstawie tych 10 lasów są bliżej siebie niż pozostałe możliwe min/max'y.

    OdpowiedzUsuń
    Odpowiedzi
    1. oczywiście biorę min i max lasów w próbie, nie chcę wychodzić z ekstrapolacją poza badany zakres, bo jest to ryzykowne (np. zależność może przestać być liniowa).
      Ależ range jest tym samym estymatorem co B w równaniu! Czy zmieniając jednostkę z metrów na kilometry zmienia Ci się obciążenie estymatora w równaniu? Moim zdaniem niesłusznie zakładasz, że range ma jakieś dodatkowe obciążenia, których nie ma B, ale może czegoś nie rozumiem.

      Usuń
  6. IMHO odpowiedź na pytanie jak standaryzować, zależy od tego co chce się uzyskać. Ja tutaj widzę przynajmniej dwie możliwości.

    1. Użycie nowego znormalizowanego współczynnika aby lepiej ocenić, które czynniki w modelu mają większy wpływ na Y
    2. Użycie nowego znormalizowanego współczynnika aby pokazać za jaki maksymalny rozrzut Y on odpowiada.

    W symulacjach losowałeś powierzchnię z rozkładu jednostajnego więc oba powyższe podejścia wydają się podobne, ale jeżeli będziesz losował np. z rozkładu wykładniczego to okaże się że się różnią.

    Przykładowo, załóżmy że chodzi o wpływ wieku na koszt leczenia (na potrzebę liczenia składki) i załóżmy że jest on liniowy - im starsze oosby tym wyższy koszt, 10 pln na rok życia.

    Jakie są możliwości?

    A. zakładamy że 99% naszych klientów ma wiek od 20 do 40 lat, więc wpływ wieku na składkę jest w granicach intercept + 20*10 - intercept + 40*10
    B. zakładamy, że najmłodszy mozliwy klient może mieć 18 lat ale najstarszy możliwy może mieć 100 więc maksymalna możliwa rozpiętość to intercept + 18*10 - intercept + 100*10
    C. budujemy model na grupie osób, które akurat miały wiek od 19 do 75 lat, więc 'range' z próby odpowiada rozpiętości intercept + 19*10 - intercept + 75*10

    Które z tych podejść jest lepsze?

    Jeżeli chcemy normować po to by porównywać parametry w modelu, wtedy ważny jest rozkład wieku naszych klientów i % wyjaśnionej wariancji będzie pewnie bardzij podążął za pierwszą z ww opcji (A)

    Jeżeli chcemy znać maksymalny rozrzut za jaki odpowiada parametr, to mamy opcję (B), wiek może się maksymalnie zmieniać o 82 lata. Choć narazamy się na zarzut, że zależność nie jest liniowa tam gdzie jej nie widzimy, czyli skoro najstarsza osoba w grupie ma 75 lat to czemu za mak wieku klienta brać 100 lat.


    Ja akurat jestem zwolennikiem metod graficznych, więc rozwiązałbym ten problem wykresem.
    Rysując rozkład wartości hat(beta)*x_i. Czyli pokazując jak wygląda udział czynnika beta*x dla takiego rozkładu x'ów jaki jest w próbie.
    Z takiego wykresu można by odczytać i min*beta - max*beta, ale równiez byłoby widać jaki jest typowy wpływ czynnika beta_x na model.

    OdpowiedzUsuń