„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….
