Wie kann ich mit ArcGIS, Python und SPSS / R eine geografisch gewichtete Hauptkomponentenanalyse durchführen?

Ich bin nach einer Beschreibung / Methodik für die Durchführung einer geografisch gewichteten Hauptkomponentenanalyse (GWPCA). Ich bin glücklich, Python für einen Teil davon zu verwenden, und ich stelle mir vor, dass SPSS oder R verwendet werden, um die PCA für die geografisch gewichteten Variablen auszuführen.

Mein Datensatz besteht aus ungefähr 30 unabhängigen Variablen, die in ~ 550 Zensus-Traktaten (Vektorgeometrie) gemessen werden.

Ich weiß, dass dies eine geladene Frage ist. Aber während ich suche und suche, scheint es keine Lösungen zu geben. Was mir begegnet ist, sind mathematische Gleichungen, die die grundlegende Zusammensetzung von GWPCA (und GWR) erklären. Was ich anstrebe, ist eher in dem Sinne anzuwenden, dass ich nach den wichtigsten Schritten suche, die ich ausführen muss, um von den Rohdaten zu den GWPCA-Ergebnissen zu gelangen.

Aufgrund der folgenden Kommentare möchte ich den ersten Teil mit dieser Bearbeitung erweitern.

Um Paul anzusprechen ...

Ich begründe mein Interesse an GWPCA mit dem folgenden Artikel:

Lloyd, CD, (2010). Analyse der Bevölkerungsmerkmale mit Hilfe der geografisch gewichteten Hauptkomponentenanalyse: Eine Fallstudie von Nordirland aus dem Jahr 2001. Computer, Umwelt und städtische Systeme, 34 (5), S. 389-399.

Für diejenigen, die keinen Zugang zur Literatur haben, habe ich Screenshots der einzelnen Abschnitte beigefügt, die die Mathematik unten erläutern:

Und um whuber anzusprechen ...

Ohne ins Detail zu gehen (Vertraulichkeit), versuchen wir, die 30 Variablen, von denen wir glauben, dass sie alle sehr gute Indikatoren sind (wenn auch global), auf die Menge der Komponenten mit Eigenwerten größer als 1 zu reduzieren um die lokalen Abweichungen zu verstehen, die durch diese Komponenten erklärt werden.

Meines Erachtens wird es unser Hauptziel sein, das Konzept von GWPCA zu beweisen, dh die räumliche Eindeutigkeit unserer Daten aufzuzeigen, und dass wir nicht alle unabhängigen Variablen als global erklärend betrachten können. Vielmehr hilft uns die lokale Skala (Nachbarschaften), die jede Komponente identifiziert, beim Verständnis der mehrdimensionalen Natur unserer Daten (wie Variablen miteinander kombiniert werden können, um bestimmte Nachbarschaften in unserem Untersuchungsgebiet zu erklären).

Wir hoffen, den prozentualen Anteil der Varianz, der auf jede Komponente entfällt (separat), abbilden zu können, um das Ausmaß der durch die betreffende Komponente erklärten Nachbarschaft zu verstehen (um das Verständnis der lokalen Räumlichkeit unserer Komponenten zu erleichtern). Vielleicht kommen Ihnen einige andere Mapping-Beispiele in den Sinn, aber derzeit keine.

Zusätzlich:

Die Mathematik hinter dem GWPCA geht über das hinaus, was ich aufgrund meines Hintergrunds in der geografischen Analyse und der Sozialstatistik verstehe. Die Anwendung der Mathematik ist am wichtigsten, das heißt, was verbinde ich mit diesen Variablen / Formeln?

"Geographisch gewichtetes PCA" ist sehr beschreibend: In Rschreibt sich das Programm praktisch von selbst. (Es benötigt mehr Kommentarzeilen als tatsächliche Codezeilen.)

Beginnen wir mit den Gewichten, denn hier befindet sich die geografisch gewichtete PCA-Teilefirma von PCA selbst. Der Begriff "geografisch" bedeutet, dass die Gewichte von den Entfernungen zwischen einem Basispunkt und den Datenpositionen abhängen. Die Standardgewichtung - aber keineswegs nur - ist eine Gaußsche Funktion; exponentieller Zerfall mit quadratischer Entfernung. Der Benutzer muss die Abklingrate oder - intuitiver - eine charakteristische Entfernung angeben, über die eine feste Abklingmenge auftritt.

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 gilt entweder für eine Kovarianz oder eine Korrelationsmatrix (die von einer Kovarianz abgeleitet ist). Hier ist also eine Funktion, um gewichtete Kovarianzen numerisch stabil zu berechnen.

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))  
}

Die Korrelation wird auf übliche Weise unter Verwendung der Standardabweichungen für die Maßeinheiten jeder Variablen abgeleitet:

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

Jetzt können wir die PCA machen:

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))
}

(Das sind bisher netto 10 Zeilen ausführbaren Codes. Nach der Beschreibung eines Rasters, über das die Analyse durchgeführt werden soll, wird nur noch eine Zeile benötigt.)

Veranschaulichen wir anhand einiger Zufallsdaten, die mit den in der Frage beschriebenen Daten vergleichbar sind: 30 Variablen an 550 Standorten.

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)

Geografisch gewichtete Berechnungen werden häufig an ausgewählten Orten durchgeführt, z. B. entlang eines Schnitts oder an Punkten eines regulären Gitters. Lassen Sie uns ein grobes Raster verwenden, um einen Überblick über die Ergebnisse zu erhalten. Später - wenn wir sicher sind, dass alles funktioniert und wir bekommen, was wir wollen - können wir das Raster verfeinern.

# 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)))))

Es ist eine Frage, welche Informationen wir von jedem PCA behalten möchten. Typischerweise gibt eine PCA für n Variablen eine sortierte Liste von n Eigenwerten und - in verschiedenen Formen - eine entsprechende Liste von n Vektoren mit jeweils der Länge n zurück . Das sind n * (n + 1) Zahlen für die Karte! Nehmen wir einige Hinweise aus der Frage und ordnen Sie die Eigenwerte zu. Diese werden aus der Ausgabe von gw.pcaüber das $sdevAttribut extrahiert , bei dem es sich um die Liste der Eigenwerte nach absteigendem Wert handelt.

# 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))

Dies ist auf diesem Computer in weniger als 5 Sekunden erledigt. Beachten Sie, dass beim Aufruf von eine charakteristische Entfernung (oder "Bandbreite") von 1 verwendet wurde gw.pca.

Der Rest ist eine Frage des Aufwischens. Lassen Sie uns die Ergebnisse mithilfe der rasterBibliothek zuordnen. (Stattdessen könnte man die Ergebnisse in einem Rasterformat für die Nachbearbeitung mit einem GIS ausgeben.)

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)})

Dies sind die ersten vier der 30 Karten, die die vier größten Eigenwerte zeigen. (Seien Sie nicht zu aufgeregt über ihre Größe, die an jedem Ort 1 übersteigt. Denken Sie daran, dass diese Daten völlig zufällig generiert wurden und daher, wenn sie überhaupt eine Korrelationsstruktur haben - was die größeren Eigenwerte in diesen Karten anzudeuten scheinen - es ist ausschließlich zufällig und spiegelt nichts "Reales" wider, was den Prozess der Datengenerierung erklärt.)

Es ist aufschlussreich, die Bandbreite zu ändern. Wenn es zu klein ist, beschwert sich die Software über Singularitäten. (Ich habe in dieser Bare-Bones-Implementierung keine Fehlerprüfung eingebaut.) Eine Reduzierung von 1 auf 1/4 (und die Verwendung der gleichen Daten wie zuvor) führt jedoch zu interessanten Ergebnissen:

Beachten Sie die Tendenz, dass die Punkte um die Grenze ungewöhnlich große Haupteigenwerte ergeben (angezeigt an den grünen Stellen der oberen linken Karte), während alle anderen Eigenwerte zur Kompensation gedrückt werden (angezeigt durch das Hellrosa in den anderen drei Karten). . Dieses Phänomen und viele andere Feinheiten der PCA und der geografischen Gewichtung müssen verstanden werden, bevor die geografisch gewichtete Version der PCA zuverlässig interpretiert werden kann. Und dann gibt es noch die anderen 30 * 30 = 900 Eigenvektoren (oder "Ladungen") zu berücksichtigen ....

Aktualisieren:

Für CRAN- GWmodel ist jetzt ein spezielles R-Paket verfügbar , das unter anderem geografisch gewichtete PCA enthält. Von der Website des Autors :

Unser neues R-Paket für geografisch gewichtete Modellierung, GWmodel, wurde kürzlich zu CRAN hochgeladen. GWmodel bietet eine Reihe von Ansätzen zur geografisch gewichteten Datenanalyse in einem einzigen Paket, darunter deskriptive Statistiken, Korrelation, Regression, allgemeine lineare Modelle und Hauptkomponentenanalyse. Die Regressionsmodelle umfassen verschiedene Daten für Gaußsche, logistische und Poisson-Strukturen sowie Kammregressionen für den Umgang mit korrelierten Prädiktoren. Ein neues Merkmal dieses Pakets ist die Bereitstellung robuster Versionen jeder Technik - diese sind resistent gegen die Auswirkungen von Ausreißern.

Orte für die Modellierung können sich entweder in einem projizierten Koordinatensystem befinden oder mithilfe von geografischen Koordinaten angegeben werden. Zu den Entfernungsmesswerten gehören Euklid, Taxis (Manhattan) und Minkowski sowie Großkreisentfernungen für Standorte, die durch Breiten- / Längengradkoordinaten angegeben werden. Verschiedene automatische Kalibrierungsmethoden sind ebenfalls verfügbar, und es stehen einige hilfreiche Werkzeuge zur Modellbildung zur Verfügung, mit deren Hilfe Sie alternative Prädiktoren auswählen können.

Es werden auch Beispieldatensätze bereitgestellt, die in der Begleitdokumentation zur Veranschaulichung der Verwendung der verschiedenen Techniken verwendet werden.

Weitere Details in einer Vorschau auf eine bevorstehende Veröffentlichung .

Ich bezweifle, dass es eine Lösung für die sofortige Verwendung Ihrer Daten gibt. Aber ich hoffe sehr, dass ich mich als falsch erweisen kann, da ich diese Methode gerne mit einigen meiner Daten testen würde.

Einige zu berücksichtigende Optionen:

Marí-Dell'Olmo und Kollegen verwendeten die Bayes-Faktor-Analyse, um den Deprivation-Index für kleine Gebiete in Spanien zu berechnen:

Bayes-Faktor-Analyse zur Berechnung eines Deprivationsindex und seiner Unsicherheit. Marí-Dell'Olmo M., Martínez-Beneito M. A., Borrell C., Zurriaga O., Nolasco A., Domínguez-Berjón MF. Epidemiologie . 2011 May; 22 (3): 356 & ndash; 64.

In dem Artikel wird eine Spezifikation für das WinBUGS-Modell angegeben, das von R ausgeführt wird und möglicherweise den Einstieg erleichtert.

Das adegenet R-Paket implementiert diespcaFunktion. Obwohl es sich auf genetische Daten konzentriert, könnte es genauso gut einer Lösung für Ihr Problem nahe kommen, wie Sie sie bekommen können. Entweder durch direkte Verwendung dieses Pakets / dieser Funktion oder durch Ändern des Codes. Es gibt eine Vignette zum Problem, die Sie zum Laufen bringen soll.

Forscher des Strategic Research Cluster scheinen aktiv an dem Thema zu arbeiten. Vor allem Paul Harris und Chris Brunsdon (hier Präsentation stieß ich auf). Die jüngste Veröffentlichung von Paul und Urska ( Volltext ) könnte auch eine nützliche Ressource sein:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Analyse der Hauptkomponenten für räumliche Daten: ein Überblick. Annalen der Association of American Geographers

Warum versuchst du nicht, sie zu kontaktieren und zu fragen, welche Lösungen genau sie verwenden? Sie könnten bereit sein, ihre Arbeit zu teilen oder Sie in eine gute Richtung zu weisen.

Cheng, Q. (2006) Räumliche und räumlich gewichtete Hauptkomponentenanalyse für die Bildverarbeitung. IGARSS 2006: 972 & ndash; 975

Papier erwähnt mit GeoDAS GIS- System. Könnte eine andere Spur sein.