¿Cómo puedo realizar un análisis de componentes principales ponderado geográficamente utilizando ArcGIS, Python y SPSS / R?

Busco una descripción / metodología para realizar un Análisis de componentes principales ponderados geográficamente (GWPCA). Estoy feliz de usar Python para cualquier parte de esto e imagino que SPSS o R se utilizan para ejecutar el PCA en las variables geográficas ponderadas.

Mi conjunto de datos se compone de aproximadamente 30 variables independientes que se miden a lo largo de ~ 550 secciones censales (geometría vectorial).

Sé que esta es una pregunta cargada. Pero, mientras busco y busco, no parece haber ninguna solución por ahí. Lo que he encontrado son ecuaciones matemáticas que explican la composición fundamental de GWPCA (y GWR). Lo que busco se aplica más en cierto sentido, es decir, estoy buscando los pasos principales que debo realizar para obtener datos sin procesar para obtener los resultados de GWPCA.

Me gustaría ampliar la primera parte con esta edición debido a los comentarios recibidos a continuación.

Para dirigirse a Paul ...

Estoy basando mi interés en GWPCA en el siguiente documento:

Lloyd, CD, (2010). Análisis de las características de la población mediante el análisis de componentes principales ponderados geográficamente: un estudio de caso de Irlanda del Norte en 2001. Computadoras, medio ambiente y sistemas urbanos, 34 (5), p.389-399.

Para aquellos que no tienen acceso a la literatura, he adjuntado capturas de pantalla de las secciones particulares que explican las matemáticas a continuación:

Y para dirigirse a whuber ...

Sin entrar en detalles (confidencialidad), estamos tratando de reducir las 30 variables, que creemos que son todos muy buenos indicadores (aunque globalmente), al conjunto de componentes con valores propios mayores que 1. Al calcular los componentes ponderados geográficamente, intentamos comprender las variaciones locales explicadas por estos componentes.

Creo que nuestro objetivo principal será probar el concepto de GWPCA, es decir, mostrar la naturaleza espacialmente explícita de nuestros datos y que no podemos considerar que todas las variables independientes sean explicativas a escala global. Más bien, la escala local (vecindarios) que identificará cada componente nos ayudará a comprender la naturaleza multidimensional de nuestros datos (cómo las variables se pueden combinar entre sí para explicar ciertos vecindarios en nuestra área de estudio).

Esperamos mapear el porcentaje de varianza explicado por cada componente (por separado), para comprender el alcance del vecindario explicado por el componente en cuestión (ayúdenos a comprender la espacialidad local de nuestros componentes). Quizás otros ejemplos de mapeo, pero ninguno viene a la mente en este momento.

Adicionalmente:

Las matemáticas detrás del GWPCA están más allá de lo que entiendo dado mi experiencia en análisis geográfico y estadísticas sociales. La aplicación de las matemáticas es lo más importante, es decir, qué conecto a estas variables / fórmulas.

"PCA ponderado geográficamente" es muy descriptivo: en R, el programa prácticamente se escribe solo. (Necesita más líneas de comentarios que líneas de código reales).

Comencemos con los pesos, porque aquí es donde la empresa de piezas PCA ponderada geográficamente de la propia PCA. El término "geográfico" significa que los pesos dependen de las distancias entre un punto base y las ubicaciones de datos. La ponderación estándar, pero de ninguna manera solamente, es una función gaussiana; es decir, disminución exponencial con distancia al cuadrado. El usuario necesita especificar la tasa de descomposición o, más intuitivamente, una distancia característica sobre la cual ocurre una cantidad fija de descomposición.

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 se aplica a una matriz de covarianza o correlación (que se deriva de una covarianza). Aquí, entonces, es una función para calcular covarianzas ponderadas de una manera numéricamente estable.

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 correlación se deriva de la forma habitual, utilizando las desviaciones estándar para las unidades de medida de cada variable:

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

Ahora podemos hacer el 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))
}

(Esa es una red de 10 líneas de código ejecutable hasta ahora. Solo se necesitará una más, a continuación, después de que describamos una cuadrícula sobre la cual realizar el análisis).

Vamos a ilustrar con algunos datos de muestra aleatorios comparables a los descritos en la pregunta: 30 variables en 550 ubicaciones.

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)

Los cálculos ponderados geográficamente a menudo se realizan en un conjunto seleccionado de ubicaciones, como a lo largo de un transecto o en puntos de una cuadrícula regular. Usemos una grilla gruesa para obtener una perspectiva de los resultados; más tarde, una vez que estamos seguros de que todo está funcionando y estamos obteniendo lo que queremos, podemos refinar la red.

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

Hay una pregunta sobre qué información deseamos retener de cada PCA. Típicamente, un PCA para n variables devuelve una lista ordenada de n valores propios y, en varias formas, una lista correspondiente de n vectores, cada uno de longitud n . ¡Eso es n * (n + 1) números para mapear! Tomando algunas pistas de la pregunta, mapeemos los valores propios. Estos se extraen de la salida de a gw.pcatravés del $sdevatributo, que es la lista de valores propios por valor descendente.

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

Esto se completa en menos de 5 segundos en esta máquina. Observe que se utilizó una distancia característica (o "ancho de banda") de 1 en la llamada a gw.pca.

El resto es una cuestión de limpieza. Vamos a mapear los resultados usando la rasterbiblioteca. (En cambio, uno podría escribir los resultados en un formato de cuadrícula para el procesamiento posterior con 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)})

Estos son los primeros cuatro de los 30 mapas, que muestran los cuatro valores propios más grandes. (No se entusiasme demasiado con sus tamaños, que exceden 1 en cada ubicación. Recuerde que estos datos se generaron totalmente al azar y, por lo tanto, si tienen alguna estructura de correlación, lo que parecen indicar los valores propios más grandes en estos mapas) - se debe únicamente al azar y no refleja nada "real" que explique el proceso de generación de datos).

Es instructivo cambiar el ancho de banda. Si es demasiado pequeño, el software se quejará de las singularidades. (No incluí ninguna comprobación de errores en esta implementación básica). Pero reducirlo de 1 a 1/4 (y usar los mismos datos que antes) da resultados interesantes:

Tenga en cuenta la tendencia de los puntos alrededor del límite a dar valores propios principales inusualmente grandes (que se muestran en las ubicaciones verdes del mapa superior izquierdo), mientras que todos los otros valores propios se deprimen para compensar (como se muestra en rosa claro en los otros tres mapas) . Este fenómeno, y muchas otras sutilezas de PCA y ponderación geográfica, deberán entenderse antes de que uno pueda esperar interpretar de manera confiable la versión ponderada geográficamente de PCA. Y luego están los otros 30 * 30 = 900 vectores propios (o "cargas") a considerar ....

Actualizar:

Ahora hay un paquete R especializado disponible en CRAN - GWmodel que incluye PCA ponderada geográficamente entre otras herramientas. Del sitio web del autor :

Nuestro nuevo paquete R para modelado ponderado geográficamente, GWmodel, se subió recientemente a CRAN. GWmodel proporciona una gama de enfoques de análisis de datos ponderados geográficamente dentro de un solo paquete, que incluyen estadísticas descriptivas, correlación, regresión, modelos lineales generales y análisis de componentes principales. Los modelos de regresión incluyen varios para datos con estructuras gaussianas, logísticas y de Poisson, así como regresión de cresta para tratar con predictores correlacionados. Una nueva característica de este paquete es la provisión de versiones robustas de cada técnica, que son resistentes a los efectos atípicos.

Las ubicaciones para el modelado pueden ser en un sistema de coordenadas proyectadas o especificarse usando coordenadas geográficas. Las métricas de distancia incluyen Euclidean, taxi (Manhattan) y Minkowski, así como distancias del Gran Círculo para ubicaciones especificadas por coordenadas de latitud / longitud. También se proporcionan varios métodos de calibración automática, y hay algunas herramientas útiles de construcción de modelos disponibles para ayudar a seleccionar entre predictores alternativos.

También se proporcionan conjuntos de datos de ejemplo, y se usan en la documentación adjunta en ilustraciones del uso de las diversas técnicas.

Más detalles en una vista previa de un próximo trabajo .

Dudo que exista una solución 'lista para usar, conecte sus datos'. Pero espero que se demuestre que estoy equivocado, ya que me encantaría probar este método con algunos de mis datos.

Algunas opciones a considerar:

Marí-Dell'Olmo y colegas utilizaron el análisis factorial bayesiano para calcular el índice de privación para áreas pequeñas en España:

Análisis factorial bayesiano para calcular un índice de privación y su incertidumbre. Marí-Dell'Olmo M, Martínez-Beneito MA, Borrell C, Zurriaga O, Nolasco A, Domínguez-Berjón MF. Epidemiología . Mayo de 2011; 22 (3): 356-64.

En el artículo, proporcionan especificaciones para el modelo WinBUGS ejecutado desde R que puede ayudarlo a comenzar.

El paquete adegenet R implementa laspcafunción. Aunque se centra en los datos genéticos, podría ser lo más parecido a una solución para su problema que pueda obtener. Ya sea usando este paquete / función directamente, o modificando su código. Hay una viñeta sobre el problema que debería ponerlo en funcionamiento.

Los investigadores del Cluster de Investigación Estratégica parecen estar trabajando activamente en el tema. Especialmente Paul Harris y Chris Brunsdon (aquí me encontré con la presentación ). La publicación reciente de Paul y Urska ( texto completo ) también podría ser un recurso útil:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Análisis de componentes principales en datos espaciales: una visión general. Anales de la Asociación de Geógrafos Americanos

¿Por qué no intentas contactarlos y preguntarles qué soluciones están usando exactamente? Pueden estar dispuestos a compartir su trabajo o señalarle en una buena dirección.

Cheng, Q. (2006) Análisis de componentes principales espaciales y espacialmente ponderados para el procesamiento de imágenes. IGARSS 2006: 972-975

El documento menciona el uso del sistema GeoDAS GIS . Podría ser otra pista.