1## ----setup, include=FALSE----------------------------------------------------- 2knitr::opts_chunk$set(message = FALSE, warning = FALSE) 3 4## ----echo=TRUE,eval=TRUE,results='hide'--------------------------------------- 5library(spdep) 6if (require(rgdal, quietly=TRUE)) { 7 eire <- readOGR(system.file("shapes/eire.shp", package="spData")[1]) 8} else { 9 require(maptools, quietly=TRUE) 10 eire <- readShapeSpatial(system.file("shapes/eire.shp", package="spData")[1]) 11} 12row.names(eire) <- as.character(eire$names) 13proj4string(eire) <- CRS("+proj=utm +zone=30 +ellps=airy +units=km") 14 15## ----echo=TRUE,eval=TRUE------------------------------------------------------ 16class(eire) 17names(eire) 18 19## ----echo=TRUE,eval=TRUE------------------------------------------------------ 20fn <- system.file("etc/misc/geary_eire.txt", package="spdep")[1] 21ge <- read.table(fn, header=TRUE) 22names(ge) 23 24## ----echo=TRUE,eval=TRUE------------------------------------------------------ 25row.names(ge) <- as.character(ge$county) 26all.equal(row.names(ge), row.names(eire)) 27library(maptools) 28eire_ge <- spCbind(eire, ge) 29 30## ----echo=TRUE,eval=TRUE------------------------------------------------------ 31eire_ge1 <- eire_ge[!(row.names(eire_ge) %in% "Dublin"),] 32length(row.names(eire_ge1)) 33 34## ----echo=TRUE,eval=TRUE------------------------------------------------------ 35skewness <- function(z) {z <- scale(z, scale=FALSE); ((sum(z^3)/length(z))^2)/((sum(z^2)/length(z))^3)} 36kurtosis <- function(z) {z <- scale(z, scale=FALSE); (sum(z^4)/length(z))/((sum(z^2)/length(z))^2)} 37 38## ----echo=TRUE,eval=TRUE------------------------------------------------------ 39print(sapply(as(eire_ge1, "data.frame")[13:24], skewness), digits=3) 40print(sapply(as(eire_ge1, "data.frame")[13:24], kurtosis), digits=4) 41print(sapply(as(eire_ge1, "data.frame")[c(13,16,18,19)], function(x) skewness(log(x))), digits=3) 42print(sapply(as(eire_ge1, "data.frame")[c(13,16,18,19)], function(x) kurtosis(log(x))), digits=4) 43 44## ----echo=TRUE,eval=TRUE------------------------------------------------------ 45fn <- system.file("etc/misc/unstand_sn.txt", package="spdep")[1] 46unstand <- read.table(fn, header=TRUE) 47summary(unstand) 48 49## ----echo=TRUE,eval=TRUE------------------------------------------------------ 50class(unstand) <- c("spatial.neighbour", class(unstand)) 51of <- ordered(unstand$from) 52attr(unstand, "region.id") <- levels(of) 53unstand$from <- as.integer(of) 54unstand$to <- as.integer(ordered(unstand$to)) 55attr(unstand, "n") <- length(unique(unstand$from)) 56 57## ----echo=TRUE,eval=TRUE------------------------------------------------------ 58lw_unstand <- sn2listw(unstand) 59lw_unstand$style <- "B" 60lw_unstand 61 62## ----echo=TRUE,eval=TRUE------------------------------------------------------ 63nb <- poly2nb(eire_ge1) 64nb 65 66## ----echo=TRUE,eval=TRUE------------------------------------------------------ 67xx <- diffnb(nb, lw_unstand$neighbours, verbose=TRUE) 68 69## ----echo=TRUE,eval=FALSE,results='hide'-------------------------------------- 70# plot(eire_ge1, border="grey60") 71# plot(nb, coordinates(eire_ge1), add=TRUE, pch=".", lwd=2) 72# plot(xx, coordinates(eire_ge1), add=TRUE, pch=".", lwd=2, col=3) 73 74## ----eval=TRUE,echo=FALSE, fig.cap="County boundaries and contiguities"------- 75par(mfrow=c(1,2)) 76plot(eire_ge1, border="grey40") 77title(xlab="25 Irish counties") 78text(coordinates(eire_ge1), labels=as.character(eire_ge1$serlet), cex=0.8) 79plot(eire_ge1, border="grey60") 80title(xlab="Contiguities") 81plot(nb, coordinates(eire_ge1), add=TRUE, pch=".", lwd=2) 82plot(xx, coordinates(eire_ge1), add=TRUE, pch=".", lwd=2, col=3) 83legend("topleft", legend=c("Contiguous", "Ferry"), lwd=2, lty=1, col=c(1,3), bty="n", cex=0.7) 84par(mfrow=c(1,1)) 85 86## ----echo=FALSE,eval=TRUE----------------------------------------------------- 87load(system.file("etc/misc/raw_grass_borders.RData", package="spdep")[1]) 88 89## ----echo=TRUE,eval=FALSE,results='hide'-------------------------------------- 90# library(maptools) 91# SG <- Sobj_SpatialGrid(eire_ge1)$SG 92# library(spgrass6) 93# grass_home <- "/home/rsb/topics/grass/g64/grass-6.4.0svn" 94# initGRASS(grass_home, home=tempdir(), SG=SG, override=TRUE) 95# writeVECT6(eire_ge1, "eire", v.in.ogr_flags=c("o", "overwrite")) 96# res <- vect2neigh("eire", ID="serlet") 97 98## ----echo=TRUE,eval=TRUE------------------------------------------------------ 99grass_borders <- sn2listw(res) 100raw_borders <- grass_borders$weights 101int_tot <- attr(res, "total") - attr(res, "external") 102prop_borders <- lapply(1:length(int_tot), function(i) raw_borders[[i]]/int_tot[i]) 103dlist <- nbdists(grass_borders$neighbours, coordinates(eire_ge1)) 104inv_dlist <- lapply(dlist, function(x) 1/(x/1.609344)) 105combo_km <- lapply(1:length(inv_dlist), function(i) inv_dlist[[i]]*prop_borders[[i]]) 106combo_km_lw <- nb2listw(grass_borders$neighbours, glist=combo_km, style="B") 107summary(combo_km_lw) 108 109## ----echo=TRUE,eval=TRUE------------------------------------------------------ 110red_lw_unstand <- lw_unstand 111Clare <- which(attr(lw_unstand, "region.id") == "C") 112Kerry <- which(attr(lw_unstand, "region.id") == "H") 113Kerry_in_Clare <- which(lw_unstand$neighbours[[Clare]] == Kerry) 114Clare_in_Kerry <- which(lw_unstand$neighbours[[Kerry]] == Clare) 115red_lw_unstand$neighbours[[Clare]] <- red_lw_unstand$neighbours[[Clare]][-Kerry_in_Clare] 116red_lw_unstand$neighbours[[Kerry]] <- red_lw_unstand$neighbours[[Kerry]][-Clare_in_Kerry] 117red_lw_unstand$weights[[Clare]] <- red_lw_unstand$weights[[Clare]][-Kerry_in_Clare] 118red_lw_unstand$weights[[Kerry]] <- red_lw_unstand$weights[[Kerry]][-Clare_in_Kerry] 119summary(red_lw_unstand) 120cor(unlist(red_lw_unstand$weights), unlist(combo_km_lw$weights)) 121 122## ----echo=TRUE,eval=TRUE------------------------------------------------------ 123flatten <- function(x, digits=3, statistic="I") { 124 res <- c(format(x$estimate, digits=digits), 125 format(x$statistic, digits=digits), 126 format.pval(x$p.value, digits=digits)) 127 res <- matrix(res, ncol=length(res)) 128 colnames(res) <- paste(c("", "E", "V", "SD_", "P_"), "I", sep="") 129 rownames(res) <- deparse(substitute(x)) 130 res 131} 132`reconstructed weights` <- moran.test(eire_ge1$ocattlepacre, combo_km_lw) 133`original weights` <- moran.test(eire_ge1$ocattlepacre, red_lw_unstand) 134print(rbind(flatten(`reconstructed weights`), flatten(`original weights`)), quote=FALSE) 135 136## ----echo=TRUE,eval=TRUE------------------------------------------------------ 137eire_ge1$ln_pagval2_10 <- log(eire_ge1$pagval2_10) 138eire_ge1$ln_cowspacre <- log(eire_ge1$cowspacre) 139eire_ge1$ln_pigspacre <- log(eire_ge1$pigspacre) 140eire_ge1$ln_sheeppacre <- log(eire_ge1$sheeppacre) 141vars <- c("pagval2_10", "ln_pagval2_10", "pagval10_50", "pagval50p", 142 "cowspacre", "ln_cowspacre", "ocattlepacre", "pigspacre", 143 "ln_pigspacre", "sheeppacre", "ln_sheeppacre", "townvillp", 144 "carspcap", "radiopcap", "retailpcap", "psinglem30_34") 145nb_B <- nb2listw(lw_unstand$neighbours, style="B") 146nb_B 147lw_std <- nb2listw(lw_unstand$neighbours, glist=lw_unstand$weights, style="W") 148lw_std 149 150## ----echo=TRUE,eval=TRUE------------------------------------------------------ 151system.time({ 152MoranN <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=nb_B, randomisation=FALSE)) 153MoranR <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=nb_B, randomisation=TRUE)) 154GearyN <- lapply(vars, function(x) geary.test(eire_ge1[[x]], listw=nb_B, randomisation=FALSE)) 155GearyR <- lapply(vars, function(x) geary.test(eire_ge1[[x]], listw=nb_B, randomisation=TRUE)) 156Prop_unstdN <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=lw_unstand, randomisation=FALSE)) 157Prop_unstdR <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=lw_unstand, randomisation=TRUE)) 158Prop_stdN <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=lw_std, randomisation=FALSE)) 159Prop_stdR <- lapply(vars, function(x) moran.test(eire_ge1[[x]], listw=lw_std, randomisation=TRUE)) 160}) 161res <- sapply(c("MoranN", "MoranR", "GearyN", "GearyR", "Prop_unstdN", "Prop_unstdR", "Prop_stdN", "Prop_stdR"), function(x) sapply(get(x), "[[", "statistic")) 162rownames(res) <- vars 163ores <- res[,c(1,2,5:8)] 164 165## ----echo=FALSE,eval=TRUE------------------------------------------------------------------------- 166options("width"=100) 167 168## ----echo=TRUE,eval=TRUE-------------------------------------------------------------------------- 169print(formatC(res, format="f", digits=4), quote=FALSE) 170 171## ----echo=FALSE,eval=TRUE--------------------------------------------------------------- 172options("width"=90) 173 174## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 175wc_unstd <- spweights.constants(lw_unstand) 176wrong_N_sqVI <- sqrt((wc_unstd$nn*wc_unstd$S1 - wc_unstd$n*wc_unstd$S2 + 3*wc_unstd$S0*wc_unstd$S0)/((wc_unstd$nn-1)*wc_unstd$S0*wc_unstd$S0)-((-1/(wc_unstd$n-1))^2)) 177raw_data <- grep("^ln_", vars, invert=TRUE) 178I <- sapply(Prop_stdN, function(x) x$estimate[1])[raw_data] 179EI <- sapply(Prop_stdN, function(x) x$estimate[2])[raw_data] 180res <- (I - EI)/wrong_N_sqVI 181names(res) <- vars[raw_data] 182print(formatC(res, format="f", digits=4), quote=FALSE) 183 184## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 185res <- lapply(c("MoranR", "GearyR", "Prop_unstdR", "Prop_stdR"), function(x) sapply(get(x), function(y) c(y$estimate[1], sqrt(y$estimate[3])))) 186res <- t(do.call("rbind", res)) 187colnames(res) <- c("I", "sigma_I", "C", "sigma_C", "unstd_r", "sigma_r", "std_r", "sigma_r") 188rownames(res) <- vars 189print(formatC(res, format="f", digits=4), quote=FALSE) 190 191## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 192oMoranf <- function(x, nb) { 193 z <- scale(x, scale=FALSE) 194 n <- length(z) 195 glist <- lapply(1:n, function(i) {ii <- nb[[i]]; ifelse(ii > i, 1, 0)}) 196 lw <- nb2listw(nb, glist=glist, style="B") 197 wz <- lag(lw, z) 198 I <- (sum(z*wz)/sum(z*z)) 199 I 200} 201res <- sapply(vars, function(x) oMoranf(eire_ge1[[x]], nb=lw_unstand$neighbours)) 202print(formatC(res, format="f", digits=4), quote=FALSE) 203 204## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 205MoranI.boot <- function(var, i, ...) { 206 var <- var[i] 207 return(moran(x=var, ...)$I) 208} 209Nsim <- function(d, mle) { 210 n <- length(d) 211 rnorm(n, mle$mean, mle$sd) 212} 213f_bperm <- function(x, nsim, listw) { 214 boot(x, statistic=MoranI.boot, R=nsim, sim="permutation", listw=listw, 215 n=length(x), S0=Szero(listw)) 216} 217f_bpara <- function(x, nsim, listw) { 218 boot(x, statistic=MoranI.boot, R=nsim, sim="parametric", ran.gen=Nsim, 219 mle=list(mean=mean(x), sd=sd(x)), listw=listw, n=length(x), 220 S0=Szero(listw)) 221} 222nsim <- 4999 223set.seed(1234) 224 225## ----echo=TRUE,eval=FALSE--------------------------------------------------------------- 226# system.time({ 227# MoranNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=nb_B)) 228# MoranRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=nb_B)) 229# Prop_unstdNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=lw_unstand)) 230# Prop_unstdRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=lw_unstand)) 231# Prop_stdNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=lw_std)) 232# Prop_stdRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=lw_std)) 233# }) 234 235## ----echo=FALSE,eval=FALSE-------------------------------------------------------------- 236# zzz <- system.time({ 237# MoranNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=nb_B)) 238# MoranRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=nb_B)) 239# Prop_unstdNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=lw_unstand)) 240# Prop_unstdRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=lw_unstand)) 241# Prop_stdNb <- lapply(vars, function(x) f_bpara(x=eire_ge1[[x]], nsim=nsim, listw=lw_std)) 242# Prop_stdRb <- lapply(vars, function(x) f_bperm(x=eire_ge1[[x]], nsim=nsim, listw=lw_std)) 243# }) 244# res <- lapply(c("MoranNb", "MoranRb", "Prop_unstdNb", "Prop_unstdRb", "Prop_stdNb", "Prop_stdRb"), function(x) sapply(get(x), function(y) (y$t0 - mean(y$t))/sd(y$t))) 245# res <- t(do.call("rbind", res)) 246# colnames(res) <- c("MoranNb", "MoranRb", "Prop_unstdNb", "Prop_unstdRb", "Prop_stdNb", "Prop_stdRb") 247# rownames(res) <- vars 248# save(zzz, res, file="backstore/boot_res.RData") 249 250## ----echo=FALSE,eval=FALSE-------------------------------------------------------------- 251# bsfn <- system.file("etc/backstore/boot_res.RData", package="spdep") 252# load(bsfn) 253# zzz 254 255## ----echo=TRUE,eval=FALSE--------------------------------------------------------------- 256# res <- lapply(c("MoranNb", "MoranRb", "Prop_unstdNb", "Prop_unstdRb", "Prop_stdNb", "Prop_stdRb"), function(x) sapply(get(x), function(y) (y$t0 - mean(y$t))/sd(y$t))) 257# res <- t(do.call("rbind", res)) 258# colnames(res) <- c("MoranNb", "MoranRb", "Prop_unstdNb", "Prop_unstdRb", "Prop_stdNb", "Prop_stdRb") 259# rownames(res) <- vars 260 261## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 262print(formatC(res, format="f", digits=4), quote=FALSE) 263oores <- ores - res 264apply(oores, 2, mad) 265alpha_0.05 <- qnorm(0.05, lower.tail=FALSE) 266all((res >= alpha_0.05) == (ores >= alpha_0.05)) 267 268## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 269lm_objs <- lapply(vars, function(x) lm(formula(paste(x, "~1")), data=eire_ge1)) 270system.time({ 271MoranSad <- lapply(lm_objs, function(x) lm.morantest.sad(x, listw=nb_B)) 272Prop_unstdSad <- lapply(lm_objs, function(x) lm.morantest.sad(x, listw=lw_unstand)) 273Prop_stdSad <- lapply(lm_objs, function(x) lm.morantest.sad(x, listw=lw_std)) 274}) 275res <- sapply(c("MoranSad", "Prop_unstdSad", "Prop_stdSad"), function(x) sapply(get(x), "[[", "statistic")) 276rownames(res) <- vars 277 278## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 279print(formatC(res, format="f", digits=4), quote=FALSE) 280oores <- res - ores[,c(1,3,5)] 281apply(oores, 2, mad) 282all((res >= alpha_0.05) == (ores[,c(1,3,5)] >= alpha_0.05)) 283 284## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 285system.time({ 286MoranEx <- lapply(lm_objs, function(x) lm.morantest.exact(x, listw=nb_B)) 287Prop_unstdEx <- lapply(lm_objs, function(x) lm.morantest.exact(x, listw=lw_unstand)) 288Prop_stdEx <- lapply(lm_objs, function(x) lm.morantest.exact(x, listw=lw_std)) 289}) 290res <- sapply(c("MoranEx", "Prop_unstdEx", "Prop_stdEx"), function(x) sapply(get(x), "[[", "statistic")) 291rownames(res) <- vars 292 293## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 294print(formatC(res, format="f", digits=4), quote=FALSE) 295oores <- res - ores[,c(1,3,5)] 296apply(oores, 2, mad) 297all((res >= alpha_0.05) == (ores[,c(1,3,5)] >= alpha_0.05)) 298 299## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 300vars_scaled <- lapply(vars, function(x) scale(eire_ge1[[x]], scale=FALSE)) 301nb_W <- nb2listw(lw_unstand$neighbours, style="W") 302pre <- spdep:::preAple(0, listw=nb_W) 303MoranAPLE <- sapply(vars_scaled, function(x) spdep:::inAple(x, pre)) 304pre <- spdep:::preAple(0, listw=lw_std, override_similarity_check=TRUE) 305Prop_stdAPLE <- sapply(vars_scaled, function(x) spdep:::inAple(x, pre)) 306res <- cbind(MoranAPLE, Prop_stdAPLE) 307colnames(res) <- c("APLE W", "APLE Gstd") 308rownames(res) <- vars 309 310## ----echo=TRUE,eval=TRUE---------------------------------------------------------------- 311print(formatC(res, format="f", digits=4), quote=FALSE) 312 313## ----results='asis',eval=TRUE,echo=FALSE, fig.cap="Three contrasted spatial weights definitions"---- 314pal <- grey.colors(9, 1, 0.5, 2.2) 315oopar <- par(mfrow=c(1,3), mar=c(1,1,3,1)+0.1) 316z <- t(listw2mat(nb_B)) 317brks <- c(0,0.1,1) 318image(1:25, 1:25, z[,ncol(z):1], breaks=brks, col=pal[c(1,9)], 319 main="Binary", axes=FALSE) 320box() 321z <- t(listw2mat(lw_unstand)) 322brks <- c(0,quantile(c(z)[c(z) > 0], seq(0,1,1/8))) 323image(1:25, 1:25, z[,ncol(z):1], breaks=brks, col=pal, main="General", axes=FALSE) 324box() 325z <- t(listw2mat(lw_std)) 326brks <- c(0,quantile(c(z)[c(z) > 0], seq(0,1,1/8))) 327image(1:25, 1:25, z[,ncol(z):1], breaks=brks, col=pal, 328 main="Std. general", axes=FALSE) 329box() 330par(oopar) 331 332## ----results='asis',eval=TRUE,echo=FALSE------------------------------------------------ 333eire_ge1$nb_B <- sapply(nb_B$weights, sum) 334eire_ge1$lw_unstand <- sapply(lw_unstand$weights, sum) 335library(lattice) 336trellis.par.set(sp.theme()) 337p1 <- spplot(eire_ge1, c("nb_B"), main="Binary") 338p2 <- spplot(eire_ge1, c("lw_unstand"), main="General") 339print(p1, split=c(1,1,2,1), more=TRUE) 340print(p2, split=c(2,1,2,1), more=FALSE) 341 342