Jak mogę przeprowadzić ważoną geograficznie analizę głównych składników za pomocą ArcGIS, Python i SPSS / R?

Jestem za opisem / metodologią przeprowadzania analizy ważonych geograficznie głównych składników (GWPCA). Z przyjemnością używam Pythona do dowolnej części tego i wyobrażam sobie, że SPSS lub R są używane do uruchamiania PCA na zmiennych ważonych geograficznie.

Mój zestaw danych składa się z około 30 niezależnych zmiennych, które są mierzone w ~ 550 obszarach spisu powszechnego (geometria wektorowa).

Wiem, że to załadowane pytanie. Ale kiedy szukam i szukam, wydaje się, że nie ma żadnych rozwiązań. Spotkałem równania matematyczne, które wyjaśniają podstawowy skład GWPCA (i GWR). To, czego szukam, jest bardziej stosowane w pewnym sensie, że szukam, jakie główne kroki muszę wykonać, aby przejść od surowych danych do wyników GWPCA.

Chciałbym rozwinąć pierwszą część tej edycji ze względu na uwagi otrzymane poniżej.

Zwracając się do Paula ...

Zainteresowanie GWPCA opieram na następującym artykule:

Lloyd, CD, (2010). Analiza cech populacji za pomocą analizy ważonych geograficznie głównych składników: studium przypadku Irlandii Północnej w 2001 r. Komputery, środowisko i systemy miejskie, 34 (5), s. 383–399.

Dla tych, którzy nie mają dostępu do literatury, załączyłem zrzuty ekranu poszczególnych sekcji, które wyjaśniają matematykę poniżej:

I żeby zająć się whuber ...

Nie wchodząc w szczegóły (poufność), staramy się zredukować 30 zmiennych, które naszym zdaniem są bardzo dobrymi wskaźnikami (choć globalnie), do zestawu składników o wartościach własnych większych niż 1. Obliczając składniki ważone geograficznie, staramy się zrozumieć lokalne wariancje wyjaśnione przez te elementy.

Myślę, że naszym głównym celem będzie udowodnienie koncepcji GWPCA, to znaczy ukazanie przestrzennie wyraźnego charakteru naszych danych i że nie możemy uznać wszystkich zmiennych niezależnych za wyjaśniające w skali globalnej. Raczej skala lokalna (dzielnice), które zidentyfikuje każdy komponent, pomoże nam zrozumieć wielowymiarowy charakter naszych danych (w jaki sposób zmienne można łączyć ze sobą w celu wyjaśnienia określonych dzielnic na naszym obszarze badań).

Mamy nadzieję zmapować procent wariancji uwzględniony przez każdy komponent (osobno), aby zrozumieć zasięg sąsiedztwa wyjaśniony przez dany komponent (pomóc nam zrozumieć lokalną przestrzenność naszych komponentów). Być może kilka innych przykładów mapowania, ale w tej chwili nie przychodzą mi na myśl żadne.

Dodatkowo:

Matematyka GWPCA wykracza poza to, co rozumiem, biorąc pod uwagę moje doświadczenie w analizie geograficznej i statystyce społecznej. Najważniejsze jest zastosowanie matematyki, czyli co mogę podłączyć do tych zmiennych / formuł.

„Geograficznie ważona PCA” jest bardzo opisowa: w Rprogramie praktycznie się pisze. (Potrzebuje więcej wierszy komentarza niż rzeczywistych wierszy kodu).

Zacznijmy od odważników, ponieważ w tym miejscu ważona geograficznie firma produkuje części PCA od samego PCA. Termin „geograficzny” oznacza, że ​​wagi zależą od odległości między punktem bazowym a lokalizacjami danych. Standardowa - ale w żadnym wypadku nie tylko - waga jest funkcją Gaussa; to jest rozkład wykładniczy z kwadratową odległością. Użytkownik musi określić szybkość zaniku lub - bardziej intuicyjnie - charakterystyczną odległość, na której występuje ustalona wielkość zaniku.

distance.weight <- function(x, xy, tau) {
  # x is a vector location
  # xy is an array of locations, one per row
  # tau is the bandwidth
  # Returns a vector of weights
  apply(xy, 1, function(z) exp(-(z-x) %*% (z-x) / (2 * tau^2)))
}

PCA ma zastosowanie do macierzy kowariancji lub korelacji (która jest pochodną kowariancji). Tutaj jest zatem funkcja do obliczania ważonych kowariancji w sposób stabilny numerycznie.

covariance <- function(y, weights) {
  # y is an m by n matrix
  # weights is length m
  # Returns the weighted covariance matrix of y (by columns).
  if (missing(weights)) return (cov(y))
  w <- zapsmall(weights / sum(weights)) # Standardize the weights
  y.bar <- apply(y * w, 2, sum)         # Compute column means
  z <- t(y) - y.bar                     # Remove the means
  z %*% (w * t(z))  
}

Korelację wyprowadza się w zwykły sposób, stosując standardowe odchylenia dla jednostek miary każdej zmiennej:

correlation <- function(y, weights) {
  z <- covariance(y, weights)
  sigma <- sqrt(diag(z))       # Standard deviations
  z / (sigma %o% sigma)
}

Teraz możemy zrobić PCA:

gw.pca <- function(x, xy, y, tau) {
  # x is a vector denoting a location
  # xy is a set of locations as row vectors
  # y is an array of attributes, also as rows
  # tau is a bandwidth
  # Returns a `princomp` object for the geographically weighted PCA
  # ..of y relative to the point x.
  w <- distance.weight(x, xy, tau)
  princomp(covmat=correlation(y, w))
}

(Do tej pory jest to 10 wierszy kodu wykonywalnego netto. Tylko jeden dodatkowy będzie potrzebny, poniżej, po opisaniu siatki, na której należy przeprowadzić analizę.)

Zilustrujmy przykładowymi danymi losowymi porównywalnymi do tych opisanych w pytaniu: 30 zmiennych w 550 lokalizacjach.

set.seed(17)
n.data <- 550
n.vars <- 30
xy <- matrix(rnorm(n.data * 2), ncol=2)
y <- matrix(rnorm(n.data * n.vars), ncol=n.vars)

Obliczenia ważone geograficznie są często wykonywane na wybranym zestawie lokalizacji, na przykład wzdłuż transektu lub w punktach regularnej siatki. Użyjmy grubej siatki, aby uzyskać pewne spojrzenie na wyniki; później - gdy będziemy pewni, że wszystko działa i otrzymujemy to, czego chcemy - możemy udoskonalić siatkę.

# Create a grid for the GWPCA, sweeping in rows
# from top to bottom.
xmin <- min(xy[,1]); xmax <- max(xy[,1]); n.cols <- 30
ymin <- min(xy[,2]); ymax <- max(xy[,2]); n.rows <- 20
dx <- seq(from=xmin, to=xmax, length.out=n.cols)
dy <- seq(from=ymin, to=ymax, length.out=n.rows)
points <- cbind(rep(dx, length(dy)),
                as.vector(sapply(rev(dy), function(u) rep(u, length(dx)))))

Jest pytanie, jakie informacje chcemy zachować z każdego PCA. Zazwyczaj PCA dla n zmiennych zwraca posortowaną listę n wartości własnych i - w różnych formach - odpowiednią listę n wektorów, każdy o długości n . To n * (n + 1) liczb do zmapowania! Biorąc kilka wskazówek z pytania, zmapujmy wartości własne. Są one wyodrębniane z danych wyjściowych gw.pcaza pomocą $sdevatrybutu, który jest listą wartości własnych według wartości malejącej.

# Illustrate GWPCA by obtaining all eigenvalues at each grid point.
system.time(z <- apply(points, 1, function(x) gw.pca(x, xy, y, 1)$sdev))

Na tym komputerze trwa to mniej niż 5 sekund. Zauważ, że w wywołaniu do użyto charakterystycznej odległości (lub „szerokości pasma”) 1 gw.pca.

Reszta to kwestia zmywania. Zmapujmy wyniki przy użyciu rasterbiblioteki. (Zamiast tego można zapisać wyniki w formacie siatki do przetwarzania końcowego za pomocą GIS).

library("raster")
to.raster <- function(u) raster(matrix(u, nrow=n.cols), 
                                xmn=xmin, xmx=xmax, ymn=ymin, ymx=ymax)
maps <- apply(z, 1, to.raster)
par(mfrow=c(2,2))
tmp <- lapply(maps, function(m) {plot(m); points(xy, pch=19)})

Są to pierwsze cztery z 30 map pokazujących cztery największe wartości własne. (Nie ekscytuj się zbytnio ich wielkościami, które przekraczają 1 w każdej lokalizacji. Przypomnij sobie, że dane te zostały wygenerowane całkowicie losowo, a zatem, jeśli mają w ogóle jakąkolwiek strukturę korelacji - co wydają się wskazywać na duże wartości własne na tych mapach - wynika to wyłącznie z przypadku i nie odzwierciedla niczego „rzeczywistego”, które wyjaśnia proces generowania danych).

Zmiana pasma jest pouczająca. Jeśli jest za mały, oprogramowanie będzie narzekać na osobliwości. (Nie wbudowałem żadnego sprawdzania błędów w tej implementacji od podstaw.) Ale zmniejszenie jej z 1 do 1/4 (i użycie tych samych danych jak poprzednio) daje interesujące wyniki:

Zwróć uwagę na tendencję, aby punkty wokół granicy dawały niezwykle duże główne wartości własne (pokazane w zielonych miejscach mapy po lewej stronie u góry), podczas gdy wszystkie inne wartości własne są obniżone w celu kompensacji (pokazane przez jasnoróżowy na pozostałych trzech mapach) . Zjawisko to oraz wiele innych subtelności PCA i ważenia geograficznego trzeba będzie zrozumieć, zanim będzie można mieć nadzieję na wiarygodną interpretację geograficznie ważonej wersji PCA. Są też inne 30 * 30 = 900 wektorów własnych (lub „ładunków”) do rozważenia….

Aktualizacja:

Obecnie dostępny jest specjalistyczny pakiet R na CRAN - GWmodel, który zawiera między innymi ważone geograficznie PCA. Od autora strony :

Nasz nowy pakiet R do modelowania ważonego geograficznie, GWmodel, został niedawno przesłany do CRAN. GWmodel zapewnia szereg podejść do analizy danych ważonych geograficznie w ramach jednego pakietu, które obejmują statystyki opisowe, korelację, regresję, ogólne modele liniowe i analizę głównych składników. Modele regresji obejmują różne dane dla struktur Gaussa, logistyki i Poissona, a także regresję grzbietową do obsługi skorelowanych predyktorów. Nową funkcją tego pakietu jest dostarczenie solidnych wersji każdej techniki - są one odporne na działanie wartości odstających.

Lokalizacje do modelowania mogą być albo w rzutowanym układzie współrzędnych, albo określone za pomocą współrzędnych geograficznych. Mierniki odległości obejmują euklidesowe, taksówki (Manhattan) i Minkowskiego, a także odległości Wielkiego Kręgu dla lokalizacji określonych przez współrzędne szerokości / długości geograficznej. Dostępne są również różne metody automatycznej kalibracji oraz dostępne są pomocne narzędzia do budowania modeli, które pomagają w wyborze alternatywnych predyktorów.

Podane są również przykładowe zestawy danych, które są wykorzystywane w towarzyszącej dokumentacji w ilustracjach użycia różnych technik.

Więcej szczegółów w podglądzie kolejnego artykułu .

Wątpię, czy istnieje rozwiązanie „gotowe do użycia, podłącz swoje dane”. Ale mam wielką nadzieję, że okaże się, że się mylę, ponieważ chciałbym przetestować tę metodę z niektórymi moimi danymi.

Niektóre opcje do rozważenia:

Marí-Dell'Olmo i współpracownicy wykorzystali analizę czynnikową Bayesa do obliczenia wskaźnika deprywacji dla małych obszarów w Hiszpanii:

Analiza czynnikowa Bayesa w celu obliczenia wskaźnika deprywacji i jego niepewności. Marí-Dell'Olmo M, Martínez-Beneito MA, Borrell C, Zurriaga O, Nolasco A, Domínguez-Berjón MF. Epidemiologia . 2011 maja; 22 (3): 356–64.

W artykule podano specyfikację modelu WinBUGS wykonanego z R, który może pomóc Ci zacząć.

Pakiet adegenet R implementujespcafunkcję. Chociaż koncentruje się na danych genetycznych, równie dobrze może być tak blisko rozwiązania twojego problemu, jak to tylko możliwe. Albo bezpośrednio używając tego pakietu / funkcji, albo modyfikując jego kod. Problemdotyczy winiety, która powinna sprawić, że zaczniesz działać.

Naukowcy Strategicznych Badań Klastra wydaje się być aktywnie pracuje nad tym tematem. Zwłaszcza Paul Harris i Chris Brunsdon (tutaj prezentacja natknąłem się). Niedawna publikacja Paula i Urskiej ( pełny tekst ) może być również użytecznym źródłem:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Analiza głównych składników danych przestrzennych: przegląd. Kroniki Stowarzyszenia Amerykańskich Geografów

Dlaczego nie spróbujesz się z nimi skontaktować i zapytać o rozwiązania, z których dokładnie korzystają? Mogą chcieć podzielić się swoją pracą lub skierować cię w dobrym kierunku.

Cheng, Q. (2006) Przestrzenna i przestrzennie ważona analiza głównych elementów do przetwarzania obrazów. IGARSS 2006: 972-975

w gazetach wspomina się o systemie GeoDAS GIS . Może to być kolejny trop.