Como posso conduzir uma análise de componentes principais ponderada geograficamente usando ArcGIS, Python e SPSS / R?

Estou buscando uma descrição / metodologia para realizar uma Análise de Componentes Principais Ponderada Geograficamente (GWPCA). Estou feliz em usar o Python para qualquer parte disso e imagino que o SPSS ou R esteja sendo usado para executar o PCA nas variáveis ​​geograficamente ponderadas.

Meu conjunto de dados é composto por aproximadamente 30 variáveis ​​independentes que são medidas ao longo de ~ 550 setores censitários (geometria vetorial).

Eu sei que esta é uma pergunta carregada. Mas, enquanto eu procuro e busco, não parece haver nenhuma solução por aí. O que encontrei são equações matemáticas que explicam a composição fundamental do GWPCA (e GWR). O que eu busco é mais aplicado, em certo sentido, que estou procurando as principais etapas que preciso realizar para obter dados brutos e resultados da GWPCA.

Gostaria de expandir a primeira parte desta edição devido aos comentários recebidos abaixo.

Para abordar Paul ...

Estou baseando meu interesse na GWPCA no seguinte artigo:

Lloyd, CD, (2010). Analisando características da população usando análise de componentes principais ponderados geograficamente: Um estudo de caso da Irlanda do Norte em 2001. Computers, Environment and Urban Systems, 34 (5), p.389-399.

Para aqueles que não têm acesso à literatura, anexei capturas de tela das seções específicas que explicam a matemática abaixo:

E para abordar whuber ...

Sem entrar em detalhes (confidencialidade), estamos tentando reduzir as 30 variáveis, que acreditamos serem todos indicadores muito bons (embora globalmente), para o conjunto de componentes com autovalores maiores que 1. Ao calcular os componentes geograficamente ponderados, tentamos para entender as variações locais explicadas por esses componentes.

Acho que nosso objetivo principal será provar o conceito de GWPCA, ou seja, mostrar a natureza espacialmente explícita de nossos dados e que não podemos considerar todas as variáveis ​​independentes como explicativas em escala global. Em vez disso, a escala local (bairros) que cada componente identificará nos ajudará a entender a natureza multidimensional de nossos dados (como as variáveis ​​podem ser combinadas entre si para explicar determinados bairros em nossa área de estudo).

Esperamos mapear a porcentagem de variação contabilizada por cada componente (separadamente), para entender a extensão da vizinhança explicada pelo componente em questão (nos ajuda a entender a espacialidade local de nossos componentes). Talvez alguns outros exemplos de mapeamento, mas nenhum venha à mente no momento.

Além disso:

A matemática por trás da GWPCA está além do que eu entendo, dada minha formação em análise geográfica e estatística social. A aplicação da matemática é mais importante, ou seja, o que eu insiro nessas variáveis ​​/ fórmulas.

"PCA geograficamente ponderado" é muito descritivo: em R, o programa praticamente se escreve. (Ele precisa de mais linhas de comentário do que linhas de código reais.)

Vamos começar com os pesos, porque é aqui que a empresa de peças de PCA geograficamente ponderada da própria PCA. O termo "geográfico" significa que os pesos dependem das distâncias entre um ponto base e os locais dos dados. A ponderação padrão - mas de nenhuma maneira somente - é uma função gaussiana; isto é, decaimento exponencial com distância ao quadrado. O usuário precisa especificar a taxa de decaimento ou - mais intuitivamente - uma distância característica sobre a qual ocorre uma quantidade fixa de decaimento.

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

O PCA se aplica a uma matriz de covariância ou correlação (que é derivada de uma covariância). Aqui, então, é uma função para calcular covariâncias ponderadas de maneira numericamente estável.

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

A correlação é derivada da maneira usual, usando os desvios padrão para as unidades de medida de cada variável:

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

Agora podemos fazer o 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))
}

(Até agora, são 10 linhas líquidas de código executável. Somente mais uma será necessária, abaixo, depois de descrevermos uma grade sobre a qual executar a análise.)

Vamos ilustrar com alguns dados de amostra aleatória comparáveis ​​aos descritos na pergunta: 30 variáveis ​​em 550 locais.

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)

Os cálculos geograficamente ponderados são geralmente realizados em um conjunto selecionado de locais, como ao longo de uma seção transversal ou em pontos de uma grade regular. Vamos usar uma grade grossa para ter uma perspectiva dos resultados; mais tarde - quando estivermos confiantes de que tudo está funcionando e conseguimos o que queremos - podemos refinar a grade.

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

Há uma pergunta de quais informações queremos reter de cada PCA. Normalmente, um PCA para n variáveis ​​retorna uma lista classificada de n autovalores e - de várias formas - uma lista correspondente de n vetores, cada um com comprimento n . São n * (n + 1) números para mapear! Tomando algumas dicas da pergunta, vamos mapear os autovalores. Eles são extraídos da saída de gw.pcavia $sdevatributo, que é a lista de valores próprios por valor decrescente.

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

Isso é concluído em menos de 5 segundos nesta máquina. Observe que uma distância característica (ou "largura de banda") de 1 foi usada na chamada para gw.pca.

O resto é uma questão de limpar. Vamos mapear os resultados usando a rasterbiblioteca. (Em vez disso, pode-se escrever os resultados em um formato de grade para pós-processamento com um 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)})

Estes são os quatro primeiros dos 30 mapas, mostrando os quatro maiores valores próprios. (Não fique muito empolgado com seus tamanhos, que excedem 1 em cada local. Lembre-se de que esses dados foram gerados totalmente aleatoriamente e, portanto, se eles possuem alguma estrutura de correlação - os quais os autovalores amplos nesses mapas parecem indicar - é apenas devido ao acaso e não reflete nada "real" que explique o processo de geração de dados.)

É instrutivo alterar a largura de banda. Se for muito pequeno, o software irá reclamar de singularidades. (Não construí nenhuma verificação de erro nesta implementação básica). Mas reduzi-lo de 1 para 1/4 (e usar os mesmos dados de antes) fornece resultados interessantes:

Observe a tendência dos pontos ao redor do limite de fornecer autovalores principais invulgarmente grandes (mostrados nas localizações verdes do mapa superior esquerdo), enquanto todos os outros autovalores são pressionados para compensar (mostrados em rosa claro nos outros três mapas) . Esse fenômeno, e muitas outras sutilezas do PCA e da ponderação geográfica, precisarão ser entendidos antes que se possa esperar com segurança interpretar a versão ponderada geograficamente do PCA. E existem os outros 30 * 30 = 900 autovetores (ou "cargas") a serem considerados ....

Atualizar:

Agora existe um pacote R especializado disponível no modelo CRAN - GW que inclui PCA geograficamente ponderado, entre outras ferramentas. Do site do autor :

Nosso novo pacote R para modelagem geograficamente ponderada, GWmodel, foi recentemente carregado no CRAN. O GWmodel fornece diversas abordagens de análise de dados com ponderação geográfica em um único pacote, incluindo estatísticas descritivas, correlação, regressão, modelos lineares gerais e análise de componentes principais. Os modelos de regressão incluem vários dados para estruturas gaussianas, logísticas e de Poisson, bem como regressão de cume para lidar com preditores correlacionados. Um novo recurso deste pacote é o fornecimento de versões robustas de cada técnica - elas são resistentes aos efeitos de outliers.

Os locais para modelagem podem estar em um sistema de coordenadas projetadas ou especificados usando coordenadas geográficas. As métricas de distância incluem Euclidean, táxi (Manhattan) e Minkowski, bem como distâncias do Grande Círculo para locais especificados pelas coordenadas de latitude / longitude. Vários métodos de calibração automática também são fornecidos, e existem algumas ferramentas úteis de criação de modelos disponíveis para ajudar a selecionar preditores alternativos.

Conjuntos de dados de exemplo também são fornecidos e são usados ​​na documentação anexa nas ilustrações do uso das várias técnicas.

Mais detalhes em uma prévia de um próximo artigo .

Duvido que exista uma solução 'pronta para usar, conecte seus dados'. Mas espero muito estar errado, porque eu adoraria testar esse método com alguns dos meus dados.

Algumas opções a serem consideradas:

Marí-Dell'Olmo e colegas usaram a análise fatorial bayesiana para calcular o índice de privação de pequenas áreas na Espanha:

Análise fatorial bayesiana para calcular um índice de privação e sua incerteza. Marí-Dell'Olmo M, Martínez-Beneito MA, Borrell C, Zurriaga O, Nolasco A, Domínguez-Berjón MF. Epidemiologia . 2011 maio; 22 (3): 356-64.

No artigo, eles fornecem especificações para o modelo WinBUGS executado a partir do R que pode ajudar você a começar.

O pacote adegenet R implementa aspcafunção. Embora se concentre nos dados genéticos, pode ser o mais próximo possível de uma solução para o seu problema. Ou usando este pacote / função diretamente ou modificando seu código. Há uma vinheta sobre o problema que deve colocá-lo em funcionamento.

Pesquisadores do Strategic Research Cluster parecem estar trabalhando ativamente no assunto. Especialmente Paul Harris e Chris Brunsdon (aqui apresentação eu tropecei em cima). A publicação recente de Paul e Urska ( texto completo ) também pode ser um recurso útil:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Análise de componentes principais em dados espaciais: uma visão geral. Anais da Associação de Geógrafos Americanos

Por que você não tenta contatá-los e perguntar sobre quais soluções exatamente eles estão usando? Eles podem estar dispostos a compartilhar seu trabalho ou apontá-lo em uma boa direção.

Cheng, Q. (2006) Análise de Componentes Principais Espacial e Espacialmente Ponderada para Processamento de Imagens. IGARSS 2006: 972-975

O artigo menciona o uso do sistema GeoDAS GIS . Pode ser outra pista.