Как я могу проводить анализ главных компонентов с географической взвешенностью, используя ArcGIS, Python и SPSS / R?

Я после описания / методологии для проведения анализа географически взвешенных основных компонентов (GWPCA). Я счастлив использовать Python для любой части этого, и я представляю, что SPSS или R используются для запуска PCA с географически взвешенными переменными.

Мой набор данных состоит из примерно 30 независимых переменных, которые измеряются по ~ 550 участкам переписи (векторная геометрия).

Я знаю, что это загруженный вопрос. Но, как я ищу, так и не вижу никаких решений. Я столкнулся с математическими уравнениями, которые объясняют фундаментальный состав GWPCA (и GWR). То, что мне нужно, в большей степени применимо в том смысле, что я ищу, какие основные шаги мне нужно выполнить, чтобы перейти от необработанных данных к результатам GWPCA.

Я хотел бы расширить первую часть этого редактирования из-за комментариев, полученных ниже.

Чтобы обратиться к Павлу ...

Я основываю свой интерес к GWPCA из следующей статьи:

Ллойд, CD, (2010). Анализ характеристик населения с использованием географически взвешенного анализа основных компонентов: тематическое исследование Северной Ирландии в 2001 году. Компьютеры, окружающая среда и городские системы, 34 (5), с.389-399.

Для тех, кто не имеет доступа к литературе, я приложил скриншоты отдельных разделов, которые объясняют математику ниже:

И по адресу ...

Не вдаваясь в детали (конфиденциальность), мы пытаемся свести 30 переменных, которые, по нашему мнению, являются очень хорошими показателями (хотя и глобально), к набору компонентов с собственными значениями, превышающими 1. Путем вычисления географически взвешенных компонентов мы пытаемся чтобы понять местные различия, объясняемые этими компонентами.

Я думаю, что наша основная цель будет состоять в том, чтобы доказать концепцию GWPCA, то есть показать пространственно явный характер наших данных и что мы не можем считать все независимые переменные пояснительными в глобальном масштабе. Скорее, локальный масштаб (окрестности), который будет идентифицировать каждый компонент, поможет нам понять многомерный характер наших данных (как переменные могут быть объединены друг с другом, чтобы объяснить определенные окрестности в нашей области исследования).

Мы надеемся отобразить процент дисперсии, учитываемый каждым компонентом (отдельно), чтобы понять степень соседства, объясняемого рассматриваемым компонентом (помочь нам в понимании локальной пространственности наших компонентов). Возможно, некоторые другие примеры картирования, но ни один из них не приходит на ум в данный момент.

Дополнительно:

Математика, лежащая в основе GWPCA, выходит за рамки того, что я понимаю, учитывая мои знания в области географического анализа и социальной статистики. Применение математики является наиболее важным, то есть то, что я подключить к этим переменным / формулам.

«Географически взвешенный PCA» очень нагляден: Rпрограмма практически сама себя пишет. (Требуется больше строк комментариев, чем реальных строк кода.)

Начнем с весов, потому что именно здесь географически взвешенная компания PCA разделяет саму PCA. Термин «географический» означает, что веса зависят от расстояний между базовой точкой и местоположениями данных. Стандарт - но не только - взвешивание - это гауссовская функция; то есть экспоненциальный спад с квадратом расстояния. Пользователь должен указать скорость затухания или - более интуитивно - характерное расстояние, на котором происходит фиксированное количество затухания.

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 применяется либо к ковариационной, либо к корреляционной матрице (которая выводится из ковариации). Здесь, затем, функция для вычисления взвешенных ковариаций численно устойчивым способом.

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

Корреляция получается обычным способом, используя стандартные отклонения для единиц измерения каждой переменной:

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

Теперь мы можем сделать СПС:

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

(На данный момент это всего 10 строк исполняемого кода. Ниже потребуется еще одна, после того, как мы опишем сетку, по которой будет выполняться анализ.)

Давайте проиллюстрируем некоторые случайные данные выборки, сопоставимые с описанными в вопросе: 30 переменных в 550 местах.

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)

Географически взвешенные вычисления часто выполняются в выбранном наборе мест, например, вдоль разреза или в точках регулярной сетки. Давайте использовать грубую сетку, чтобы получить представление о результатах; позже - когда мы уверены, что все работает и мы получаем то, что хотим, - мы можем улучшить сетку.

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

Возникает вопрос, какую информацию мы хотим сохранить от каждого СПС. Как правило, PCA для n переменных возвращает отсортированный список из n собственных значений и - в различных формах - соответствующий список из n векторов, каждый из которых имеет длину n . Это n * (n + 1) номеров на карте! Если взять некоторые подсказки из этого вопроса, давайте сопоставим собственные значения. Они извлекается из выходного сигнала gw.pcaчерез $sdevатрибут, который является списком собственных значений по убыванию значения.

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

Это завершается менее чем за 5 секунд на этой машине. Обратите внимание, что характерное расстояние (или «полоса пропускания»), равное 1, было использовано при вызове gw.pca.

Остальное дело зачистки. Давайте сопоставим результаты, используя rasterбиблиотеку. (Вместо этого можно записать результаты в формате сетки для последующей обработки с помощью ГИС.)

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

Это первые четыре из 30 карт, показывающие четыре самых больших собственных значения. (Не будьте слишком взволнованы их размерами, которые превышают 1 в каждом месте. Напомним, что эти данные были сгенерированы совершенно случайным образом и, следовательно, если они вообще имеют какую-либо структуру корреляции - что, по-видимому, указывают на большие собственные значения в этих картах - это исключительно случайно и не отражает ничего «реального», объясняющего процесс генерации данных.)

Поучительно изменить пропускную способность. Если он слишком маленький, программное обеспечение будет жаловаться на особенности. (Я не включал проверку ошибок в этой простой реализации.) Но уменьшение ее с 1 до 1/4 (и использование тех же данных, что и раньше) действительно дает интересные результаты:

Обратите внимание на тенденцию для точек вокруг границы давать необычно большие главные собственные значения (показанные в зеленых местах на верхней левой карте), в то время как все другие собственные значения подавлены, чтобы компенсировать (показанный светло-розовым на других трех картах) , Это явление и многие другие тонкости PCA и географического взвешивания необходимо будет понять, прежде чем можно будет надежно интерпретировать географически взвешенную версию PCA. И затем есть другие 30 * 30 = 900 собственных векторов (или «нагрузок»), которые необходимо учитывать ...

Обновить:

В настоящее время в CRAN доступен специальный пакет R - GWmodel, который включает в себя географически взвешенный PCA и другие инструменты. С сайта автора :

Наш новый пакет R для географически взвешенного моделирования, GWmodel, был недавно загружен в CRAN. GWmodel предоставляет ряд подходов к анализу географически взвешенных данных в рамках одного пакета, которые включают описательную статистику, корреляцию, регрессию, общие линейные модели и анализ главных компонентов. Модели регрессии включают в себя различные данные для гауссовых, логистических и пуассоновских структур, а также регрессию гребня для работы с коррелированными предикторами. Новая особенность этого пакета - предоставление надежных версий каждой техники - они устойчивы к воздействию выбросов.

Места для моделирования могут быть либо в проекционной системе координат, либо заданы с использованием географических координат. Метрики расстояний включают евклидову, такси (Манхэттен) и Минковского, а также расстояния по Большому кругу для местоположений, указанных координатами широты / долготы. Также предусмотрены различные методы автоматической калибровки, и есть несколько полезных инструментов построения моделей, которые помогут выбрать один из альтернативных предикторов.

Примеры наборов данных также предоставляются, и они используются в сопроводительной документации в иллюстрациях использования различных методов.

Больше подробностей в превью предстоящей статьи .

Я сомневаюсь, что существует решение «готово к использованию, подключи свои данные». Но я очень надеюсь, что окажется неправым, так как хотел бы проверить этот метод на некоторых моих данных.

Некоторые варианты для рассмотрения:

Мари-Делл'Ольмо и его коллеги использовали байесовский факторный анализ для расчета индекса депривации для небольших территорий в Испании:

Байесовский факторный анализ для расчета индекса депривации и его неопределенности. Мари-Делл'Ольмо М, Мартинес-Бенейто М.А., Боррелл С, Сурриага О, Ноласко А, Домингес-Берхон М.Ф. Эпидемиология . 2011 май; 22 (3): 356-64.

В этой статье они предоставляют спецификацию для модели WinBUGS, выполненной из R, которая может помочь вам начать работу.

Пакет adegenet R реализуетspcaфункцию. Хотя он сфокусирован на генетических данных, он может быть как можно ближе к решению вашей проблемы. Либо используя этот пакет / функцию напрямую, либо изменив его код. Существует виньетка на проблемекоторая должна получить вас и работает.

Исследователи из Кластера стратегических исследований, похоже, активно работают над этой темой. Особенно Пол Харрис и Крис Брансдон (здесь я наткнулся на презентацию ). Недавняя публикация Пола и Урски ( полный текст ) также может быть полезным ресурсом:

Демшар У, Харрис П., Брунсдон С., Фотерингем А.С., Маклоун С. (2012) Анализ основных компонентов пространственных данных: обзор. Летопись Ассоциации американских географов

Почему бы вам не попробовать связаться с ними и спросить, какие именно решения они используют? Они могут быть готовы поделиться своей работой или указать вам правильное направление.

Cheng, Q. (2006) Пространственный и пространственно-взвешенный анализ главных компонентов для обработки изображений. IGARSS 2006: 972-975

в статье упоминается использование ГИС системы GeoDAS . Может быть, еще один повод.