Bagaimana saya bisa melakukan Analisis Komponen Utama yang Tertimbang Secara Geografis menggunakan ArcGIS, Python, dan SPSS / R?

Saya mencari penjelasan / metodologi untuk melakukan Analisis Komponen Utama yang Tertimbang Secara Geografis (GWPCA). Saya senang menggunakan Python untuk setiap bagian dari ini dan saya membayangkan SPSS atau R digunakan untuk menjalankan PCA pada variabel-variabel geografis tertimbang.

Dataset saya terdiri dari sekitar 30 variabel independen yang diukur di seluruh ~ 550 trus sensus (vektor geometri).

Saya tahu ini adalah pertanyaan yang dimuat. Tapi, saat saya mencari dan mencari, sepertinya tidak ada solusi di luar sana. Apa yang saya temui adalah persamaan matematika yang menjelaskan komposisi dasar GWPCA (dan GWR). Apa yang saya kejar lebih banyak diterapkan dalam arti, bahwa saya mencari langkah-langkah utama apa yang perlu saya selesaikan untuk mendapatkan dari data mentah ke hasil GWPCA.

Saya ingin memperluas pada bagian pertama dengan edit ini karena komentar yang diterima di bawah ini.

Untuk mengatasi Paul ...

Saya mendasarkan minat saya pada GWPCA dari makalah berikut:

Lloyd, CD, (2010). Menganalisa karakteristik populasi menggunakan analisis komponen utama yang berbobot geografis: Studi kasus Irlandia Utara pada tahun 2001. Komputer, Lingkungan dan Sistem Perkotaan, 34 (5), hal.389-399.

Bagi mereka yang tidak memiliki akses ke literatur, saya telah melampirkan tangkapan layar bagian-bagian tertentu yang menjelaskan matematika di bawah ini:

Dan untuk mengatasi ...

Tanpa merinci (kerahasiaan), kami berusaha mengurangi 30 variabel, yang kami yakini semua indikator yang sangat baik (meskipun secara global), ke sekumpulan komponen dengan nilai eigen lebih besar dari 1. Dengan menghitung komponen yang ditimbang secara geografis, kami berupaya untuk memahami varian lokal yang dijelaskan oleh komponen-komponen ini.

Saya pikir tujuan utama kami adalah untuk membuktikan konsep GWPCA, yaitu menunjukkan sifat eksplisit spasial dari data kami dan bahwa kami tidak dapat menganggap semua variabel independen sebagai penjelas pada skala global. Sebaliknya, skala lokal (lingkungan) yang akan diidentifikasi oleh masing-masing komponen akan membantu kita dalam memahami sifat multi-dimensi data kami (bagaimana variabel dapat digabungkan satu sama lain untuk menjelaskan lingkungan tertentu di wilayah studi kami).

Kami berharap dapat memetakan persentase varians yang diperhitungkan oleh masing-masing komponen (secara terpisah), untuk memahami tingkat lingkungan yang dijelaskan oleh komponen yang dimaksud (membantu kami dalam memahami spasial lokal dari komponen kami). Mungkin beberapa contoh pemetaan lain tetapi tidak ada yang terlintas dalam pikiran saat ini.

Selain itu:

Matematika di balik GWPCA melampaui apa yang saya pahami mengingat latar belakang saya dalam analisis geografis dan statistik sosial. Penerapan matematika adalah yang paling penting, yaitu, apa yang saya hubungkan ke variabel / formula ini.

"PCA tertimbang secara geografis" sangat deskriptif: dalam R, program ini praktis menulis sendiri. (Perlu lebih banyak baris komentar daripada baris kode yang sebenarnya.)

Mari kita mulai dengan bobot, karena ini adalah di mana perusahaan suku cadang PCA secara geografis tertimbang dari PCA itu sendiri. Istilah "geografis" berarti bobot tergantung pada jarak antara titik dasar dan lokasi data. Standar - tetapi tidak berarti hanya - pembobotan adalah fungsi Gaussian; yaitu, peluruhan eksponensial dengan jarak kuadrat. Pengguna perlu menentukan tingkat peluruhan atau - lebih intuitif - jarak karakteristik di mana jumlah peluruhan tetap terjadi.

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 berlaku untuk matriks kovarians atau korelasi (yang berasal dari kovarians). Di sini, kemudian, adalah fungsi untuk menghitung kovarian tertimbang dengan cara yang stabil secara numerik.

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

Korelasi diturunkan dengan cara biasa, dengan menggunakan standar deviasi untuk unit pengukuran setiap variabel:

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

Sekarang kita bisa melakukan 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))
}

(Sejauh ini, 10 baris net dari kode yang dapat dieksekusi. Hanya satu lagi yang akan diperlukan, di bawah ini, setelah kami menggambarkan kisi yang akan digunakan untuk melakukan analisis.)

Mari kita ilustrasikan dengan beberapa data sampel acak yang sebanding dengan yang dijelaskan dalam pertanyaan: 30 variabel di 550 lokasi.

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)

Perhitungan berbobot geografis sering dilakukan pada set lokasi yang dipilih, seperti sepanjang transek atau pada titik-titik grid biasa. Mari kita gunakan kisi kasar untuk mendapatkan perspektif tentang hasilnya; nanti - setelah kami yakin semuanya bekerja dan kami mendapatkan apa yang kami inginkan - kami dapat memperbaiki 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)))))

Ada pertanyaan tentang informasi apa yang ingin kami simpan dari setiap PCA. Biasanya, PCA untuk n variabel return daftar diurutkan n nilai eigen dan - dalam berbagai bentuk - daftar yang sesuai dari n vektor, masing-masing dengan panjang n . Itu n * (n +1) angka untuk dipetakan! Dengan mengambil beberapa isyarat dari pertanyaan, mari kita petakan nilai eigen. Ini diekstraksi dari output gw.pcamelalui $sdevatribut, yang merupakan daftar nilai eigen dengan nilai menurun.

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

Ini selesai dalam waktu kurang dari 5 detik pada mesin ini. Perhatikan bahwa jarak karakteristik (atau "bandwidth") 1 digunakan dalam panggilan ke gw.pca.

Sisanya adalah masalah pembersihan. Mari kita petakan hasilnya menggunakan rasterperpustakaan. (Sebagai gantinya, orang mungkin menuliskan hasilnya dalam format kisi untuk pasca-pemrosesan dengan 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)})

Ini adalah empat dari 30 peta pertama, yang menunjukkan empat nilai eigen terbesar. (Jangan terlalu senang dengan ukurannya, yang melebihi 1 di setiap lokasi. Ingat bahwa data ini dihasilkan secara acak dan oleh karena itu, jika mereka memiliki struktur korelasi sama sekali - yang nilai eigen besar dalam peta ini tampaknya mengindikasikan - itu semata-mata karena kebetulan dan tidak mencerminkan sesuatu yang "nyata" yang menjelaskan proses pembuatan data.)

Ini instruktif untuk mengubah bandwidth. Jika terlalu kecil, perangkat lunak akan mengeluh tentang singularitas. (Saya tidak membuat kesalahan saat memeriksa implementasi kosong ini.) Tetapi menguranginya dari 1 menjadi 1/4 (dan menggunakan data yang sama seperti sebelumnya) memang memberikan hasil yang menarik:

Perhatikan kecenderungan titik-titik di sekitar batas untuk memberikan nilai eigen utama yang luar biasa besar (ditunjukkan di lokasi hijau peta kiri atas), sementara semua nilai eigen lainnya ditekan untuk mengkompensasi (ditunjukkan oleh warna merah muda terang di tiga peta lainnya) . Fenomena ini, dan banyak seluk-beluk lainnya dari PCA dan pembobotan geografis, perlu dipahami sebelum orang dapat berharap untuk menafsirkan versi PCA yang ditimbang secara geografis. Dan kemudian ada 30 * 30 = 900 vektor eigen lainnya (atau "memuat") untuk dipertimbangkan ....

Memperbarui:

Sekarang ada paket R khusus yang tersedia di CRAN - GWmodel yang mencakup PCA yang ditimbang secara geografis di antara alat lainnya. Dari situs web penulis :

Paket R baru kami untuk Pemodelan Berbobot Geografis, GWmodel, baru-baru ini diunggah ke CRAN. GWmodel menyediakan berbagai pendekatan analisis data tertimbang secara Geografis dalam satu paket, termasuk statistik deskriptif, korelasi, regresi, model linier umum, dan analisis komponen utama. Model regresi mencakup berbagai data dengan struktur Gaussian, logistik dan Poisson, serta regresi ridge untuk berurusan dengan prediktor berkorelasi. Fitur baru dari paket ini adalah penyediaan versi yang kuat dari setiap teknik - ini tahan terhadap efek outlier.

Lokasi untuk pemodelan dapat berupa sistem koordinat yang diproyeksikan, atau ditentukan menggunakan koordinat geografis. Metrik jarak mencakup Euclidean, taxicab (Manhattan) dan Minkowski, serta jarak Lingkaran Besar untuk lokasi yang ditentukan oleh koordinat lintang / bujur. Berbagai metode kalibrasi otomatis juga disediakan, dan ada beberapa alat pembangun model yang bermanfaat tersedia untuk membantu memilih dari pemrediksi alternatif.

Kumpulan data contoh juga disediakan, dan mereka digunakan dalam dokumentasi yang menyertai dalam ilustrasi penggunaan berbagai teknik.

Lebih detail dalam pratinjau makalah yang akan datang .

Saya ragu apakah ada solusi 'siap pakai, tancapkan data Anda'. Tetapi saya sangat berharap untuk terbukti salah karena saya ingin menguji metode ini dengan beberapa data saya.

Beberapa opsi untuk dipertimbangkan:

Marí-Dell'Olmo dan rekannya menggunakan analisis faktor Bayesian untuk menghitung indeks kekurangan untuk area kecil di Spanyol:

Analisis faktor Bayesian untuk menghitung indeks perampasan dan ketidakpastiannya. Marí-Dell'Olmo M, Martínez-Beneito MA, Borrell C, Zurriaga O, Nolasco A, Domínguez-Berjón MF. Epidemiologi . 2011 Mei; 22 (3): 356-64.

Dalam artikel tersebut mereka memberikan spesifikasi untuk model WinBUGS yang dijalankan dari R yang mungkin bisa Anda mulai.

adegenet R package mengimplementasikanspcafungsi. Meskipun berfokus pada data genetik, itu mungkin juga dekat dengan solusi untuk masalah Anda. Baik dengan menggunakan paket / fungsi ini secara langsung, atau memodifikasi kodenya. Ada sketsa tentang masalah yang seharusnya membuat Anda bangun dan berjalan.

Para peneliti di Cluster Penelitian Strategis tampaknya secara aktif bekerja pada topik tersebut. Terutama Paul Harris dan Chris Brunsdon (di sini presentasi yang saya temukan ). Publikasi terbaru Paul dan Urska ( teks lengkap ) juga dapat menjadi sumber yang berguna:

Demšar U, Harris P, Brunsdon C, Fotheringham AS, McLoone S (2012) Analisis komponen utama pada data spasial: ikhtisar. Annals of Association of American Geographers

Mengapa Anda tidak mencoba menghubungi mereka dan bertanya tentang solusi apa yang sebenarnya mereka gunakan? Mereka mungkin bersedia membagikan pekerjaan mereka atau mengarahkan Anda ke arah yang baik.

Cheng, Q. (2006) Analisis Komponen Utama Spasial dan Tertimbang Secara Spasial untuk Pemrosesan Gambar. IGARSS 2006: 972-975

makalah menyebutkan menggunakan sistem GeoDAS GIS . Mungkin menjadi petunjuk lain.