czwartek, 2 kwietnia 2015

Autokorelacja przestrzenna w modelowaniu


Problem statystyczny, może ktoś pomoże?
Mam dane z liczeń punktowych z 5 kompleksów łąk. W każdym z nich jest ok. 20-30 punktów, w których liczono ptaki. Najbliższe punkty oddalone są od siebie o ok. 300m, a całe kompleksy łąk oddalone o setki kilometrów. 
Chcę uwzględnić autokorelację przestrzenną, ale wydaje mi się, że zwykłe wrzucenie współrzędnych każdego z punktów do jakiejś poprawki w modelach GLM nie będzie chyba dobrym rozwiązaniem, bo efekt lokalizacji geograficznej jest napędzany głównie oddaleniem poszczególnych kompleksów łąk. A może się mylę?? Standaryzacja w obrębie kompleksu również odpada (zależy mi na skonstruowaniu jednego modelu dla wszystkich 5 kompleksów łąk).

będę wdzięczny za podpowiedzi,
michał żmihorski

19 komentarzy:

  1. Cześć,
    Po pierwsze sprawdziłbym, czy istnieje autokorelacja osobno dla każdego z pięciu kompleksów łąk. Jeśli nie - użyłbym zwykłego modelu mieszanego z id powierzchni jako czynnikiem losowym (który, obejmuje w pewien sposób autokorelacje przestrzenna w tej większej skali, między kompleksami). Jeśli jednak autokorelacja jest obecna w obrębie każdego z tych 5 kompleksów łąk, to wtedy użyłbym jakiegoś modelu przestrzennego, który obliczyłby autokowariatę opisującą zmiany zagęszczenia(?) w obrębie każdego kompleksu. I taką sumaryczną kowariatę włączyłbym do modelu mieszanego (być może wraz z interakcją między tą kowariatą a id kompleksu). Ale może są jakieś lepsze rozwiązania tego problemu.
    Pozdrawiam
    Piotrek

    OdpowiedzUsuń
    Odpowiedzi
    1. Aha, tą autokowariatę policzyłbym z jakiejś spatial autoregression, np. z lagged model.

      Usuń
  2. sprawdziłem na początku tak jak powiedziałeś, Moran's I istotny dla każdego kompleksu

    OdpowiedzUsuń
  3. zasadniczo problem jest prosty. w modelu liniowym zamiast wariancji resztkowej (residual) takiej ze poszczególne residuale są i.i.d. (czyli de facto mają rozkład o wariancji I{i,j} x sigma_e) wystarczy zdefiniować strukturę wariancji gdzie zamiast I{i,j} jest macierz kowariancji z elementami off-diagonal zależnymi od odległości (takich zależności jest kilka rodzajów, konkretne modele można porównać za pomocą AIC). łatwo w R można to zrobić używając opcji correlation w nlme - ALE to pozwoli na modelowanie autokorelacji dla danej powierzchni. Żeby można to było zrobić dla wszystkich na raz macierz (ko)wariancji resztkowych (residual) musi mieć strukturę blokową czyli na przekątnej mamy w/w macierze (z autokorelacją przestrzenną) dla każdej powierzchni o poza przekątną zera (bo (?) zakładamy że autokorelacja punktów znajdujących się w 2 różnych z 5 kompleksów jest zero - nie wiem czy o to Ci chodzi?; oczywiście jeśli w skali tych setek km też ma być autokorelacja można ją także zdefiniować). i tutaj obawiam się, że nlme nie da rady, takie modele są dość trudne, powinien dać radę ASReml ale nie jest free. uwzględnienie jakiejkolwiek "poprawki" w modelu (po stronie efektów ustalonych) nie bardzo rozwiąże cokolwiek bo efekty ustalone nie "leczą" korelacji residuali więc nawet po uwzględnieniu czegoś takiego residuale mogą ciągle być z dużym prawdopodobieństwem skorelowane przestrzennie, czego oczywiście nie chcemy.

    OdpowiedzUsuń
    Odpowiedzi
    1. jest jak napisałeś: autokorelacja między pięcioma kompleksami łąk jest równa zero (tak zakładam), interesuje mnie jedynie poprawka na autokorelację w obrębie kompleksów, ale w jednym modelu, uwzględniającym tych 5 kompleksów.
      Twój komentarz, że "nlme nie da rady" odnosi się do tej struktury blokowej właśnie?
      dzięki!

      Usuń
    2. tak - problem polega na tym, że większość procedur w R nie radzi sobie z takimi strukturami (możesz googlować "block-diagonal matrix"). potrafi to na 100% MCMCglmm - ale z kolei jeszcze nie ma w nim opcji definiowania autokorelacji residuali. Znam funkcję pdBlocked która tworzy taką (positive-definite) macierz w formie w jakiej można to wsadzić do nlme ale nie używałem tego więc nie pomogę ;) nie bardzo wiem czy pdBlocked potrafi współpracować z corXXX (czyli funkcjami które obsługują różne postaci autokorelacji; ?corClasses) - z reguły używam ASReml do podobnych potworków.

      Usuń
  4. Szkoda, że nie napisałeś, co jest zmienną zależną :) Najprościej wyrzucić część punktów, tak żeby pozbyć się autokorelacji. Ale po tym zrobisz się smutny. Możesz więc (w każdej lokalizacji osobno) odsiewać stopniowo punkty (np. co 300 m, co 600m itd), aż do momentu, gdy autokorelacja zmiennej zależnej będzie słaba. A w drugim kroku zrobić bootstrap modelu, w każdej iteracji używając wylosowanego zestawu punktów spełniających kryteria zdefiniowane w pierwszym kroku. Użyjesz wtedy wszystkich danych (bo np. kryterium '> 600 m' da za każdym razem inny zestaw, w zależności od tego, gdzie zaczynasz wybierać) pomijając wpływ autokorelacji w przestrzeni. KH

    OdpowiedzUsuń
    Odpowiedzi
    1. ok, przepraszam - zmienną zależną jest liczebność ptaków w danym punkcie.
      A podrzuciłbyś jakieś namiary na literaturę albo przykładową pracę wykorzystującą tę metodę?
      Dzięki za sugestię

      Usuń
  5. Myślałem jeszcze o wrzuceniu liczebności badanego gatunku w punktach sąsiednich jako dodatkowej kowariaty do modelu. Wtedy uciekam od problemu przestrzennego skupienia punktów w pięciu kompleksach łąkowych, ale to zabieg dość prymitywny...

    OdpowiedzUsuń
    Odpowiedzi
    1. I niestety prawie na pewno nie rozwiąże problemu autokorelacji całkowicie ;)

      Usuń
  6. Taka kowariata na pewno nie pomoże, chyba że autokorelacja jest istotna tylko dla najbliższych sąsiadów. Mi przychodzi jeszcze inny, też dosyć prymitywny pomysł. Można by sprawdzić, dla każdego kompleksu łąk, na wykresach wskaźnika autokorelacji Morana, przy jakich dystansach jego wartość spada do zera. Następnie można podzielić kompleks łąk na bloki, o wielkości odpowiadającej temu dystansowi (odpowiednio zmienionego dla każdego kompleksu). Uzyskujesz wtedy układ hierarchiczny, z blokami zagnieżdżonymi w kompleksie. Oczywiście, coś takiego można zrobić, kiedy autokorelacja przestrzenna ma "klasyczną" postać (podobieństwo spada wraz z odległością do zera). Kiedyś, czytałem prackę, gdzie autorzy stosowali coś podobnego w GLMM, ale nie potrafię sobie przypomnieć jej tytułu.
    Może znajdziesz coś ciekawego w tym artykule: http://www.whoi.edu/cms/files/DormannEcography30_57164.pdf

    OdpowiedzUsuń
  7. > A podrzuciłbyś jakieś namiary na literaturę
    > albo przykładową pracę wykorzystującą tę metodę

    Nie przypominam sobie pracy, w której w taki właśnie sposób ktoś obchodził problem niechcianej autokorelacji w przestrzeni. Jestem jednak przekonany, że można by było tak zrobić i formalnie będzie to ok. Ale no guarantees. Jak ktoś potrafi myśleć na wyższym poziomie abstrakcji, to pewnie będzie mu się taka numeryczna sieczka wydawała toporna.

    W podręczniku Foxa "Applied regression analysis and generalized linear models'" jest rozdział "Bootstrapping regression models", na 20 stron. KH

    OdpowiedzUsuń
  8. Być może "autokorelacja" zmiennej zależnej jest indukowana przez jakieś zautokorelowane czynniki zewnętrzne (którąś ze zmiennych niezależnych?). Należałoby zbudować model bez uwzględniania autokorelacji, a następnie sprawdzić, czy resztki z modelu są zautokorelowane w obrębie kompleksów łąk. Jeśli nie, problem nie istnieje. Robiłeś tak?

    OdpowiedzUsuń
    Odpowiedzi
    1. dzięki - tak, sprawdzam autokorelację Morana dla residuali z pełnego modelu.
      Przy okazji pytanie: zakładam, ze nie ma sensu sprawdzać tej korelacji dla kompleksów łąk, na których w ogóle nie było gatunku??

      Usuń
  9. to jeszcze mam pytanie: czy wartość statystyki Morana można jakoś rozsądnie interpretować? Twierdzi się na przykład, że korelacje powyżej 0.7 utrudniają modelowanie wielowymiarowe, albo VIF>10 jest już nieakceptowalny itp. Czy z Moranem jest podobnie?? Mam dla kilku gatunków I Morana istotnego ale bezwzględne wartości są bardzo małe, rzędu 0.1
    Nie mogę nic znaleźć w necie na szybko...

    OdpowiedzUsuń
    Odpowiedzi
    1. Można interpretować tak samo jak współczynnik korelacji. Teraz widzę, że to co obliczyłeś jest globalnym współczynnikiem Morana, który opisuje przeciętną wartość autokorelacji w terenie. Jeżeli jest bliski wartości 1, to oznacza, że wszystkie wartości w sąsiednich punktach są podobne (takie same). Jeżeli wartość wynosi - 1, to oznacza, że sąsiednie wartości są kompletnie niepodobne. Jeżeli wartość statystyki Morana wynosi 0, to sąsiednie wartości badanej zmiennej w sąsiednich punktach przyjmują wartości losowe (równie dobrze mogą być podobne, jak i inne). Analogicznie do zwykłej korelacji - jeśli liczebność próby jest dostatecznie duża, to i niewielka wartość autokorelacji może być istotna. Tak na marginesie, wydaje mi się, że globalny współczynnik Morana (o ile taki był liczony), nie jest najlepszym rozwiązaniem (tzn. mówi niewiele). Być może lepiej byłoby konstruować semi-wariogramy albo liczyć współczynnik Morana dla różnych dystansów. Wynika to z tego, że wartość autokorelacji może być wysoka przy niewielkich dystansach od badanego punktu, i maleć wraz z odległością. Jest to tzw. lokalna autokorelacja. Niektórzy uważają, że dlatego lepiej jest używać współczynnika "C" Geary-ego.

      Usuń
    2. dzięki Piotrek! tak chyba właśnie robiłem, bo liczyłem Morana dla 5 klas dystansu (liczyłem też dla innej liczby klas, ale to niewiele zmienia). Generalnie jest tak jak piszesz, że dla pierwszej klasy jest najwyższy, a dla kolejnych już płaski w okolicy zera.
      To jeszcze pytanie: czy wartość współczynnika Morana daje mi coś zbliżonego do współczynnika determinacji (tj. czy mogę skwadratowanym Moranem określić jaką część zmienności wektora jest wyjaśniana przez jego wartość w sąsiadującej klasie odległości)?
      dzięki

      Usuń