如何使用ArcGIS，Python和SPSS / R进行地理加权主成分分析？

#1 楼

“地理位置加权的PCA”非常具有描述性：在R中，该程序实际上编写了自己。 （与实际的代码行相比，它需要更多的注释行。）让我们从权重开始，因为这是PCA本身在地理上加权的PCA零件公司的所在地。术语“地理”是指权重取决于基点和数据位置之间的距离。标准-绝不是唯一-加权是高斯函数；也就是说，指数衰减与平方距离成正比。用户需要指定衰减率，或者更直观地说，是要在其上发生固定衰减量的特征距离。

qA12078q或相关矩阵（从协方差派生）。那么，这里是一个以数值稳定方式计算加权协方差的函数。

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

相关性是通过通常的方式得出的，即使用每个变量的计量单位：

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

现在我们可以执行PCA：

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

（到目前为止，总共有10行可执行代码。下面，在我们描述要执行分析的网格之后，仅需要一行。）


让我们用一些随机示例进行说明与问题中描述的数据可比的数据：550个位置的30个变量。

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

地理加权计算通常是在选定的一组位置上执行的，例如沿着横断面或在常规网格的点处。让我们使用一个粗略的网格对结果进行一些透视；稍后-一旦我们确信一切正常，就可以得到想要的-我们可以优化网格。

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)

我们希望从每个PCA中保留哪些信息，这是一个问题。通常，用于n个变量的PCA返回n个特征值的排序列表，并以各种形式返回n个向量的相应列表，每个向量的长度为n。要映射的是n *（n + 1）个数字！从问题中得到一些提示，让我们映射特征值。这些是通过gw.pca属性从$sdev的输出中提取的，该属性是按降序值排列的特征值列表。 。请注意，在对gw.pca的调用中使用了特征距离（或“带宽”）1。


剩下的就是擦了擦。让我们使用raster库映射结果。 （相反，可以将结果以网格格式写出，以便使用GIS进行后期处理。）

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

30张地图的前四张，显示了四个最大的特征值。 （不要对它们的大小感到兴奋，每个地方的大小都超过1。回想一下，这些数据是完全随机生成的，因此，如果它们完全具有任何相关结构，则这些图中的特征值似乎表明了这一点-这完全是出于偶然，并不反映任何解释数据生成过程的“真实”信息。）

更改带宽很有启发性。如果太小，该软件将抱怨奇异之处。 （我没有在此基本实现中进行任何错误检查。）但是将其从1减少到1/4（并使用与以前相同的数据）确实会得到有趣的结果：



请注意，边界周围的点趋向于给出异常大的主特征值（显示在左上图的绿色位置），而所有其他特征值都被压低以进行补偿（在其他三个图中显示为浅粉红色） 。在人们希望可靠地解释PCA的地理加权版本之前，需要先了解这一现象以及PCA和地理加权的许多其他细微之处。然后还有其他30 * 30 = 900个特征向量（或“载荷”）要考虑...。

#2 楼

更新：

CRAN-GWmodel上现在有一个专用的R包，其中包括地理加权PCA以及其他工具。来自作者的网站：


我们用于地理加权建模的新R包GWmodel最近被上传到CRAN。 GWmodel在单个程序包中提供了地理上的加权数据分析方法，这些方法包括描述性统计，相关性，回归，一般线性模型和主成分分析。回归模型包括具有高斯，逻辑和泊松结构的各种数据，以及用于处理相关预测变量的岭回归。该软件包的新功能是提供每种技术的强健的
版本-这些版本可以抵抗
异常值的影响。

建模位置可以在投影中坐标
系统，或使用地理坐标指定。距离度量标准
包括欧几里得，出租车（曼哈顿）和明可夫斯基，以及由纬度/经度
坐标指定的位置的大
圆距离。还提供了各种自动校准方法，
，并且有一些有用的模型构建工具可用来帮助
从替代的预测变量中进行选择。

还提供了示例数据集，并将其用于在
随附文档中说明了各种
技术的用法。


在即将发表的论文预览中有更多详细信息。


我怀疑是否存在“易于使用，插入数据”的解决方案。
但是我非常希望能被证明是错误的，因为我很想用自己的一些数据来测试这种方法。

要考虑的一些选项： >
Marí-Dell'Olmo及其同事使用贝叶斯因子分析来计算西班牙小区域的贫困指数：


贝叶斯因子分析以计算剥夺指数及其不确定性。
Marí-Dell'OlmoM，Martínez-BeneitoMA，Borrell C，Zurriaga O，Nolasco
A，Domínguez-BerjónMF。
流行病学。 2011年5月； 22（3）：356-64。


在本文中，他们提供了从R执行的WinBUGS模型的规范，可能会帮助您入门。


adegenet R软件包实现了spca功能。尽管它着重于遗传数据，但它可能也尽可能地接近您的问题的解决方案。直接使用此程序包/功能或修改其代码。这个问题上有一个小插图，应该可以让您正常运行。


战略研究小组的研究人员似乎正在积极研究该主题。尤其是Paul Harris和Chris Brunsdon（我偶然发现的演示文稿）。保罗和乌尔斯卡（Paul and Urska）的最新出版物（全文）可能也是有用的资源：


DemšarU，Harris P，Brunsdon C，Fotheringham AS，McLoone S（2012）
主要组成部分空间数据分析：概述。
美国地理学家协会纪事


为什么不尝试与他们联系并询问他们到底使用了什么解决方案？他们可能愿意分享他们的工作或为您指明正确的方向。



Cheng，Q.（2006）空间和空间加权主成分
分析用于图像处理。 IGARSS 2006：972-975



论文提到使用GeoDAS GIS系统。可能会成为另一个线索。