Commit 243ed26b authored by Andreas Tille's avatar Andreas Tille

New upstream version 0.7-7

parent 82b13a5d
Package: spdep
Version: 0.7-4
Date: 2017-11-12
Version: 0.7-7
Date: 2018-04-03
Title: Spatial Dependence: Weighting Schemes, Statistics and Models
Encoding: UTF-8
Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bivand@nhh.no", comment=c(ORCID="0000-0003-2392-6140")), person("Micah", "Altman", role = "ctb"),
......@@ -31,6 +31,7 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bi
person("Gianfranco", "Piras", role = "ctb"),
person("Markus", "Reder", role = "ctb"),
person("Michael", "Tiefelsdorf", role = "ctb"),
person("René", "Westerholt", role="ctb"),
person("Danlin", "Yu", role = "ctb"))
Depends: R (>= 3.0.0), methods, sp (>= 1.0), Matrix (>= 1.0.12), spData
(>= 0.2.6.0)
......@@ -49,8 +50,9 @@ Description: A collection of functions to create spatial weights matrix
including global 'Morans I', 'APLE', 'Gearys C', 'Hubert/Mantel' general cross
product statistic, Empirical Bayes estimates and 'Assunção/Reis' Index,
'Getis/Ord' G and multicoloured join count statistics, local 'Moran's I'
and 'Getis/Ord' G, 'saddlepoint' approximations and exact tests for global
and local 'Moran's I'; and functions for estimating spatial simultaneous
and 'Getis/Ord' G, 'saddlepoint' approximations, exact tests for global
and local 'Moran's I' and 'LOSH' local indicators of spatial heteroscedasticity;
and functions for estimating spatial simultaneous
'autoregressive' ('SAR') lag and error models, impact measures for lag
models, weighted and 'unweighted' 'SAR' and 'CAR' spatial regression models,
semi-parametric and Moran 'eigenvector' spatial filtering, 'GM SAR' error
......@@ -58,8 +60,8 @@ Description: A collection of functions to create spatial weights matrix
License: GPL (>= 2)
VignetteBuilder: knitr
NeedsCompilation: yes
Packaged: 2017-11-20 08:33:45 UTC; rsb
Author: Roger Bivand [cre, aut] (0000-0003-2392-6140),
Packaged: 2018-04-03 08:14:15 UTC; rsb
Author: Roger Bivand [cre, aut] (<https://orcid.org/0000-0003-2392-6140>),
Micah Altman [ctb],
Luc Anselin [ctb],
Renato Assunção [ctb],
......@@ -88,7 +90,8 @@ Author: Roger Bivand [cre, aut] (0000-0003-2392-6140),
Gianfranco Piras [ctb],
Markus Reder [ctb],
Michael Tiefelsdorf [ctb],
René Westerholt [ctb],
Danlin Yu [ctb]
Maintainer: Roger Bivand <Roger.Bivand@nhh.no>
Repository: CRAN
Date/Publication: 2017-11-22 10:04:21 UTC
Date/Publication: 2018-04-03 21:37:35 UTC
2428ff9c8fc1c5390035b6f09bcc9b13 *ChangeLog
4139082680c932812dc6a0ef7fc5b5c8 *DESCRIPTION
d1d829561706282fb0752406c39081c4 *NAMESPACE
16a78bd3b242db52a6fef6ac6d0687c4 *DESCRIPTION
66845965ba6a0335e7f5047a2e6335b0 *NAMESPACE
ca3af94aff175d2e462ffb33bd9ac568 *R/AAA.R
b0e1b90e65dd7b62bdd3948361cb5237 *R/EBI.R
8010b6a2e0a46f013ad55ec541812c8e *R/LOSH.R
94be548ba7c37f33ef8076b74533092f *R/LOSH.cs.R
d605828e91b5a8f9e61b7788949cecc5 *R/LOSH.mc.R
0fbd6d0a73bd6f6b204b4a5fadc045c9 *R/LR.spsarlm.R
779dcd05ca5740e251c19f7d0aae899a *R/ME.R
61a73d8c7ae411b1fd3bbc08074ce580 *R/SpatialFiltering.R
......@@ -24,15 +27,15 @@ d6b8860ef2337d9746ff5d1119f36590 *R/droplinks.R
1c7c266113c4f079d2f0c6d9ccc5141f *R/edit.nb.R
f95f17856064f449cf1603d083f8ec87 *R/eigenw.R
6b22c57e33fb394b821be58faf1bb1d2 *R/error.hessian.R
8cd1a56343aef6870ff7d24c4dd8f2fd *R/error.spsarlm.R
3b4d8cea68dde978d297833e7247f82c *R/error.spsarlm.R
cbe45866b710921490b6eea92e0fa7a2 *R/gabrielneigh.R
5c7fc104be2f03783a4193fbf2351b19 *R/geary.R
c419021fe036e04b29c3cc9c38cda874 *R/globalG.R
494162cf871ebd36494283ccde5dc1df *R/globalG.R
d6320c56e705e0f9e81e46a851a74fa0 *R/graph2nb.R
7d90621f339f34ea56997a6f9d30d661 *R/impacts.R
d37880c837d66c5167434da6858c1c40 *R/jacobian.R
f842c1af5dcab63b5f1eb0bd3099341a *R/jacobian_setup.R
4eac610dd1a2b7e624f3cc47696e1d91 *R/jc.R
e9042b5d26783b2224f601afca856bbf *R/jc.R
a22b5a5326ae3f1ab5fc2122212ad5b0 *R/knearneigh.R
c71eef240bf2e5614276b28ab0f01af7 *R/knn2nb.R
6226042eae5ffc30b37891544a60507a *R/kpgm_new.R
......@@ -48,12 +51,12 @@ d41377ef1fcf3c0d8001eb6f853f2337 *R/listw2Matrix.R
4a9458975190059999b49e22ba3f4c05 *R/listw2sn.R
3f616a97e845028a86cecded418390a1 *R/lm.LMtests.R
6642336a28d6b6a1dc3dad6f81923f19 *R/lm.morantest.R
2d9f29e96c8851ab8619b5c4d4d85771 *R/localG.R
6be3f46cfa21d316513b9dc802f10990 *R/localmoran.R
e01b338dd1d95a4d4576af188c98ef8b *R/localG.R
788edea8fed2b1f37238cf57c384d094 *R/localmoran.R
955f0e2897bec106a72a13cf0d4cb83b *R/localmoran.exact.R
03998c2eafe6a603e5301f2377bf06d7 *R/mcmcsamp.R
ac6f3cf2d5c1f7145ae54ba0b1135176 *R/mess.R
f1632cb1c51eea17fa2c5993cdf953f0 *R/moran.R
6b136a9a9cb241c672b2e156de899233 *R/moran.R
f583b597ad81dd8f1468a73779acce11 *R/moran.exact.R
acc3f60a11d6251a64157e3b4802b3f6 *R/moran.exact.alt.R
13c1ddfb0473badd8c20d92c7c5aa3e0 *R/moran.plot.R
......@@ -76,7 +79,7 @@ a73ab332d82c29834dce2b5b650822b1 *R/plot.mst.R
59a5f3980e8585d96ac511d4d23d669f *R/prunecost.R
58fc4f67d40bb82e24a32af9ee5ee0c2 *R/prunemst.R
72519954bdafc8521aa950d7d3b67855 *R/read.gal.R
ef6f3184db69b4f7cfe934baf99b57c6 *R/read.gwt2nb.R
77b2d8c576414616e13c5ff987388b15 *R/read.gwt2nb.R
64633a828f14988e6d3aa97d6ce8ea2d *R/relneigh.R
40534663fb14a40fc659768038487304 *R/rotation.R
b8ae91f8e8af0ad2ddb651840f879a3d *R/s2sls.R
......@@ -96,7 +99,7 @@ afc9a6a96713d8c63d8459cf7a6e4ad5 *R/summary.spsarlm.R
79530fa1378ad438ea0695b44a3b134f *R/tri2nb.R
e4c7a7d4e336723499fe98a9e13fd67f *R/utils.R
e0483061c02957fd53635f7762ce617b *R/weights-utils.R
b46a4e2c35b6dfc039a9fb74eafe7110 *build/vignette.rds
23cf9e550c30850e6b988b4e2981ab3f *build/vignette.rds
b818650dc1a3e0335a7767fe78f721f5 *data/columbus.rda
5e7e97b208f3d4a44393e668799ce382 *data/eire.rda
3c410289cf12787247a95ff858bf4a97 *data/oldcol.rda
......@@ -105,24 +108,24 @@ eaa320807543a7d0b75ca403b68a5397 *inst/CITATION
2428ff9c8fc1c5390035b6f09bcc9b13 *inst/ChangeLog
99c4d523114b0e2769464b65a9cf1150 *inst/doc/CO69.R
709cb9502be2c6176e03747da963bc15 *inst/doc/CO69.Rnw
8e7e52e058773cb6553b68ff569a13c0 *inst/doc/CO69.pdf
1456076e416657a0edfdab25123d0e57 *inst/doc/CO69.pdf
27fb91633d64bff1470c754abdc09ac1 *inst/doc/SpatialFiltering.R
00a3b786362437b0d02c578910457a5e *inst/doc/SpatialFiltering.Rnw
0f1034c5ddb9f21ddacc1ae554efdc22 *inst/doc/SpatialFiltering.pdf
584ef55153e6b402eb660e80b38ffd35 *inst/doc/SpatialFiltering.pdf
68a52b4a10480b4a2adedb18494c7deb *inst/doc/backstore/boot_res.RData
3889914095150cc9ed13c578b5f51a73 *inst/doc/backstore/nyME_res.RData
e12b7288630644ee28894e5b04cea06c *inst/doc/nb.R
96bd583488e35a9e980cbe5b79f57d46 *inst/doc/nb.Rnw
e0e0f9baf855e52d15d20fba7901a3b3 *inst/doc/nb.pdf
143b31deb7697a9cbac1ad15b5397085 *inst/doc/nb.pdf
52619007cbaa4f7cb7e64423347d720c *inst/doc/nb_igraph.R
0a4d2c611361f8e9242b91d09a7d7354 *inst/doc/nb_igraph.Rmd
3e8f0ec5cdb61bddb75ba4deb710eabe *inst/doc/nb_igraph.html
a4aaaca85143eeb52b01b0d3c826f8b0 *inst/doc/nb_igraph.html
ee80ddbeeca2900bbece5e203a657817 *inst/doc/nb_sf.R
da968c529f783ba884f37f27ac6b90fd *inst/doc/nb_sf.Rmd
11951d6c2ea963b802051ac536d4317f *inst/doc/nb_sf.html
9d599c66add8a55cfeae25a96e949cb6 *inst/doc/nb_sf.html
98ba2b15ccd8b39b51dfd25ed6a25d15 *inst/doc/sids.R
af6bf80e9a912760fa14ff55f82d2d00 *inst/doc/sids.Rnw
cedaeb77fc27a5e4722cfc1e6cff343d *inst/doc/sids.pdf
fbc53850e2e408387f838a64cc059223 *inst/doc/sids.pdf
53178ac43dd150bbb835b63e018f0cd4 *inst/etc/misc/geary_eire.txt
e3fb11232e6cb27770d15aaf4f632aa4 *inst/etc/misc/nc_xtra.dbf
2591ac95bd9191d0aee07450f7caf9bb *inst/etc/misc/nyadjwts.dbf
......@@ -148,6 +151,9 @@ d903c1e405c15c03ab1a04d2ad56b0fd *man/EBImoran.mc.Rd
2e798ab52cf6472d8e176adef1fc1e9a *man/EBest.Rd
a28e4db93466660e68b436bc624d8bb3 *man/EBlocal.Rd
aac6a73fab3a3e0ffdc4847561b589ee *man/GMerrorsar.Rd
bfc736ca4756763541c5eb7e47188842 *man/LOSH.Rd
213bee2cea27ac466d2e2c06d6610be8 *man/LOSH.cs.Rd
8fddf844d045ccde4bc4d29786329e82 *man/LOSH.mc.Rd
baf7aebaead23be5de6cc920ee49a638 *man/LR.sarlm.Rd
f9e8b18b2939f12f8bc39a5d5489a43b *man/MCMCsamp.Rd
7b8d9a985f1c3f434e82ea8c80f267c1 *man/ME.Rd
......@@ -183,10 +189,10 @@ e876fa32337f7b7e9dbe701e84dc33c3 *man/graphneigh.Rd
fe5102da5d3290cdbb47395805511ec2 *man/gstsls.Rd
6a5e7a5a456e45fa51ac5121911b7b87 *man/impacts.sarlm.Rd
3ebb36e6b2d2b8aea7d122dcfdddb527 *man/include.self.Rd
91829bb2d17e32b1ceafa52906c0fedc *man/invIrM.Rd
9a65838851fd6b5e78cc853b3eab8dd1 *man/invIrM.Rd
3a5aa076bb3c12d7c608a5d57eae9fb7 *man/joincount.mc.Rd
3ef157b0a852acb019081ea67d9b3c0c *man/joincount.multi.Rd
e68b6065a95c7fbb05ccebcc9a11f449 *man/joincount.test.Rd
3ca8dd3ff7a52d1c2fa8ce590e66fe70 *man/joincount.test.Rd
ff43f248f85add1589a1f8c97f3b459d *man/knearneigh.Rd
3666112819e0fb6937f4684147755588 *man/knn2nb.Rd
530a137b1bad26cbbc658a52e2f48562 *man/lag.listw.Rd
......@@ -202,14 +208,14 @@ bf22c7dc637db3a45d7a63839a7269a8 *man/lm.morantest.Rd
752bd7f99bcb86600a677a854e7de9b5 *man/lm.morantest.exact.Rd
ae4f18b438f4ecc640d75a696628e021 *man/lm.morantest.sad.Rd
6270f17997be10a78ef6fce75d4b31e5 *man/localG.Rd
6fcbc451d3e3bfe4c608509d078df98f *man/localmoran.Rd
306d794d53252687ea9b4cf3f02289ba *man/localmoran.Rd
70ed836efe7eef4514ef1749cd9b4680 *man/localmoran.exact.Rd
f508c2ce9a3c4cf78e5d60894927fb14 *man/localmoran.sad.Rd
ea6173d1d51344b6ff67f3e6d8652d40 *man/mat2listw.Rd
9bb3b49ce5a324c53c2fd9e6c408f4bd *man/moran.Rd
463bf13f2d8711e0b052c31afba57473 *man/moran.mc.Rd
6170704c9f5410a63268d8ece463d0a9 *man/moran.plot.Rd
a3fec8c2c6b261954cc743a53974fec9 *man/moran.test.Rd
5e9fb24a502ee01751cd2fd14f897773 *man/moran.test.Rd
9cadacbbbcebab440f54cc4d6ced36d8 *man/mstree.Rd
efe41773577381cd5b9db056701674e7 *man/nb2INLA.Rd
ee1cfe442c1f97b5a1d6d19079721d11 *man/nb2WB.Rd
......@@ -231,7 +237,7 @@ a422bc9ac4822cc263438f83fba90d5e *man/predict.sarlm.Rd
09b2776bb47f80cdbe5562d88b36ef18 *man/prunecost.Rd
f592261ce51e820587da92afaa7cd853 *man/prunemst.Rd
f8dc0b0518d2badebbd7f93e552a3f73 *man/read.gal.Rd
fa145967d30f1fcc2e680a045fe52852 *man/read.gwt2nb.Rd
fe86bf830f291a16d60102c52bfdf01e *man/read.gwt2nb.Rd
22587da32f258cb04ac771791ff2a5f4 *man/residuals.sarlm.Rd
e3f28cd38c69d847efdf6143af31637d *man/rotation.Rd
7395261393ac6e81d3523ebe35582603 *man/sacsarlm.Rd
......@@ -267,7 +273,7 @@ ab3174ddddfeeb9f99bb8395975d9e31 *src/eminmaxC.c
18eb40c98c79a76c1857e5d5bf904c55 *src/insiders.c
a34376d431d495f08af5623b7132b71c *src/jc.c
219854edc76315acc25b0e01e74b8e94 *src/knn.c
ef5c9e567ee7400345ab2915060aa99e *src/lagw.c
2f356573b774ef51511fcfba1ea00aaf *src/lagw.c
33aadfbe6cbd73fd38ee803d9b1cc670 *src/listw2Matrix.c
9d048d14d51b07c40be4ccd015366869 *src/listw2sn.c
d994bae2c2a98827cb490b8ec2748cae *src/ml_sse.c
......
......@@ -114,6 +114,8 @@ export(lee.mc, lee, lee.test)
export(spBreg_lag)
export(LOSH, LOSH.cs, LOSH.mc)
#import(maptools)
S3method(print, nb)
......
LOSH <- function(x, listw, a = 2, var_hi = TRUE, zero.policy = NULL, na.action = na.fail, spChk = NULL) {
if (is.null(zero.policy))
zero.policy <- get("zeroPolicy", envir = .spdepOptions)
stopifnot(is.logical(zero.policy))
n <- length(listw$neighbours)
#a <- 2 ## "a" could be any other positive value, but the chi-square-based inference is then no longer possible
if (n != length(x))
stop("Different numbers of observations")
NAOK <- deparse(substitute(na.action)) == "na.pass"
if(var_hi) {
res <- matrix(nrow = n, ncol = 6)
colnames(res) <- c("Hi", "E.Hi", "Var.Hi", "Z.Hi", "x_bar_i", "ei")
} else {
res <- matrix(nrow = n, ncol = 3)
colnames(res) <- c("Hi", "x_bar_i", "ei")
}
### calculation of the row sums of the spatial weights
Wi <- vapply(listw$weights, sum, FUN.VALUE = 0.0)
### calculation of x_bar_i
res[,"x_bar_i"] <- lag.listw(listw, x, zero.policy = zero.policy, NAOK = NAOK) / Wi
### calculation of ej
res[,"ei"] <- abs(x - res[,"x_bar_i"])^a
### calculation of 1/(h1 * Wi)
denom_hi <- mean(res[,"ei"]) * Wi
### calculation of Hi
res[,"Hi"] <- lag.listw(listw, res[,"ei"], zero.policy = zero.policy, NAOK = NAOK) / denom_hi
if(var_hi) {
### calculation of the variance of the ei terms (a global value)
var_ei <- (sum(res[,"ei"]^2) / n) - mean(res[,"ei"])^2
### calculation of the variance of Hi
res[,"Var.Hi"] <- (n-1)^(-1) * denom_hi^(-2) * var_ei * (n * vapply(listw$weights, function(y) sum(y^2), FUN.VALUE = 0.0) - Wi^2)
### calculation of Zi
res[,"Z.Hi"] <- (2 * res[,"Hi"]) / res[,"Var.Hi"]
### storing the expectations
res[,"E.Hi"] <- 1
}
res
}
LOSH.cs <- function(x, listw, zero.policy = NULL, na.action = na.fail,
p.adjust.method = "none", spChk = NULL) {
stopifnot(is.vector(x))
if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
if (is.null(zero.policy))
zero.policy <- get("zeroPolicy", envir = .spdepOptions)
stopifnot(is.logical(zero.policy))
if (!is.null(attr(listw$neighbours, "self.included")) &&
attr(listw$neighbours, "self.included"))
stop("Self included among neighbours")
if (is.null(spChk))
spChk <- get.spChkOption()
if (spChk && !chkIDs(x, listw))
stop("Check of data and weights ID integrity failed")
if (!is.numeric(x))
stop(paste(deparse(substitute(x)), "is not a numeric vector"))
# NAOK <- deparse(substitute(na.action)) == "na.pass"
x <- na.action(x)
na.act <- attr(x, "na.action")
rn <- attr(listw, "region.id")
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
excl <- class(na.act) == "exclude"
}
### calculate LOSH for all locations
res <- LOSH(x, listw, 2, TRUE, zero.policy, na.action, spChk) #### a = 2 is used, because the chi-square-based inference is only valid for a measure of variance.
### add a column for the chi-square-based p-values
res <- cbind(res, "Pr()" = 1)
### p-values using the non-central chi-square distribution
res[,"Pr()"] <- pchisq(res[,"Z.Hi"], df = 2/res[,"Var.Hi"], lower.tail = FALSE)
### correction for multiple hypothesis testing (if p.adjust.method != "none")
res[,"Pr()"] <- p.adjustSP(res[,"Pr()"], listw$neighbours, method = p.adjust.method)
if (!is.null(na.act) && excl) {
res <- naresid(na.act, res)
}
if (!is.null(rn))
rownames(res) <- rn
attr(res, "call") <- match.call()
if (!is.null(na.act))
attr(res, "na.action") <- na.act
class(res) <- c("LOSH", class(res))
res
}
LOSH.mc <- function(x, listw, a = 2, nsim = 99, zero.policy = NULL, na.action = na.fail,
spChk = NULL, adjust.n = TRUE, p.adjust.method = "none") {
stopifnot(is.vector(x))
if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
if (is.null(zero.policy))
zero.policy <- get("zeroPolicy", envir = .spdepOptions)
stopifnot(is.logical(zero.policy))
if (!is.null(attr(listw$neighbours, "self.included")) &&
attr(listw$neighbours, "self.included"))
stop("Self included among neighbours")
if (is.null(spChk))
spChk <- get.spChkOption()
if (spChk && !chkIDs(x, listw))
stop("Check of data and weights ID integrity failed")
if (!is.numeric(x))
stop(paste(deparse(substitute(x)), "is not a numeric vector"))
if (missing(nsim))
stop("nsim must be given")
cards <- card(listw$neighbours)
if (!zero.policy && any(cards == 0))
stop("regions with no neighbours found")
if (deparse(substitute(na.action)) == "na.pass")
stop("na.pass not permitted")
x <- na.action(x)
na.act <- attr(x, "na.action")
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
n <- length(listw$neighbours)
rn <- attr(listw, "region.id")
if (n != length(x))
stop("objects of different length")
gamres <- suppressWarnings(nsim > gamma(n + 1))
if (gamres)
stop("nsim too large for this number of observations")
if (nsim < 1)
stop("nsim too small")
if (adjust.n)
n <- n - sum(cards == 0L)
res <- LOSH(x, listw, a, FALSE, zero.policy, na.action, spChk)
res <- cbind(res, "Pr()" = 1)
losh_boot <- function(data, indices, curr_i, ...) {
var <- data[indices]
var[curr_i] <- data[curr_i]
return(LOSH(x = var, ...)[curr_i,"Hi"])
}
cores <- get.coresOption()
if (is.null(cores)) {
parallel <- "no"
}
else {
parallel <- ifelse(get.mcOption(), "multicore", "snow")
}
ncpus <- ifelse(is.null(cores), 1L, cores)
cl <- NULL
if (parallel == "snow") {
cl <- parallel::makeCluster(get.coresOption())
parallel::clusterExport(cl, list("LOSH", "lag.listw"), envir = environment())
on.exit(parallel::stopCluster(cl))
if (is.null(cl)) {
parallel <- "no"
warning("no cluster in ClusterOption, parallel set to no")
}
}
pvals <- numeric(length = length(x))
for(curr_i in 1:length(x)) {
boot_obj <- boot(x, statistic = losh_boot, curr_i = curr_i, R = nsim, sim = "permutation",
listw = listw, a = a, var_hi = FALSE, zero.policy = zero.policy,
na.action = na.action, spChk = spChk, parallel = parallel, ncpus = ncpus, cl = cl)
resampled <- append(boot_obj$t, res[curr_i, "Hi"])
rankboot <- rank(unlist(resampled))
xrank <- rankboot[length(resampled)]
diff <- nsim - xrank
diff <- ifelse(diff > 0, diff, 0)
pval <- punif((diff + 1)/(nsim + 1))
if (!is.finite(pval) || pval < 0 || pval > 1)
warning("Out-of-range p-value: reconsider test arguments")
pvals[curr_i] <- pval
}
res[,"Pr()"] <- pvals
res[,"Pr()"] <- p.adjustSP(res[,"Pr()"], listw$neighbours, method = p.adjust.method)
if (!is.null(rn))
rownames(res) <- rn
if (!is.null(na.act))
attr(res, "na.action") <- na.act
class(res) <- c("htest", "mc.sim", "LOSH", class(res))
res
}
......@@ -459,11 +459,16 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, zero.p
WX <- create_WX(x, listw, zero.policy=zero.policy, prefix="lag")
x <- cbind(x, WX)
lm.model <- lm(y ~ x - 1, weights=weights)
# 180128 Mark L. Burkey summary.lm error for SLX object
colnames(x) <- make.names(colnames(x))
if (attr(mt, "intercept") == 1L) {
lm.model <- lm(formula(paste("y ~ ", paste(colnames(x)[-1], collapse="+"))), data=as.data.frame(x), weights=weights)
} else {
lm.model <- lm(formula(paste("y ~ 0 + ", paste(colnames(x), collapse="+"))), data=as.data.frame(x), weights=weights)
}
sum_lm_model <- summary.lm(lm.model, correlation = FALSE)
mixedImps <- NULL
K <- ifelse(isTRUE(grep("\\(Intercept\\)",
K <- ifelse(isTRUE(grep("Intercept",
names(coefficients(lm.model))[1]) == 1L), 2, 1)
m <- length(coefficients(lm.model))
odd <- (m%/%2) > 0
......
# Copyright 2002-2017 by Hisaji ONO and Roger Bivand
# Copyright 2002-2018 by Hisaji ONO and Roger Bivand
#
# General G Statistics
#
......
......@@ -15,12 +15,13 @@ joincount <- function(dums, listw) {
}
joincount.test <- function(fx, listw, zero.policy=NULL,
alternative="greater", #adjust.n=TRUE,
alternative="greater", sampling="nonfree",
spChk=NULL, adjust.n=TRUE) {
if (is.null(zero.policy))
zero.policy <- get("zeroPolicy", envir = .spdepOptions)
stopifnot(is.logical(zero.policy))
alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
sampling <- match.arg(sampling, c("nonfree", "free"))
if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
"is not a listw object"))
if (!is.factor(fx)) stop(paste(deparse(substitute(x)),
......@@ -57,13 +58,21 @@ joincount.test <- function(fx, listw, zero.policy=NULL,
# wc$n2 <- N-2
# wc$n3 <- N-3
# }
Ejc <- (wc$S0*(ntab*(ntab-1))) / (2*N*wc$n1)
if (sampling == "nonfree") {
Ejc <- (wc$S0*(ntab*(ntab-1))) / (2*N*wc$n1) # CO 1981 p. 20 (1.31)
Vjc <- (wc$S1*(ntab*(ntab-1))) / (N*wc$n1)
Vjc <- Vjc + (((wc$S2 - 2*wc$S1)*ntab*(ntab-1)*(ntab-2)) /
(N*wc$n1*wc$n2))
Vjc <- Vjc + (((S02 + wc$S1 - wc$S2)*ntab*(ntab-1)*(ntab-2)*
(ntab-3)) / (N*wc$n1*wc$n2*wc$n3))
Vjc <- (0.25 * Vjc) - Ejc^2
Vjc <- (0.25 * Vjc) - Ejc^2 # CO 1981 p. 20 (1.32)
} else if (sampling == "free") {
p <- ntab/n
Ejc <- (wc$S0*(p^2))/2 # CO 1981 p. 20 (1.25)
Vjc <- ((wc$S1*(p^2)) +
((wc$S2-2*wc$S1)*(p^3)) +
((wc$S1-wc$S2)*(p^4)))/4 # CO 1981 p. 20 (1.26)
} else stop("sampling must be nonfree or free")
for (i in 1:nBB) {
estimate <- c(BB5[i], Ejc[i], Vjc[i])
names(estimate) <- c("Same colour statistic",
......@@ -80,7 +89,8 @@ joincount.test <- function(fx, listw, zero.policy=NULL,
if (!is.finite(p.value) || p.value < 0 || p.value > 1)
warning("Out-of-range p-value: reconsider test arguments")
}
method <- "Join count test under nonfree sampling"
method <- paste0("Join count test under ", sampling,
" sampling")
data.name <- paste(deparse(substitute(fx)), "\nweights:",
deparse(substitute(listw)), "\n")
res[[i]] <- list(statistic=statistic, p.value=p.value,
......
# Copyright 2001-3 by Roger Bivand
# Copyright 2001-18 by Roger Bivand
#
localG <- function(x, listw, zero.policy=NULL, spChk=NULL, return_internals=FALSE, GeoDa=FALSE) {
......
# Copyright 2001-12 by Roger Bivand
# Copyright 2001-18 by Roger Bivand
#
localmoran <- function(x, listw, zero.policy=NULL, na.action=na.fail,
alternative = "greater", p.adjust.method="none", mlvar=TRUE,
spChk=NULL) {
spChk=NULL, adjust.x=FALSE) {
stopifnot(is.vector(x))
if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
......@@ -34,19 +34,45 @@ localmoran <- function(x, listw, zero.policy=NULL, na.action=na.fail,
else if (alternative == "greater") Prname <- "Pr(z > 0)"
else Prname <- "Pr(z < 0)"
colnames(res) <- c("Ii", "E.Ii", "Var.Ii", "Z.Ii", Prname)
if (adjust.x) {
nc <- card(listw$neighbours) > 0L
xx <- mean(x[nc], na.rm=NAOK)
} else {
xx <- mean(x, na.rm=NAOK)
}
z <- x - xx
lz <- lag.listw(listw, z, zero.policy=zero.policy, NAOK=NAOK)
if (mlvar) s2 <- sum(z^2, na.rm=NAOK)/n
else s2 <- sum(z^2, na.rm=NAOK)/(n-1)
if (mlvar) {
if (adjust.x) {
s2 <- sum(z[nc]^2, na.rm=NAOK)/sum(nc)
} else {
s2 <- sum(z^2, na.rm=NAOK)/n
}
} else {
if (adjust.x) {
s2 <- sum(z[nc]^2, na.rm=NAOK)/(sum(nc)-1)
} else {
s2 <- sum(z^2, na.rm=NAOK)/(n-1)
}
}
res[,1] <- (z/s2) * lz
Wi <- sapply(listw$weights, sum)
res[,2] <- -Wi / (n-1)
if (mlvar) b2 <- (sum(z^4, na.rm=NAOK)/n)/(s2^2)
else b2 <- (sum(z^4, na.rm=NAOK)/(n-1))/(s2^2)
if (mlvar) {
if (adjust.x) {
b2 <- (sum(z[nc]^4, na.rm=NAOK)/sum(nc))/(s2^2)
} else {
b2 <- (sum(z^4, na.rm=NAOK)/n)/(s2^2)
}
} else {
if (adjust.x) {
b2 <- (sum(z[nc]^4, na.rm=NAOK)/(sum(nc)-1))/(s2^2)
} else {
b2 <- (sum(z^4, na.rm=NAOK)/(n-1))/(s2^2)
}
}
Wi2 <- sapply(listw$weights, function(x) sum(x^2))
A <- (n-b2) / (n-1)
......
# Copyright 2001-5 by Roger Bivand
# Copyright 2001-18 by Roger Bivand
#
moran <- function(x, listw, n, S0, zero.policy=NULL, NAOK=FALSE) {
......@@ -21,7 +21,7 @@ moran <- function(x, listw, n, S0, zero.policy=NULL, NAOK=FALSE) {
moran.test <- function(x, listw, randomisation=TRUE, zero.policy=NULL,
alternative="greater", rank = FALSE, na.action=na.fail, spChk=NULL,
adjust.n=TRUE) {
adjust.n=TRUE, drop.EI2=FALSE) {
alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
"is not a listw object"))
......@@ -60,14 +60,12 @@ moran.test <- function(x, listw, randomisation=TRUE, zero.policy=NULL,
tmp <- K*(wc$S1*(wc$nn - wc$n) - 2*wc$n*wc$S2 + 6*S02)
if (tmp > VI) warning("Kurtosis overflow,\ndistribution of variable does not meet test assumptions")
VI <- (VI - tmp) / (wc$n1*wc$n2*wc$n3*S02)
tmp <- (VI - EI^2)
if (tmp < 0) warning("Negative variance,\ndistribution of variable does not meet test assumptions")
VI <- tmp
if (!drop.EI2) VI <- (VI - EI^2)
if (VI < 0) warning("Negative variance,\ndistribution of variable does not meet test assumptions")
} else {
VI <- (wc$nn*wc$S1 - wc$n*wc$S2 + 3*S02) / (S02*(wc$nn - 1))
tmp <- (VI - EI^2)
if (tmp < 0) warning("Negative variance,\ndistribution of variable does not meet test assumptions")
VI <- tmp
if (!drop.EI2) VI <- (VI - EI^2)
if (VI < 0) warning("Negative variance,\ndistribution of variable does not meet test assumptions")
}
ZI <- (I - EI) / sqrt(VI)
statistic <- ZI
......@@ -86,7 +84,10 @@ moran.test <- function(x, listw, randomisation=TRUE, zero.policy=NULL,
data.name <- paste(xname, ifelse(rank,
"using rank correction",""), "\nweights:",
wname, ifelse(is.null(na.act), "", paste("\nomitted:",
paste(na.act, collapse=", "))), "\n")
paste(na.act, collapse=", "))),
ifelse(adjust.n && isTRUE(any(sum(card(listw$neighbours) == 0L))),
"n reduced by no-neighbour observations\n", ""),
ifelse(drop.EI2, "EI^2 term dropped in VI", ""), "\n")
res <- list(statistic=statistic, p.value=PrI, estimate=vec,
alternative=alternative, method=method, data.name=data.name)
if (!is.null(na.act)) attr(res, "na.action") <- na.act
......
......@@ -31,7 +31,10 @@ read.gwt2nb <- function(file, region.id=NULL) {
odij <- read.table(file, skip=1)
# convert region.id to order
regodij <- match(odij[,1], region.id)
# 7 anmolter 2018-01-26
stopifnot(!anyNA(regodij))
regddij <- match(odij[,2], region.id)
stopifnot(!anyNA(regddij))
odij <- cbind(regodij, regddij, odij[,3])
qorder <- order(odij[,1],odij[,2])
odij <- odij[qorder,]
......
No preview for this file type
No preview for this file type
No preview for this file type
......@@ -257,7 +257,7 @@ if (dothis) {
</code></pre>
<pre><code>## GDAL GDAL_with_GEOS PROJ.4 sp
## &quot;2.2.2&quot; &quot;TRUE&quot; &quot;4.9.3&quot; &quot;1.2-5&quot;
## &quot;2.2.3&quot; &quot;TRUE&quot; &quot;5.0.0&quot; &quot;1.2-8&quot;
</code></pre>
<pre><code class="r">library(rgdal)
......
......@@ -258,8 +258,8 @@ st_write(sf_bna2_utm18, &quot;NY8_bna_utm18.gpkg&quot;)
if (dothis) sf_extSoftVersion()
</code></pre>
<pre><code>## GEOS GDAL proj.4 lwgeom GDAL_with_GEOS
## &quot;3.6.2&quot; &quot;2.2.2&quot; &quot;4.9.3&quot; NA &quot;true&quot;
<pre><code>## GEOS GDAL proj.4 GDAL_with_GEOS
## &quot;3.6.2&quot; &quot;2.2.3&quot; &quot;5.0.0&quot; &quot;true&quot;
</code></pre>
<pre><code class="r">if (!suppressPackageStartupMessages(require(rgdal, quietly=TRUE))) {
......@@ -273,7 +273,7 @@ if (dothis) {
</code></pre>
<pre><code>## GDAL GDAL_with_GEOS PROJ.4 sp
## &quot;2.2.2&quot; &quot;TRUE&quot; &quot;4.9.3&quot; &quot;1.2-5&quot;
## &quot;2.2.3&quot; &quot;TRUE&quot; &quot;5.0.0&quot; &quot;1.2-8&quot;
</code></pre>
<pre><code class="r">if (!suppressPackageStartupMessages(require(rgeos, quietly=TRUE))) {
......@@ -321,7 +321,7 @@ system.time(for(i in 1:reps) NY8_sp_1_nb &lt;- poly2nb(NY8_sp, queen=TRUE, snap=
</code></pre>
<pre><code>## user system elapsed
## 0.0578 0.0002 0.0582
## 0.0574 0.0002 0.0578
</code></pre>
<pre><code class="r">NY8_sp_1_nb
......@@ -341,7 +341,7 @@ system.time(for(i in 1:reps) NY8_sp_1_fB_nb &lt;- poly2nb(NY8_sp, queen=TRUE, sn
</code></pre>
<pre><code>## user system elapsed
## 0.0207 0.0001 0.0209
## 0.0209 0.0000 0.0211
</code></pre>
<pre><code class="r">all.equal(NY8_sp_1_fB_nb, NY8_sp_1_nb, check.attributes=FALSE)
......@@ -361,7 +361,7 @@ system.time(for(i in 1:reps) NY8_sp_1_fB_nb &lt;- poly2nb(NY8_sp, queen=TRUE, sn
</code></pre>
<pre><code>## user system elapsed
## 0.4699 0.0000 0.4710
## 0.4565 0.0000 0.4575
</code></pre>
<p>As we can see, the sf-based contiguity test is an order of magnitude slower than spdep::poly2nb; fortunately, it also scales linearly in the number of observations. spdep::poly2nb uses two heuristics, first to find candidate neighbours from intersecting bounding boxes, and second to use the symmetry of the relationship to halve the number of remaining tests. This means that performance is linear in n, but with overhead for identifying candidates, and back-filling symmetric neighbours. Further, spdep::poly2nb stops searching for queen contiguity as soon as the first neighbour point is found within snap distance (if not identical, which is tested first); second neighbour point for rook contiguities. spdep::poly2nb was heavily optimised when written, as processor speed was a major constraint at that time. </p>
......@@ -405,7 +405,7 @@ all.equal(NY8_sf_1_nb, NY8_sp_1_nb, check.attributes=FALSE)
</code></pre>
<pre><code>## user system elapsed
## 0.0511 0.0002 0.0515
## 0.0505 0.0000 0.0507
</code></pre>
<p>It is the use of GEOS functionality that costs time, as we can see by using rgeos::gTouches:</p>
......@@ -414,7 +414,7 @@ all.equal(NY8_sf_1_nb, NY8_sp_1_nb, check.attributes=FALSE)
</code></pre>
<pre><code>## user system elapsed
## 0.4591 0.0000 0.4614
## 0.4604 0.0000 0.4613
</code></pre>
<pre><code class="r">all.equal(sp_touches, NY8_sp_1_nb, check.attributes=FALSE)
......@@ -429,7 +429,7 @@ all.equal(NY8_sf_1_nb, NY8_sp_1_nb, check.attributes=FALSE)
</code></pre>
<pre><code>## user system elapsed
## 0.4580 0.0000 0.4591
## 0.4606 0.0000 0.4616
</code></pre>
<pre><code class="r">class(NY8_sf_1_touch) &lt;- &quot;sgbp&quot;
......@@ -579,7 +579,7 @@ all.equal(NY8_sf_old_1_nb_buf, NY8_sf_1_nb, check.attributes=FALSE)
</code></pre>
<pre><code>## user system elapsed
## 7e-04 0e+00 7e-04
## 8e-04 0e+00 8e-04
</code></pre>
<p>&#39;row.names&#39; gives the IDs (from feature FIDs if read with &#39;rgdal::readOGR&#39;):</p>
......@@ -601,7 +601,7 @@ is.projected(NY8_sp)
</code></pre>
<pre><code>## user system elapsed
## 0.5000 0.0000 0.5014
## 0.0250 0.0000 0.0252
</code></pre>
<pre><code class="r">if (!isTRUE(all.equal(st_crs(NY8_pt_sf), st_crs(NY8_sf)))) st_crs(NY8_pt_sf) &lt;- st_crs(NY8_sf)
......@@ -634,7 +634,7 @@ if (zm %in% c(&quot;XYZ&quot;))
</code></pre>
<pre><code>## user system elapsed
## 0.001 0.000 0.001
## 0.0010 0.0000 0.0011
</code></pre>
<p>Unfortunately, the coordinate matrices differ:</p>
......@@ -680,7 +680,7 @@ points(coords_sf[diffs,], pch=4, col=&quot;blue&quot;)
</code></pre>
<pre><code>## user system elapsed
## 0.0056 0.0000 0.0055
## 0.0055 0.0000 0.0055
</code></pre>