execute defunct

Author: rsbivand

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: spatialreg
-Version: 1.1-9
-Date: 2021-09-28
+Version: 1.2-1
+Date: 2021-11-10
 Title: Spatial Regression Analysis
 Encoding: UTF-8
 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "[email protected]", comment=c(ORCID="0000-0003-2392-6140")),
@@ -15,6 +15,7 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bi
              person("Rein", "Halbersma", role = "ctb"),
              person("James", "LeSage", role = "ctb"),
              person("Angela", "Li", role = "ctb"),
+	     person("Hongfei", "Li", role = "ctb"),
              person("Jielai", "Ma", role = "ctb"),
              person("Abhirup", "Mallik", role = c("ctb", "trl")),
              person("Giovanni", "Millo", role = "ctb"),
@@ -24,11 +25,11 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bi
              person(given = "Mauricio", family = "Sarrias", role = c("ctb"), email = "[email protected]"),
              person(given = "JuanTomas", family = "Sayago", role = c("ctb"), email = "[email protected]"),
              person("Michael", "Tiefelsdorf", role = "ctb"))
-Depends: R (>= 3.3.0), spData, Matrix
+Depends: R (>= 3.3.0), spData, Matrix, sf
 Imports: spdep, expm, coda, methods, MASS, boot, splines, LearnBayes,
         nlme, gmodels
-Suggests: parallel, RSpectra, sf, tmap, foreign, spam, knitr, lmtest,
-        sandwich, rmarkdown
+Suggests: parallel, RSpectra, tmap, foreign, spam, knitr, lmtest,
+        sandwich, rmarkdown, igraph
 Description: A collection of all the estimation functions for spatial cross-sectional models (on lattice/areal data using spatial weights matrices) contained up to now in 'spdep', 'sphet' and 'spse'. These model fitting functions include maximum likelihood methods for cross-sectional models proposed by 'Cliff' and 'Ord' (1973, ISBN:0850860369) and (1981, ISBN:0850860814), fitting methods initially described by 'Ord' (1975) <doi:10.1080/01621459.1975.10480272>. The models are further described by 'Anselin' (1988) <doi:10.1007/978-94-015-7799-1>. Spatial two stage least squares and spatial general method of moment models initially proposed by 'Kelejian' and 'Prucha' (1998) <doi:10.1023/A:1007707430416> and (1999) <doi:10.1111/1468-2354.00027> are provided. Impact methods and MCMC fitting methods proposed by 'LeSage' and 'Pace' (2009) <doi:10.1201/9781420064254> are implemented for the family of cross-sectional spatial regression models. Methods for fitting the log determinant term in maximum likelihood and MCMC fitting are compared by 'Bivand et al.' (2013) <doi:10.1111/gean.12008>, and model fitting methods by 'Bivand' and 'Piras' (2015) <doi:10.18637/jss.v063.i18>; both of these articles include extensive lists of references. 'spatialreg' >= 1.1-* correspond to 'spdep' >= 1.1-1, in which the model fitting functions are deprecated and pass through to 'spatialreg', but will mask those in 'spatialreg'. From versions 1.2-*, the functions will be made defunct in 'spdep'.
 License: GPL-2
 URL: https://github.com/r-spatial/spatialreg/, https://r-spatial.github.io/spatialreg/

NAMESPACE

@@ -15,10 +15,11 @@ importFrom("stats", "lm.fit", "as.formula", "coef", "lm", "model.extract",
 
 importFrom("spdep", "is.symmetric.glist", "lag.listw", "card",
  "nb2listw", "n.comp.nb", "listw2U", "listw2mat", "Szero", "listw2sn",
- "invIrW", "mat2listw", "is.symmetric.nb", "chkIDs", "get.spChkOption",
+ "mat2listw", "is.symmetric.nb", "chkIDs", "get.spChkOption",
  "subset.listw")
 
 import(Matrix, except=c("expm"))
+import(sf)
 importFrom(expm, expAtv, expm)
 importFrom(coda, as.mcmc, HPDinterval)
 importFrom(MASS, mvrnorm)
@@ -150,7 +151,7 @@ exportMethods(coerce)
 export(as_dgRMatrix_listw, as_dsTMatrix_listw, as_dsCMatrix_I,
 	as_dsCMatrix_IrW, Jacobian_W, listw2U_Matrix)
 
-export(powerWeights)
+export(powerWeights, invIrW, invIrM)
 
 export(as.spam.listw, listw2U_spam)
 
@@ -175,4 +176,6 @@ export(set.mcOption, get.mcOption, set.coresOption, get.coresOption,
 export(set.VerboseOption, get.VerboseOption, set.ZeroPolicyOption, 
         get.ZeroPolicyOption)
 
+export(aple.plot, localAple, aple, aple.mc)
+
 

NEWS.md

@@ -1,7 +1,9 @@
-# Version 1.1-9 (development)
+# Version 1.2-1 (development)
 
 * Add Fortran character handling USE_FC_LEN_T WRE §6.6.1
 
+* Add **spdep** split-out functionality
+
 # Version 1.1-8 (2021-05-03)
 
 * #18 standardize use of `coef()` methods for (some) fitted model summary objects

R/AAA.R

@@ -8,7 +8,7 @@ assign("verbose", spdep::get.VerboseOption(), envir = .spatialregOptions)
 assign("mc", spdep::get.mcOption(), envir = .spatialregOptions)
 assign("cores", spdep::get.coresOption(), envir = .spatialregOptions)
 assign("cluster", spdep::get.ClusterOption(), envir = .spatialregOptions)
-#assign("rlecuyerSeed", rep(12345, 6), envir = .spatialregOptions)
+assign("rlecuyerSeed", rep(12345, 6), envir = .spatialregOptions)
 #assign("listw_is_CsparseMatrix", FALSE, envir = .spatialregOptions)
 
 setOldClass(c("listw"))

R/aple.R

@@ -0,0 +1,46 @@
+preAple <- function(x, listw, override_similarity_check=FALSE, useTrace=TRUE) {
+    stopifnot(isTRUE(all.equal(mean(x), 0.0)))
+    stopifnot(is.vector(x))
+#    if (!requireNamespace("spatialreg", quietly=TRUE))
+#        stop("install the spatialreg package")
+    W <- as(listw, "CsparseMatrix")
+    n <- dim(W)[1]
+    if (useTrace) {
+        trWW <- sum(diag(W %*% W))
+    } else {
+        if (listw$style %in% c("W", "S") && !override_similarity_check) {
+            can.sim <- can.be.simmed(listw)
+            eig <- eigenw(similar.listw(listw))
+        } else {
+            can.sim <- FALSE
+            eig <- eigenw(listw)
+        }
+        if (is.complex(eig)) trWW <- Re(crossprod(eig))
+        else trWW <- crossprod(eig)
+# modified 110414 RSB
+#            eig <- Re(eig)
+#        trWW <- crossprod(eig)
+    }
+    corterm <- (trWW/n) * Diagonal(n)
+    corterm <- as(corterm, "CsparseMatrix")
+    WU <- ((W + t(W))/2)
+    W2 <- crossprod(W) + corterm
+    res <- list(W=W, corterm=corterm, W2=W2, WU=WU, n=n)
+    res
+}
+
+inAple <- function(x, pre) {
+    xwx <- crossprod(x, (pre$WU %*% x))
+    xwwx <- crossprod(x, (pre$W2 %*% x))
+    res <- c(as.matrix(xwx/xwwx))
+    res
+}
+
+aple <- function(x, listw, override_similarity_check=FALSE, useTrace=TRUE) {
+    pre <- preAple(x=x, listw=listw,
+        override_similarity_check=override_similarity_check, useTrace=useTrace)
+    res <- inAple(x=x, pre=pre)
+    res
+}
+
+

R/aple.mc.R

@@ -0,0 +1,60 @@
+aple.mc <- function(x, listw, nsim, override_similarity_check=FALSE,
+    useTrace=TRUE) {
+    aple.boot <- function(var, i, ...) {
+        var <- var[i]
+        return(inAple(x=var, ...))
+    }
+    pre <- preAple(x=x, listw=listw,
+        override_similarity_check=override_similarity_check, useTrace=useTrace)
+    
+    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 <- get.ClusterOption()
+        if (is.null(cl)) {
+            parallel <- "no"
+            warning("no cluster in ClusterOption, parallel set to no")
+        }
+    }
+    res <- boot(x, statistic=aple.boot, R=nsim, sim="permutation", pre=pre,
+        parallel=parallel, ncpus=ncpus, cl=cl)
+
+    res
+}
+
+boot_wrapper_in <- function(cl, nsim) {
+      if (requireNamespace("parallel", quietly = TRUE)) {
+#        require(rlecuyer)
+        rlseed <- get("rlecuyerSeed", envir = .spatialregOptions)
+        if (storage.mode(rlseed) != "integer") rlseed <- as.integer(rlseed)
+        if (length(rlseed) != 6L) rlseed <- rep(12345L, 6)
+        parallel::clusterSetRNGStream(cl, iseed=rlseed)
+        parallel::clusterEvalQ(cl, library(spdep))
+        nnsim <- ceiling(nsim/length(cl))
+        nnsim
+      } else {
+        stop("parallel not available")
+      }
+}
+
+boot_wrapper_out <- function(lres, mcall) {
+        res <- list()
+        res$t0 <- lres[[1]]$t0
+        res$t <- matrix(c(sapply(lres, function(x) x$t)), ncol=1)
+        res$R <- sum(sapply(lres, function(x) x$R))
+        res$data <- lres[[1]]$data
+        res$seed <- c(sapply(lres, function(x) x$seed))
+        res$statistic <- lres[[1]]$statistic
+        res$sim <- lres[[1]]$sim
+        res$call <- mcall
+        res$stype <- lres[[1]]$stype
+        res$strata <- lres[[1]]$strata
+        class(res) <- "boot"
+        res
+}

R/apleplot.R

@@ -0,0 +1,22 @@
+aple.plot <- function(x, listw, override_similarity_check=FALSE, useTrace=TRUE,
+    do.plot=TRUE, ...) {
+    pre <- preAple(x=x, listw=listw,
+        override_similarity_check=override_similarity_check, useTrace=useTrace)
+    W2e <- eigen(pre$W2)
+    SQRTW2 <- W2e$vectors %*% (diag(W2e$values^(0.5)) %*% t(W2e$vectors))
+    X <- drop(SQRTW2 %*% x)
+    NSQRTW2 <- W2e$vectors %*% (diag(W2e$values^(-0.5)) %*% t(W2e$vectors))
+    Y <- drop(NSQRTW2 %*% pre$WU %*% x)
+    if (do.plot) {
+        plot(X, Y, ...)
+    }
+    list(X=X, Y=Y)
+}
+
+localAple <- function(x, listw, override_similarity_check=FALSE, useTrace=TRUE) {
+    aplepl <- aple.plot(x, listw,
+        override_similarity_check=override_similarity_check,
+        useTrace=useTrace, do.plot=FALSE)
+    res <- (length(aplepl$X) * aplepl$Y * aplepl$X) / c(crossprod(aplepl$X))
+    res
+}

R/sparse_mat.R

@@ -114,3 +114,48 @@ powerWeights <- function(W, rho, order=250, X, tol=.Machine$double.eps^(3/5)) {
     acc
 }
 
+invIrM <- function(neighbours, rho, glist=NULL, style="W", method="solve", 
+	feasible=NULL) {
+	if(class(neighbours) != "nb") stop("Not a neighbours list")
+	invIrW(nb2listw(neighbours, glist=glist, style=style), rho=rho, 
+		method=method, feasible=feasible)
+}
+
+invIrW <- function(x, rho, method="solve", feasible=NULL) {
+	if(inherits(x, "listw")) {
+	  n <- length(x$neighbours)
+	  V <- listw2mat(x)
+	} else if (inherits(x, "Matrix") || inherits(x, "matrix")) {
+	  if (method == "chol" && all(t(x) == x))
+            stop("No Cholesky method for matrix or sparse matrix object")
+          n <- dim(x)[1]
+          V <- x
+	} else stop("Not a weights list or a Sparse Matrix")
+	if (is.null(feasible) || (is.logical(feasible) && !feasible)) {
+		e <- eigen(V, only.values = TRUE)$values
+		if (is.complex(e)) feasible <- 1/(range(Re(e)))
+		else feasible <- 1/(range(e))
+		if (rho <= feasible[1] || rho >= feasible[2])
+			stop(paste("Rho", rho, "outside feasible range:",
+                        paste(feasible, collapse=":")))
+	}
+        if (method == "chol"){
+	    if (x$style %in% c("W", "S") && !(spatialreg::can.be.simmed(x)))
+			stop("Cholesky method requires symmetric weights")
+		if (x$style %in% c("B", "C", "U") && 
+			!(is.symmetric.glist(x$neighbours, x$weights)))
+			stop("Cholesky method requires symmetric weights")
+		if (x$style %in% c("W", "S")) {
+			V <- listw2mat(listw2U(spatialreg::similar.listw(x)))
+		}
+		mat <- diag(n) - rho * V
+		res <- chol2inv(chol(mat))
+	} else if (method == "solve") {
+		mat <- diag(n) - rho * V
+		res <- solve(mat)
+	} else stop("unknown method")
+	attr(res, "call") <- match.call()
+	res
+}
+
+

man/aple.Rd

@@ -0,0 +1,50 @@
+\name{aple}
+\alias{aple}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Approximate profile-likelihood estimator (APLE)}
+\description{
+  The Approximate profile-likelihood estimator (APLE) of the simultaneous autoregressive model's spatial dependence parameter was introduced in Li et al. (2007). It employs a correction term using the eigenvalues of the spatial weights matrix, and consequently should not be used for large numbers of observations. It also requires that the variable has a mean of zero, and it is assumed that it has been detrended. The spatial weights object is assumed to be row-standardised, that is using default \code{style="W"} in \code{nb2listw}.
+}
+\usage{
+aple(x, listw, override_similarity_check=FALSE, useTrace=TRUE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{a zero-mean detrended continuous variable}
+  \item{listw}{a \code{listw} object from for example \code{nb2listw}}
+  \item{override\_similarity\_check}{default FALSE, if TRUE - typically for row-standardised weights with asymmetric underlying general weights - similarity is not checked}
+  \item{useTrace}{default TRUE, use trace of sparse matrix \code{W \%*\% W} (Li et al. (2010)), if FALSE, use crossproduct of eigenvalues of \code{W} as in Li et al. (2007)}
+}
+\details{
+ This implementation has been checked with Hongfei Li's own implementation using her data; her help was very valuable.
+}
+\value{
+  A scalar APLE value.
+}
+\references{Li, H, Calder, C. A. and Cressie N. A. C. (2007) Beyond Moran's I: testing for spatial dependence based on the spatial autoregressive model. Geographical Analysis 39, 357-375; Li, H, Calder, C. A. and Cressie N. A. C. (2012) One-step estimation of spatial dependence parameters: Properties and extensions of the APLE statistic, Journal of Multivariate Analysis 105, 68-84.}
+\author{Roger Bivand \email{[email protected]}}
+
+\seealso{\code{\link{nb2listw}}, \code{\link{aple.mc}}, \code{\link{aple.plot}}
+}
+\examples{
+wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
+library(spdep)
+nbr1 <- poly2nb(wheat, queen=FALSE)
+nbrl <- nblag(nbr1, 2)
+nbr12 <- nblag_cumul(nbrl)
+cms0 <- with(as.data.frame(wheat), tapply(yield, c, median))
+cms1 <- c(model.matrix(~ factor(c) -1, data=wheat) \%*\% cms0)
+wheat$yield_detrend <- wheat$yield - cms1
+isTRUE(all.equal(c(with(as.data.frame(wheat),
+ tapply(yield_detrend, c, median))), rep(0.0, 25),
+ check.attributes=FALSE))
+moran.test(wheat$yield_detrend, nb2listw(nbr12, style="W"))
+aple(as.vector(scale(wheat$yield_detrend, scale=FALSE)), nb2listw(nbr12, style="W"))
+\dontrun{
+errorsarlm(yield_detrend ~ 1, wheat, nb2listw(nbr12, style="W"))
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{spatial}
+

man/aple.mc.Rd

@@ -0,0 +1,72 @@
+\name{aple.mc}
+\alias{aple.mc}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Approximate profile-likelihood estimator (APLE) permutation test}
+\description{
+  A permutation bootstrap test for the approximate profile-likelihood estimator (APLE).
+}
+\usage{
+aple.mc(x, listw, nsim, override_similarity_check=FALSE, useTrace=TRUE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{a zero-mean detrended continuous variable}
+  \item{listw}{a \code{listw} object from for example \code{nb2listw}}
+  \item{nsim}{number of simulations}
+  \item{override\_similarity\_check}{default FALSE, if TRUE - typically for row-standardised weights with asymmetric underlying general weights - similarity is not checked}
+  \item{useTrace}{default TRUE, use trace of sparse matrix \code{W \%*\% W} (Li et al. (2010)), if FALSE, use crossproduct of eigenvalues of \code{W} as in Li et al. (2007)}
+}
+
+\value{
+  A \code{boot} object as returned by the \code{boot} function.
+}
+\references{Li, H, Calder, C. A. and Cressie N. A. C. (2007) Beyond Moran's I: testing for spatial dependence based on the spatial autoregressive model. Geographical Analysis 39, 357-375; Li, H, Calder, C. A. and Cressie N. A. C. (2012) One-step estimation of spatial dependence parameters: Properties and extensions of the APLE statistic, Journal of Multivariate Analysis 105, 68-84.}
+\author{Roger Bivand \email{[email protected]}}
+
+\seealso{\code{\link{aple}}, \code{\link[boot]{boot}}}
+\examples{
+\dontrun{
+wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
+nbr1 <- poly2nb(wheat, queen=FALSE)
+nbrl <- nblag(nbr1, 2)
+nbr12 <- nblag_cumul(nbrl)
+wheat_g <- wheat
+st_geometry(wheat_g) <- NULL
+cms0 <- with(wheat_g, tapply(yield, c, median))
+cms1 <- c(model.matrix(~ factor(c) -1, data=wheat) \%*\% cms0)
+wheat$yield_detrend <- wheat$yield - cms1
+oldRNG <- RNGkind()
+RNGkind("L'Ecuyer-CMRG")
+set.seed(1L)
+boot_out_ser <- aple.mc(as.vector(scale(wheat$yield_detrend, scale=FALSE)),
+ nb2listw(nbr12, style="W"), nsim=500)
+plot(boot_out_ser)
+boot_out_ser
+library(parallel)
+oldCores <- set.coresOption(NULL)
+nc <- detectCores(logical=FALSE)
+# set nc to 1L here
+if (nc > 1L) nc <- 1L
+invisible(set.coresOption(nc))
+set.seed(1L)
+if (!get.mcOption()) {
+  cl <- makeCluster(nc)
+  set.ClusterOption(cl)
+} else{
+  mc.reset.stream()
+}
+boot_out_par <- aple.mc(as.vector(scale(wheat$yield_detrend, scale=FALSE)),
+    nb2listw(nbr12, style="W"), nsim=500)
+if (!get.mcOption()) {
+  set.ClusterOption(NULL)
+  stopCluster(cl)
+}
+boot_out_par
+invisible(set.coresOption(oldCores))
+RNGkind(oldRNG[1], oldRNG[2])
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{spatial}
+