invIrM {spdep} | R Documentation |
Computes the matrix used for generating simultaneous autoregressive random variables, for a given value of rho, a neighbours list object, a chosen coding scheme style, and optionally a list of general weights corresponding to neighbours.
invIrM(neighbours, rho, glist=NULL, style="W", method="solve", feasible=NULL) invIrW(listw, rho, method="solve", feasible=NULL) powerWeights(W, rho, order=250, X, tol=.Machine$double.eps^(3/5))
neighbours |
an object of class |
rho |
autoregressive parameter |
glist |
list of general weights corresponding to neighbours |
style |
|
method |
default |
feasible |
if NULL, the given value of rho is checked to see if it lies within its feasible range, if TRUE, the test is not conducted |
listw |
a |
W |
A spatial weights matrix (either a dense matrix or a CsparseMatrix) |
order |
Power series maximum limit |
X |
A numerical matrix |
tol |
Tolerance for convergence of power series |
The invIrW
function generates the full weights matrix V, checks that rho lies in its feasible range between 1/min(eigen(V)) and 1/max(eigen(V)), and returns the nxn inverted matrix
(I - ρ V)^{-1}
. With method=“chol”, Cholesky decomposition is used, thanks to contributed code by Markus Reder and Werner Mueller.
The powerWeights
function uses power series summation to cumulate the product
(I - ρ V)^{-1} \%*\% X
from
(I + ρ V + (ρ V)^2 + …) \%*\% X
, which can be done by storing only sparse V and several matrices of the same dimensions as X. This makes it possible to handle larger spatial weights matrices.
An nxn matrix with a "call" attribute; the powerWeights
function returns a matrix of the same dimensions as X which has been multipled by the power series equivalent of the dense matrix
(I - ρ V)^{-1}
.
Roger Bivand Roger.Bivand@nhh.no
Tiefelsdorf, M., Griffith, D. A., Boots, B. 1999 A variance-stabilizing coding scheme for spatial link matrices, Environment and Planning A, 31, pp. 165-180; Tiefelsdorf, M. 2000 Modelling spatial processes, Lecture notes in earth sciences, Springer, p. 76; Haining, R. 1990 Spatial data analysis in the social and environmental sciences, Cambridge University Press, p. 117; Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 152; Reder, M. and Mueller, W. (2007) An Improvement of the invIrM Routine of the Geostatistical R-package spdep by Cholesky Inversion, Statistical Projects, LV No: 238.205, SS 2006, Department of Applied Statistics, Johannes Kepler University, Linz
nb7rt <- cell2nb(7, 7, torus=TRUE) set.seed(1) x <- matrix(rnorm(500*length(nb7rt)), nrow=length(nb7rt)) res0 <- apply(invIrM(nb7rt, rho=0.0, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) res2 <- apply(invIrM(nb7rt, rho=0.2, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) res4 <- apply(invIrM(nb7rt, rho=0.4, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) res6 <- apply(invIrM(nb7rt, rho=0.6, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) res8 <- apply(invIrM(nb7rt, rho=0.8, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) res9 <- apply(invIrM(nb7rt, rho=0.9, method="chol", feasible=TRUE) %*% x, 2, function(x) var(x)/length(x)) plot(density(res9), col="red", xlim=c(-0.01, max(density(res9)$x)), ylim=range(density(res0)$y), xlab="estimated variance of the mean", main=expression(paste("Effects of spatial autocorrelation for different ", rho, " values"))) lines(density(res0), col="black") lines(density(res2), col="brown") lines(density(res4), col="green") lines(density(res6), col="orange") lines(density(res8), col="pink") legend(c(-0.02, 0.01), c(7, 25), legend=c("0.0", "0.2", "0.4", "0.6", "0.8", "0.9"), col=c("black", "brown", "green", "orange", "pink", "red"), lty=1, bty="n") ## Not run: x <- matrix(rnorm(length(nb7rt)), ncol=1) system.time(e <- invIrM(nb7rt, rho=0.9, method="chol", feasible=TRUE) %*% x) system.time(e <- invIrM(nb7rt, rho=0.9, method="chol", feasible=NULL) %*% x) system.time(e <- invIrM(nb7rt, rho=0.9, method="solve", feasible=TRUE) %*% x) system.time(e <- invIrM(nb7rt, rho=0.9, method="solve", feasible=NULL) %*% x) W <- as(as_dgRMatrix_listw(nb2listw(nb7rt)), "CsparseMatrix") system.time(ee <- powerWeights(W, rho=0.9, X=x)) all.equal(e, as(ee, "matrix"), check.attributes=FALSE) nb60rt <- cell2nb(60, 60, torus=TRUE) W <- as(as_dgRMatrix_listw(nb2listw(nb60rt)), "CsparseMatrix") set.seed(1) x <- matrix(rnorm(dim(W)[1]), ncol=1) system.time(ee <- powerWeights(W, rho=0.3, X=x)) str(as(ee, "matrix")) obj <- errorsarlm(as(ee, "matrix")[,1] ~ 1, listw=nb2listw(nb60rt), method="Matrix") coefficients(obj) ## End(Not run)