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

32

Estoy buscando una descripción / metodología para realizar un Análisis de componentes principales (GWPCA). Estoy contento de usar Python para cualquier parte de esto y me imagino que SPSS o R se utilizan para ejecutar el PCA en las variables geográficamente 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 es más aplicado en cierto sentido, que estoy buscando los pasos principales que debo realizar para obtener los datos sin procesar de los resultados de GWPCA.

Me gustaría expandir 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, C. D., (2010). Análisis de las características de la población mediante el análisis de los 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, adjunto capturas de pantalla de las secciones particulares que explican las matemáticas a continuación:

Y para abordar whuber ...

Sin entrar en detalles (confidencialidad), intentamos reducir las 30 variables, que creemos que son todos indicadores muy buenos (aunque a nivel mundial), al conjunto de componentes con valores propios mayores que 1. Al calcular los componentes geográficamente ponderados , 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 cada componente identificará nos ayudará a comprender la naturaleza multidimensional de nuestros datos (cómo se pueden combinar las variables entre sí para explicar ciertos vecindarios en nuestra área de estudio).

Esperamos hacer un mapa del porcentaje de varianza contabilizada por cada componente (por separado), para comprender el alcance del vecindario explicado por el componente en cuestión (nos ayuda a comprender la espacialidad local de nuestros componentes). Tal vez algunos otros ejemplos de mapeo, pero ninguno viene a la mente en este momento.

Además:

Las matemáticas detrás del GWPCA están más allá de lo que entiendo, dada 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é debo conectar a estas variables / fórmulas?

    
pregunta Michael Markieta 07.10.2012 - 07:49

2 respuestas

29

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

Comencemos con los pesos, porque aquí es donde la compañía de piezas de PCA geográficamente ponderada de PCA. El término "geográfico" significa que los pesos dependen de las distancias entre un punto base y las ubicaciones de los datos. La ponderación estándar, pero de ninguna manera solo, es una función gaussiana; es decir, decaimiento exponencial con distancia al cuadrado. El usuario debe especificar la tasa de caída o, de manera más intuitiva, una distancia característica sobre la cual se produce una cantidad fija de caída.

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 covarianza o matriz de correlación (que se deriva de una covarianza). Aquí, entonces, hay una función para calcular las 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))
}

(Eso 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 que realizar el análisis.)

Ilustremos 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 se realizan a menudo en un conjunto seleccionado de ubicaciones, como a lo largo de un transecto o en puntos de una cuadrícula regular. Usemos una cuadrícula gruesa para obtener una perspectiva de los resultados; más tarde, una vez que estemos seguros de que todo está funcionando y obtenemos lo que queremos, podemos refinar la cuadrícula.

# 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 queremos retener de cada PCA. Normalmente, una PCA para las variables n devuelve una lista ordenada de valores propios n y, en varias formas, una lista correspondiente de vectores n , cada uno de longitud n . ¡Eso es n * (n + 1) números para mapear! Tomando algunas pistas de la pregunta, vamos a mapear los valores propios. Estos se extraen de la salida de gw.pca a través del atributo $sdev , 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 cuestión de limpiar. Asignemos los resultados utilizando la biblioteca raster . (En su lugar, uno podría escribir los resultados en un formato de cuadrícula para su posterior procesamiento 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)})

Estossonlosprimeroscuatrodelos30mapas,quemuestranloscuatrovalorespropiosmásgrandes.(Noseentusiasmedemasiadoconsustamaños,queexcedende1encadaubicación.Recuerdequeestosdatossegenerarontotalmentealazary,porlotanto,sitienenalgunaestructuradecorrelación,loqueLosvalorespropiosenestosmapasparecenindicarquesedebeúnicamentealazarynoreflejanada"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 incorporé ninguna comprobación de errores en esta implementación básica). Pero reducirla de 1 a 1/4 (y usar los mismos datos que antes) da resultados interesantes:

Observelatendenciadelospuntosalrededordellímiteadarvalorespropiosprincipalesinusualmentegrandes(quesemuestranenlasubicacionesverdesdelmapadelapartesuperiorizquierda),mientrasquetodoslosdemásvalorespropiossepresionanparacompensar(semuestraconlaluzrosaenlaotratresmapas).Estefenómeno,ymuchasotrassutilezasdelaPCAylaponderacióngeográfica,deberánentenderseantesdequeunopuedaesperardemaneraconfiableinterpretarlaversióngeográficamenteponderadadelaPCA.Yluegoestánlosotros30*30=900vectorespropios(o"cargas") para considerar ...

    
respondido por el whuber 23.10.2012 - 22:48
14

Actualizar:

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

  

Nuestro nuevo paquete R para modelado ponderado geográficamente, GWmodel, fue   recientemente subido a CRAN. GWmodel proporciona gama de Geográficamente   Los enfoques de análisis de datos ponderados dentro de un solo paquete, estos   Incluye estadística descriptiva, correlación, regresión, general.   Análisis de modelos lineales y componentes principales. La regresión   Los modelos incluyen varios para datos con gaussiano, logístico y Poisson.   estructuras, así como la regresión de cresta para tratar con correlacionado   predictores Una nueva característica de este paquete es la provisión de robustos   Versiones de cada técnica - son resistentes a los efectos de   valores atípicos.

     

Las ubicaciones para el modelado pueden estar en una coordenada proyectada   Sistema, o especificado utilizando coordenadas geográficas. Métricas de distancia   incluyen Euclidean, taxicab (Manhattan) y Minkowski, así como Great   Distancias circulares para lugares especificados por latitud / longitud   coordenadas 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   seleccione de predictores alternativos.

     

También se proporcionan ejemplos de conjuntos de datos, y se utilizan en el   Documentación adjunta en ilustraciones del uso de los diversos   técnicas.

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

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

Algunas opciones a considerar:

Marí-Dell'Olmo y sus 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 una especificación para el modelo WinBUGS ejecutado desde R que podría ayudarte a comenzar.

adegenet el paquete R implementa la función spca . Aunque se centra en los datos genéticos, podría ser lo más cercano a una solución para su problema. Ya sea utilizando este paquete / función directamente, o modificando su código. Hay una vignette sobre el problema que debe ponerlo en funcionamiento .

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

  

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

¿Por qué no intentas contactarte con ellos y preguntarles qué soluciones están usando exactamente? Podrían estar dispuestos a compartir su trabajo o señalarle una buena dirección.

  

Cheng, Q. (2006) Componente principal espacial y espacialmente ponderado   Análisis para el procesamiento de imágenes. IGARSS 2006: 972-975

menciones en papel usando GeoDAS GIS sistema. Podría ser otra ventaja.

    
respondido por el radek 11.10.2012 - 14:33

Lea otras preguntas en las etiquetas