# How can I conduct Geographically Weighted Principal Component Analysis using ArcGIS, Python and SPSS/R?

I am after a description/methodology for conducting a Geographically Weighted Principal Components Analysis (GWPCA). I am happy using Python for any portion of this and I imagine SPSS or R being used to run the PCA on the geographically weighted variables.

My dataset is composed of roughly 30 independent variables that are measured throughout ~550 census tracts (vector geometry).

I know this is a loaded question. But, as I search and search, there does not seem to be any solutions out there. What I have come across are mathematical equations that explain the fundamental composition of GWPCA (and GWR). What I am after is more applied in a sense, that I am looking for what major steps I need to accomplish to get from raw data to the GWPCA results.

I would like to expand on the first part with this edit due to the comments received below.

I am basing my interest in GWPCA off of the following paper:

Lloyd, C. D., (2010). Analysing population characteristics using geographically weighted principal components analysis: A case study of Northern Ireland in 2001. Computers, Environment and Urban Systems, 34(5), p.389-399.

For those who do not have access to the literature, I have attached screenshots of the particular sections which explain the mathematics below:

Without going into detail (confidentiality), we are attempting to reduce the 30 variables, which we believe are all very good indicators (albeit globally), to the set of components with eigenvalues greater than 1. By computing the geographically weighted components, we attempt to understand the local variances explained by these components.

I think our primary goal will be to proof the concept of GWPCA, that is, show the spatially explicit nature of our data and that we cannot consider all independent variables to be explanatory on a global scale. Rather, the local scale (neighbourhoods) that each component will identify will aid us in understanding the multi-dimensional nature of our data (how variables can be combined with one another to explain a certain neighbourhoods in our study area).

We hope to map the percentage of variance accounted for by each component (separately), to understand the extent of the neighbourhood explained by the component in question (aid us in understanding local spatiality of our components). Perhaps some other mapping examples but none come to mind at the moment.

The mathematics behind the GWPCA is beyond what I understand given my background in geographic analysis and social statistics. The application of the maths is most important, that is, what do I plug in to these variables/formulas.

-
I don't know o an out of the box solution in R, but it should not be too difficult. Please post the relevant math if you want more feedback than : "R can probably do this". – Paul Hiemstra Oct 7 '12 at 11:09
What kinds of results are you looking for? The largest eigenvalues? An estimated number of principal components? The major steps should be clear enough--at a point, pick weights, compute the weighted covariance (or correlation) matrix, obtain the PCA from the SVD of that matrix. Repeat for a bunch of points. Are you looking for details of any of these steps? – whuber Oct 7 '12 at 19:37
Updated my question. Thank you both! – Michael Markieta Oct 8 '12 at 0:06
my pleasure, whuber. to illustrate my point. n.rows=20 n.cols=30 sq=seq(1,600) rast=raster(matrix(sq,nrow=n.rows,byrow=T)) rast2=raster(matrix(sq,nrow=n.cols)) rast2 is flipped. if you look at your maps, you'll see you actually have 20 columns instead of 30 (wide cells on the x-axis, only 20 of them). just wanted to help out. – user13335 Dec 5 '12 at 3:27
You might be interested to know that there is a new improved package of GW methods for R, including GW PCA, due out soon - it was presented at GISRUK 2013 last month. – AnserGIS May 7 '13 at 10:02

"Geographically weighted PCA" is very descriptive: in `R`, the program practically writes itself. (It needs more comment lines than actual lines of code.)

Lets begin with the weights, because this is where geographically weighted PCA parts company from PCA itself. The term "geographical" means the weights depend on distances between a base point and the data locations. The standard--but by no means only--weighting is a Gaussian function; that is, exponential decay with squared distance. The user needs to specify the decay rate or--more intuitively--a characteristic distance over which a fixed amount of decay occurs.

``````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 applies either to a covariance or correlation matrix (which is derived from a covariance). Here, then, is a function to compute weighted covariances in a numerically stable way.

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

The correlation is derived in the usual way, by using the standard deviations for the units of measurement of each variable:

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

Now we can do the 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))
}
``````

(That's a net 10 lines of executable code so far. Only one more will be needed, below, after we describe a grid over which to perform the analysis.)

Let's illustrate with some random sample data comparable to those described in the question: 30 variables at 550 locations.

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

Geographically weighted calculations are often performed on a selected set of locations, such as along a transect or at points of a regular grid. Let's use a coarse grid to get some perspective on the results; later--once we're confident everything is working and we are getting what we want--we can refine the grid.

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

There's a question of what information we wish to retain from each PCA. Typically, a PCA for n variables returns a sorted list of n eigenvalues and--in various forms--a corresponding list of n vectors, each of length n. That's n*(n+1) numbers to map! Taking some cues from the question, let's map the eigenvalues. These are extracted from the output of `gw.pca` via the `\$sdev` attribute, which is the list of eigenvalues by descending value.

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

This completes in less than 5 seconds on this machine. Notice that a characteristic distance (or "bandwidth") of 1 was used in the call to `gw.pca`.

The rest is a matter of mopping up. Let's map the results using the `raster` library. (Instead, one might write the results out in a grid format for post-processing with a 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)})
``````

These are the first four of the 30 maps, showing the four largest eigenvalues. (Don't get too excited by their sizes, which exceed 1 at every location. Recall that these data were generated totally at random and therefore, if they have any correlation structure at all--which the largish eigenvalues in these maps seem to indicate--it is solely due to chance and does not reflect anything "real" that explains the data generation process.)

It's instructive to change the bandwidth. If it's too small, the software will complain about singularities. (I did not build in any error checking in this bare-bones implementation.) But reducing it from 1 to 1/4 (and using the same data as before) does give interesting results:

Note the tendency for the points around the boundary to give unusually large principal eigenvalues (shown in the green locations of the upper left hand map), while all the other eigenvalues are depressed to compensate (shown by the light pink in the other three maps). This phenomenon, and many other subtleties of PCA and geographic weighting, will need to be understood before one can hope reliably to interpret the geographically weighted version of PCA. And then there are the other 30*30 = 900 eigenvectors (or "loadings") to consider... .

-
+1 Really excellent and thorough answer. – om_henners Oct 23 '12 at 22:18
Remarkable as usual @whuber, thank you very very much! – Michael Markieta Oct 24 '12 at 3:24
just wanted to make you aware that in the to.raster function, you need to have matrix(u, nrow=n.rows, byrow=TRUE) instead of matrix(u, nrow=n.cols). – user12782 Nov 14 '12 at 19:07
@cqh Thank you for looking at this code so carefully! You point to a legitimate concern; I recall having to deal with this issue. However, I think the code is correct as it stands. If I had mixed up the row/column ordering, the illustrations would be totally (and obviously) screwed up. (That's why I tested with different row and column counts.) I apologize for the unfortunate expression `nrow=n.cols`, but that is how it worked out (based on how `points` was created) and I did not want to go back and rename everything. – whuber Nov 14 '12 at 21:04

Update:

There is now a specialized R package available on CRAN - GWmodel that includes geographically weighted PCA among other tools. From author's website:

Our new R package for Geographically Weighted Modelling, GWmodel, was recently uploaded to CRAN. GWmodel provides range of Geographically Weighted data analysis approaches within a single package, these include descriptive statistics, correlation, regression, general linear models and principal components analysis. The regression models include various for data with Gaussian, logistic and Poisson structures, as well as ridge regression for dealing with correlated predictors. A new feature of this package is the provision of robust versions of each technique – these are resistant to the effects of outliers.

Locations for modelling can be either in a projected coordinate system, or specified using geographical coordinates. Distance metrics include Euclidean, taxicab (Manhattan) and Minkowski, as well as Great Circle distances for locations specified by latitude/longitude coordinates. Various automatic calibration methods are also provided, and there are some helpful model building tools available to help select from alternative predictors.

Example datasets are also provided, and they are used in the accompanying documentation in illustrations of the use of the various techniques.

More details ina preview of a forthcoming paper.

I doubt if a 'ready to use, plug in your data' solution exist. But I'm very much hoping to be proven wrong since I'd love to test this method with some of my data.

Some options to consider:

Marí-Dell'Olmo and colleagues used Bayesian factor analysis to calculate deprivation index for small areas in Spain:

Bayesian factor analysis to calculate a deprivation index and its uncertainty. Marí-Dell'Olmo M, Martínez-Beneito MA, Borrell C, Zurriaga O, Nolasco A, Domínguez-Berjón MF. Epidemiology. 2011 May;22(3):356-64.

In the article they provide specification for WinBUGS model executed from R that might get you started.

adegenet R package implements `spca` function. Although it focuses on genetic data, it might as well be as close to a solution for your problem as you can get. Either by using this package/function directly, or modifying its code. There is a vignette on the problem that should get you up and running.

Researchers at Strategic Research Cluster seems to be actively working on the topic. Especially Paul Harris and Chris Brunsdon (here presentation I stumbled upon). Paul and Urska's recent publication (full text) might also be useful resource:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Principal components analysis on spatial data: an overview. Annals of the Association of American Geographers

Why don't you try contacting them and asking about what solutions exactly they are using? They might be willing to share their work or point you in a good direction.

Cheng, Q. (2006) Spatial and Spatially Weighted Principal Component Analysis for Images Processing. IGARSS 2006: 972-975

paper mentions using GeoDAS GIS system. Might be another lead.

-
+1 Excellent links to information about GWPCA. – om_henners Oct 23 '12 at 3:59
+1 The Brunsdon presentation emphasizes the use of PCA as an exploratory tool to find local multivariate outliers. (This use is also featured in the `spca` vignette.) That's a powerful and legitimate use for GWPCA. (However, this method could be much improved, and be more in the spirit of exploratory spatial data analysis, if PCA were replaced by a more robust procedure.) – whuber Oct 25 '12 at 10:19
It seems like an alternative would be kernel PCA. tribesandclimatechange.org/docs/tribes_450.pdf – Jeffrey Evans Dec 5 '12 at 4:45
Thanks for the updated information--`GWmodel` looks like a package worth acquiring. – whuber Jul 19 '13 at 15:44