"PCA ponderado geográficamente" es muy descriptivo: en R El programa se escribe prácticamente solo. (Necesita más líneas de comentario que líneas de código reales).
Empecemos por las ponderaciones, porque es aquí donde el PCA ponderado geográficamente se separa del propio PCA. El término "geográfico" significa que las ponderaciones dependen de las distancias entre un punto base y las ubicaciones de los datos. La ponderación estándar, pero no la única, es una función gaussiana, es decir, un decaimiento exponencial con la distancia al cuadrado. El usuario tiene que especificar la tasa de decaimiento o -más intuitivamente- una distancia característica sobre la que se produce una cantidad fija de decaimiento.
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)))
}
El ACP se aplica a una matriz de covarianza o de correlación (que se deriva de una covarianza). He aquí, pues, una función para calcular covarianzas ponderadas de forma 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 obtiene de la forma habitual, utilizando las desviaciones estándar de 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 un neto de 10 líneas de código ejecutable hasta ahora. Sólo se necesitará una más, más adelante, después de describir una cuadrícula sobre la que realizar el análisis).
Ilustrémoslo con algunos datos de muestras aleatorias comparables a las descritas en la pregunta: 30 variables en 550 lugares.
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 suelen realizarse en un conjunto seleccionado de lugares, como a lo largo de un transecto o en puntos de una cuadrícula regular. Utilicemos una cuadrícula gruesa para tener una perspectiva de los resultados; más adelante -cuando estemos seguros de que todo funciona y obtenemos lo que queremos- podremos afinar 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)))))
Se trata de saber qué información queremos conservar de cada PCA. Normalmente, una PCA para n devuelve una lista ordenada de n valores propios y -en varias formas- una lista correspondiente de n vectores, cada uno de ellos de longitud n . Son n*(n+1) números que hay que 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 de la $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 ha utilizado 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 utilizando el raster biblioteca. (En su lugar, se podrían escribir los resultados en un formato de cuadrícula para su postprocesamiento 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)})
![Maps]()
Estos son los primeros cuatro de los 30 mapas, que muestran los cuatro mayores valores propios. (No te emociones demasiado por sus tamaños, que superan el 1 en cada lugar. Recordemos que estos datos fueron generados totalmente al azar y, por tanto, si tienen alguna estructura de correlación -lo que parecen indicar los grandes valores propios de 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 he incorporado ninguna comprobación de errores en esta implementación básica.) Pero reducirlo de 1 a 1/4 (y utilizar los mismos datos que antes) da resultados interesantes:
![Maps 2]()
Obsérvese la tendencia de los puntos alrededor de la frontera a dar valores propios principales inusualmente grandes (mostrados en las ubicaciones verdes del mapa superior izquierdo), mientras que todos los demás valores propios se reducen para compensar (mostrados por el rosa claro en los otros tres mapas). Este fenómeno, así como muchas otras sutilezas del ACP y de la ponderación geográfica, deberán comprenderse antes de poder esperar interpretar de forma fiable la versión del ACP ponderada geográficamente. Además, hay que tener en cuenta los otros 30*30 = 900 vectores propios (o "cargas")...
