Comment puis-je effectuer une analyse en composantes principales pondérée géographiquement à l'aide d'ArcGIS, de Python et de SPSS / R?

Je suis après une description / méthodologie pour mener une analyse en composantes principales pondérée géographiquement (GWPCA). Je suis heureux d’utiliser Python pour n’importe quelle partie de cela et j’imagine que SPSS ou R sont utilisés pour exécuter la PCA sur les variables pondérées géographiquement.

Mon jeu de données est composé d'environ 30 variables indépendantes qui sont mesurées dans environ 550 secteurs de recensement (géométrie vectorielle).

Je sais que c'est une question chargée. Mais, lorsque je cherche et que je cherche, il ne semble pas y avoir de solutions. Ce que j'ai rencontré sont des équations mathématiques qui expliquent la composition fondamentale de GWPCA (et GWR). Ce que je recherche est plus appliqué dans un sens, c’est-à-dire que je recherche les étapes principales que je dois accomplir pour passer des données brutes aux résultats de la GWPCA.

Je voudrais développer la première partie avec cette modification en raison des commentaires reçus ci-dessous.

Pour adresser Paul ...

Je fonde mon intérêt pour GWPCA sur le document suivant:

Lloyd, CD, (2010). Analyse des caractéristiques de la population à l'aide d'une analyse en composantes principales pondérée géographiquement: une étude de cas de l'Irlande du Nord en 2001. Computers, Environment and Urban Systems, 34 (5), p. 389-399.

Pour ceux qui n'ont pas accès à la littérature, j'ai joint des captures d'écran des sections particulières expliquant les mathématiques ci-dessous:

Et pour aborder whuber ...

Sans entrer dans les détails (confidentialité), nous essayons de réduire les 30 variables, que nous considérons tous comme de très bons indicateurs (bien que globalement), à l'ensemble des composants avec des valeurs propres supérieures à 1. En calculant les composants pondérés géographiquement, comprendre les variances locales expliquées par ces composantes.

Notre objectif principal sera de prouver le concept de GWPCA, c'est-à-dire de montrer la nature spatialement explicite de nos données et que nous ne pouvons pas considérer toutes les variables indépendantes comme explicatives à l'échelle mondiale. Plutôt, l'échelle locale (les quartiers) que chaque composante identifiera nous aidera à comprendre la nature multidimensionnelle de nos données (comment les variables peuvent être combinées entre elles pour expliquer certains quartiers de notre zone d'étude).

Nous espérons cartographier le pourcentage de variance représenté par chaque composant (séparément), pour comprendre l'étendue du voisinage expliquée par le composant en question (nous aider à comprendre la spatialité locale de nos composants). Peut-être quelques autres exemples de cartographie, mais aucun ne vient à l'esprit pour le moment.

Aditionellement:

Les mathématiques derrière la GWPCA vont au-delà de ce que je comprends compte tenu de ma formation en analyse géographique et en statistiques sociales. L’application des mathématiques est la plus importante, c’est-à-dire, que dois-je connecter à ces variables / formules?

"PCA géographiquement pondéré" est très descriptif: Rle programme s’écrit pratiquement lui-même. (Il faut plus de lignes de commentaires que de lignes de code.)

Commençons par les poids, car c’est là que PCA, géographiquement pondérée, se sépare de la société PCA. Le terme "géographique" signifie que les poids dépendent des distances entre un point de base et les emplacements de données. La pondération standard - mais pas seulement - est une fonction gaussienne; c'est-à-dire une décroissance exponentielle avec une distance au carré. L'utilisateur doit spécifier le taux de décroissance ou, plus intuitivement, une distance caractéristique sur laquelle une décroissance fixe se produit.

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

La CPA s’applique soit à une covariance, soit à une matrice de corrélation (dérivée d’une covariance). Voici donc une fonction permettant de calculer les covariances pondérées de manière numériquement stable.

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

La corrélation est dérivée de la manière habituelle, en utilisant les écarts-types pour les unités de mesure de chaque variable:

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

Maintenant nous pouvons faire la 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))
}

(Cela représente jusqu'à présent 10 lignes de code exécutable. Il ne vous en faudra plus qu'une, après avoir décrit une grille sur laquelle effectuer l'analyse.)

Illustrons avec des échantillons de données aléatoires comparables à ceux décrits dans la question: 30 variables sur 550 emplacements.

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)

Les calculs pondérés géographiquement sont souvent effectués sur un ensemble sélectionné d'emplacements, par exemple le long d'un transect ou à des points d'une grille régulière. Utilisons une grille grossière pour avoir une perspective sur les résultats; plus tard - une fois que nous sommes certains que tout fonctionne et que nous obtenons ce que nous voulons - nous pouvons affiner le réseau.

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

Il y a une question de savoir quelles informations nous souhaitons conserver de chaque APC. En règle générale, une ACP pour n variables renvoie une liste triée de n valeurs propres et - sous diverses formes - une liste correspondante de n vecteurs, chacun de longueur n . C'est n * (n + 1) nombres à mapper! En prenant quelques indices de la question, cartographions les valeurs propres. Celles-ci sont extraites de la sortie de gw.pcavia l' $sdevattribut, qui correspond à la liste des valeurs propres par valeur décroissante.

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

Cela se termine en moins de 5 secondes sur cette machine. Notez qu'une distance caractéristique (ou "bande passante") de 1 a été utilisée dans l'appel à gw.pca.

Le reste est une question de nettoyage. Mappons les résultats à l'aide de la rasterbibliothèque. (Au lieu de cela, on pourrait écrire les résultats dans un format de grille pour le post-traitement avec un SIG.)

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

Ce sont les quatre premières des 30 cartes, montrant les quatre plus grandes valeurs propres. (Ne soyez pas trop excités par leurs tailles, qui dépassent 1 à chaque endroit. Rappelez-vous que ces données ont été générées totalement au hasard et donc, si elles ont une structure de corrélation - ce que les valeurs propres approximatives dans ces cartes semblent indiquer. - C’est uniquement dû au hasard et ne reflète en rien ce qui est "réel" qui explique le processus de génération de données.)

Il est instructif de changer la bande passante. Si c'est trop petit, le logiciel se plaindra des singularités. (Je n'ai intégré aucune vérification d'erreur dans cette implémentation complète.) Mais le réduire de 1 à 1/4 (et en utilisant les mêmes données qu'auparavant) donne des résultats intéressants:

Notez la tendance des points autour de la limite à donner des valeurs propres principales inhabituellement grandes (indiquées dans les emplacements verts de la carte en haut à gauche), tandis que toutes les autres valeurs propres sont déprimées pour compenser (indiquées par le rose pâle dans les trois autres cartes). . Ce phénomène, ainsi que de nombreuses autres subtilités de la PCA et de la pondération géographique, devra être compris avant de pouvoir espérer interpréter de manière fiable la version pondérée de la PCA. Et puis il y a les 30 * 30 = 900 autres vecteurs propres (ou "charges") à considérer ....

Mise à jour:

CRAN - GWmodel propose désormais un package R spécialisé qui comprend, entre autres outils, une analyse PCA pondérée géographiquement. Du site de l'auteur :

Notre nouveau package R pour la modélisation pondérée géographiquement, GWmodel, a récemment été chargé sur CRAN. GWmodel propose une gamme d'approches d'analyse de données pondérées géographiquement au sein d'un même progiciel, notamment les statistiques descriptives, la corrélation, la régression, les modèles linéaires généraux et l'analyse en composantes principales. Les modèles de régression comprennent divers types de données avec des structures gaussiennes, logistiques et de Poisson, ainsi qu'une régression de crête permettant de traiter des prédicteurs corrélés. Une nouvelle fonctionnalité de ce package est la fourniture de versions robustes de chaque technique - celles-ci résistent aux effets des valeurs aberrantes.

Les emplacements pour la modélisation peuvent être soit dans un système de coordonnées projetées, soit spécifiés à l'aide de coordonnées géographiques. Les mesures de distance incluent Euclidean, taxicab (Manhattan) et Minkowski, ainsi que les distances du grand cercle pour les emplacements spécifiés par les coordonnées de latitude / longitude. Diverses méthodes d'étalonnage automatique sont également fournies, ainsi que des outils de création de modèle utiles permettant de sélectionner d'autres prédicteurs.

Des exemples de jeux de données sont également fournis et sont utilisés dans la documentation d'accompagnement pour illustrer l'utilisation des différentes techniques.

Plus de détails dans l'aperçu d'un prochain article .

Je doute qu'il existe une solution «prête à l'emploi, connectez vos données». Mais j'espère beaucoup avoir tort, car j'aimerais bien tester cette méthode avec certaines de mes données.

Quelques options à considérer:

Marí-Dell'Olmo et ses collègues ont utilisé une analyse factorielle bayésienne pour calculer l'indice de défavorisation pour les petites régions d'Espagne:

Analyse factorielle bayésienne pour calculer un indice de défavorisation et son incertitude. Marí-Dell'Olmo M, MA Martínez-Beneito, C Borrell, O Zurriaga, A Nolasco, MF Domínguez-Berjón. Épidémiologie . 2011 mai; 22 (3): 356-64.

Dans cet article, ils fournissent des spécifications du modèle WinBUGS exécuté à partir de R qui pourraient vous aider à démarrer.

Le package adegenet R implémente laspcafonction. Bien qu'il se concentre sur les données génétiques, il pourrait tout aussi bien être aussi proche d'une solution à votre problème que vous pouvez obtenir. Soit en utilisant directement ce paquet / cette fonction, soit en modifiant son code. Il existe une vignette sur le problème qui devrait vous permettre de démarrer.

Les chercheurs du groupe de recherche stratégique semble travailler activement sur le sujet. Surtout Paul Harris et Chris Brunsdon (la présentation que je suis tombé sur). La récente publication de Paul et Urska ( texte intégral ) pourrait également être une ressource utile:

Demšar U, Harris P, C Brunsdon, AS Fotheringham, McLoone S (2012) Analyse en composantes principales des données spatiales: un aperçu. Annales de l'Association of American Geographers

Pourquoi n'essayez-vous pas de les contacter et de leur demander quelles solutions ils utilisent exactement? Ils pourraient être disposés à partager leur travail ou à vous orienter dans la bonne direction.

Cheng, Q. (2006), Analyse en composantes principales pondérée spatialement pour le traitement d'images. IGARSS 2006: 972-975

L'article mentionne l'utilisation du système GeoDAS GIS . Peut-être une autre piste.