diff -Nru r-cran-spdep-0.7-9+dfsg/debian/changelog r-cran-spdep-0.8-1+dfsg/debian/changelog --- r-cran-spdep-0.7-9+dfsg/debian/changelog 2018-10-05 19:50:00.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/debian/changelog 2018-11-26 08:30:23.000000000 +0000 @@ -1,3 +1,9 @@ +r-cran-spdep (0.8-1+dfsg-1) unstable; urgency=medium + + * New upstream version + + -- Andreas Tille Mon, 26 Nov 2018 09:30:23 +0100 + r-cran-spdep (0.7-9+dfsg-1) unstable; urgency=medium * Team upload. diff -Nru r-cran-spdep-0.7-9+dfsg/DESCRIPTION r-cran-spdep-0.8-1+dfsg/DESCRIPTION --- r-cran-spdep-0.7-9+dfsg/DESCRIPTION 2018-09-28 11:40:03.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/DESCRIPTION 2018-11-21 11:30:03.000000000 +0000 @@ -1,6 +1,6 @@ Package: spdep -Version: 0.7-9 -Date: 2018-08-31 +Version: 0.8-1 +Date: 2018-11-21 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"), @@ -60,7 +60,7 @@ License: GPL (>= 2) VignetteBuilder: knitr NeedsCompilation: yes -Packaged: 2018-09-28 10:40:02 UTC; rsb +Packaged: 2018-11-21 09:58:45 UTC; rsb Author: Roger Bivand [cre, aut] (), Micah Altman [ctb], Luc Anselin [ctb], @@ -94,4 +94,4 @@ Danlin Yu [ctb] Maintainer: Roger Bivand Repository: CRAN -Date/Publication: 2018-09-28 11:40:03 UTC +Date/Publication: 2018-11-21 11:30:03 UTC Binary files /tmp/tmpA7Ul0C/j0LBvJ1zDU/r-cran-spdep-0.7-9+dfsg/inst/doc/CO69.pdf and /tmp/tmpA7Ul0C/csjUDQMPPe/r-cran-spdep-0.8-1+dfsg/inst/doc/CO69.pdf differ Binary files /tmp/tmpA7Ul0C/j0LBvJ1zDU/r-cran-spdep-0.7-9+dfsg/inst/doc/nb.pdf and /tmp/tmpA7Ul0C/csjUDQMPPe/r-cran-spdep-0.8-1+dfsg/inst/doc/nb.pdf differ Binary files /tmp/tmpA7Ul0C/j0LBvJ1zDU/r-cran-spdep-0.7-9+dfsg/inst/doc/sids.pdf and /tmp/tmpA7Ul0C/csjUDQMPPe/r-cran-spdep-0.8-1+dfsg/inst/doc/sids.pdf differ Binary files /tmp/tmpA7Ul0C/j0LBvJ1zDU/r-cran-spdep-0.7-9+dfsg/inst/doc/SpatialFiltering.pdf and /tmp/tmpA7Ul0C/csjUDQMPPe/r-cran-spdep-0.8-1+dfsg/inst/doc/SpatialFiltering.pdf differ diff -Nru r-cran-spdep-0.7-9+dfsg/man/cell2nb.Rd r-cran-spdep-0.8-1+dfsg/man/cell2nb.Rd --- r-cran-spdep-0.7-9+dfsg/man/cell2nb.Rd 2017-11-01 16:24:24.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/man/cell2nb.Rd 2018-11-16 08:49:19.000000000 +0000 @@ -14,7 +14,7 @@ grid is mapped onto a torus, removing edge effects. } \usage{ -cell2nb(nrow, ncol, type="rook", torus=FALSE) +cell2nb(nrow, ncol, type="rook", torus=FALSE, legacy=FALSE) mrc2vi(rowcol, nrow, ncol) rookcell(rowcol, nrow, ncol, torus=FALSE, rmin=1, cmin=1) queencell(rowcol, nrow, ncol, torus=FALSE, rmin=1, cmin=1) @@ -25,6 +25,7 @@ \item{ncol}{number of columns in the grid} \item{type}{rook or queen} \item{torus}{map grid onto torus} + \item{legacy}{default FALSE, nrow/ncol reversed, if TRUE wrong col/row directions (see \url{https://github.com/r-spatial/spdep/issues/20})} \item{rowcol}{matrix with two columns of row, column indices} \item{i}{vector of vector indices corresponding to rowcol} \item{rmin}{lowest row index} @@ -46,5 +47,18 @@ plot(nb7rt, xy) nb7rt <- cell2nb(7, 7, torus=TRUE) summary(nb7rt) +# https://github.com/r-spatial/spdep/issues/20 +GT <- GridTopology(c(1, 1), c(1,1), c(10, 50)) +SPix <- as(SpatialGrid(GT), "SpatialPixels") +nb_rook_cont <- poly2nb(as(SPix, "SpatialPolygons"), queen=FALSE) +nb_rook_dist <- dnearneigh(coordinates(SPix), 0, 1.01) +all.equal(nb_rook_cont, nb_rook_dist, check.attributes=FALSE) +## [1] TRUE +t.nb <- cell2nb(nrow=50, ncol=10, type='rook', legacy=TRUE) +isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) +## [1] FALSE +t.nb <- cell2nb(nrow=50, ncol=10, type='rook') +isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) +## [1] TRUE } \keyword{spatial} diff -Nru r-cran-spdep-0.7-9+dfsg/man/errorsarlm.Rd r-cran-spdep-0.8-1+dfsg/man/errorsarlm.Rd --- r-cran-spdep-0.7-9+dfsg/man/errorsarlm.Rd 2018-09-27 11:35:20.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/man/errorsarlm.Rd 2018-11-20 20:19:41.000000000 +0000 @@ -1,8 +1,9 @@ -% Copyright 2002-16 by Roger S. Bivand +% Copyright 2002-18 by Roger S. Bivand \name{errorsarlm} \alias{errorsarlm} \alias{lmSLX} \alias{create_WX} +\alias{spBreg_err} %\alias{sar.error.f} %\alias{sar.error.f.s} \title{Spatial simultaneous autoregressive error model estimation} @@ -18,6 +19,8 @@ errorsarlm(formula, data=list(), listw, na.action, weights=NULL, Durbin, etype, method="eigen", quiet=NULL, zero.policy=NULL, interval = NULL, tol.solve=1.0e-10, trs=NULL, control=list()) +spBreg_err(formula, data = list(), listw, na.action, Durbin, etype, + zero.policy=NULL, control=list()) lmSLX(formula, data = list(), listw, na.action, weights=NULL, Durbin=TRUE, zero.policy=NULL) create_WX(x, listw, zero.policy=NULL, prefix="") @@ -96,6 +99,34 @@ \item{pre_eig}{default NULL; may be used to pass a pre-computed vector of eigenvalues} }} +\section{Extra Bayesian control arguments}{ +\describe{ + \item{ldet_method}{default \dQuote{SE_classic}; equivalent to the \code{method} argument in \code{lagsarlm}} + \item{interval}{default \code{c(-1, 1)}; used unmodified or set internally by \code{jacobianSetup}} + \item{ndraw}{default \code{2500L}; integer total number of draws} + \item{nomit}{default \code{500L}; integer total number of omitted burn-in draws} + \item{thin}{default \code{1L}; integer thinning proportion} + \item{verbose}{default \code{FALSE}; inverse of \code{quiet} argument in \code{lagsarlm}} + \item{detval}{default \code{NULL}; not yet in use, precomputed matrix of log determinants} + \item{prior}{a list with the following components: + \describe{ + \item{lambdaMH}{default FALSE; if TRUE use Metropolis or FALSE griddy Gibbs} + \item{Tbeta}{default \code{NULL}; values of the betas variance-covariance matrix, set to \code{diag(k)*1e+12} if \code{NULL}} + \item{c_beta}{default \code{NULL}; values of the betas set to 0 if \code{NULL}} + \item{rho}{default \code{0.5}; value of the autoregressive coefficient} + \item{sige}{default \code{1}; value of the residual variance} + \item{nu}{default \code{0}; informative Gamma(nu,d0) prior on sige} + \item{d0}{default \code{0}; informative Gamma(nu,d0) prior on sige} + \item{a1}{default \code{1.01}; parameter for beta(a1,a2) prior on rho} + \item{a2}{default \code{1.01}; parameter for beta(a1,a2) prior on rho} + \item{cc}{default \code{0.2}; initial tuning parameter for M-H sampling} + \item{gG_sige}{default TRUE; include sige in lambda griddy Gibbs update} + }} +} + +} + + \value{ A list object of class \code{sarlm} \item{type}{"error"} @@ -171,15 +202,15 @@ \examples{ data(oldcol) lw <- nb2listw(COL.nb, style="W") +ev <- eigenw(similar.listw(lw)) COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", quiet=FALSE) + lw, quiet=FALSE, control=list(pre_eig=ev)) summary(COL.errW.eig) -ev <- eigenw(similar.listw(lw)) COL.errW.eig_ev <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", control=list(pre_eig=ev)) + lw, control=list(pre_eig=ev)) all.equal(coefficients(COL.errW.eig), coefficients(COL.errW.eig_ev)) COL.errB.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - nb2listw(COL.nb, style="B"), method="eigen", quiet=FALSE) + nb2listw(COL.nb, style="B")) summary(COL.errB.eig) W <- as(nb2listw(COL.nb), "CsparseMatrix") trMatc <- trW(W, type="mult") @@ -187,13 +218,13 @@ lw, method="Matrix", quiet=FALSE, trs=trMatc) summary(COL.errW.M) COL.SDEM.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", etype="emixed") + lw, etype="emixed", control=list(pre_eig=ev)) summary(COL.SDEM.eig) COL.SDEM.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", Durbin=TRUE) + lw, Durbin=TRUE, control=list(pre_eig=ev)) summary(COL.SDEM.eig) COL.SDEM.eig <- errorsarlm(CRIME ~ DISCBD + INC + HOVAL, data=COL.OLD, - lw, method="eigen", Durbin=~INC) + lw, Durbin=~INC, control=list(pre_eig=ev)) summary(COL.SDEM.eig) summary(impacts(COL.SDEM.eig)) COL.SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw) @@ -237,45 +268,83 @@ resid(COL.err.NA) \donttest{ lw <- nb2listw(COL.nb, style="W") -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen")) +print(system.time(ev <- eigenw(similar.listw(lw)))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="eigen", control=list(pre_eig=ev)))) ocoef <- coefficients(COL.errW.eig) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", control=list(LAPACK=FALSE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="eigen", control=list(compiled_sse=TRUE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix_J", control=list(super=TRUE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix_J", control=list(super=FALSE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix_J", control=list(super=as.logical(NA)))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix", control=list(super=TRUE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix", control=list(super=FALSE))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="Matrix", control=list(super=as.logical(NA)))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="spam", control=list(spamPivot="MMD"))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="spam", control=list(spamPivot="RCM"))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="spam_update", control=list(spamPivot="MMD"))) -all.equal(ocoef, coefficients(COL.errW.eig)) -system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, - lw, method="spam_update", control=list(spamPivot="RCM"))) -all.equal(ocoef, coefficients(COL.errW.eig)) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="eigen", control=list(pre_eig=ev, LAPACK=FALSE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="eigen", control=list(pre_eig=ev, compiled_sse=TRUE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix_J", control=list(super=TRUE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix_J", control=list(super=FALSE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix_J", control=list(super=as.logical(NA))))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix", control=list(super=TRUE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix", control=list(super=FALSE)))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="Matrix", control=list(super=as.logical(NA))))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="spam", control=list(spamPivot="MMD")))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="spam", control=list(spamPivot="RCM")))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="spam_update", control=list(spamPivot="MMD")))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, + lw, method="spam_update", control=list(spamPivot="RCM")))) +print(all.equal(ocoef, coefficients(COL.errW.eig))) +} +\dontrun{ +if (require(coda, quietly=TRUE)) { +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw) +print(summary(COL.err.Bayes)) +print(raftery.diag(COL.err.Bayes, r=0.01)) +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, + control=list(prior=list(lambdaMH=TRUE))) +print(summary(COL.err.Bayes)) +print(raftery.diag(COL.err.Bayes, r=0.01)) +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, + Durbin=TRUE) +print(summary(COL.err.Bayes)) +print(summary(impacts(COL.err.Bayes))) +print(raftery.diag(COL.err.Bayes, r=0.01)) +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, + Durbin=TRUE, control=list(prior=list(lambdaMH=TRUE))) +print(summary(COL.err.Bayes)) +print(summary(impacts(COL.err.Bayes))) +print(raftery.diag(COL.err.Bayes, r=0.01)) +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, + Durbin=~INC) +print(summary(COL.err.Bayes)) +print(summary(impacts(COL.err.Bayes))) +print(raftery.diag(COL.err.Bayes, r=0.01)) +set.seed(1) +COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, + Durbin=~INC, control=list(prior=list(lambdaMH=TRUE))) +print(summary(COL.err.Bayes)) +print(summary(impacts(COL.err.Bayes))) +print(raftery.diag(COL.err.Bayes, r=0.01)) +} } } \keyword{spatial} diff -Nru r-cran-spdep-0.7-9+dfsg/man/impacts.sarlm.Rd r-cran-spdep-0.8-1+dfsg/man/impacts.sarlm.Rd --- r-cran-spdep-0.7-9+dfsg/man/impacts.sarlm.Rd 2018-09-27 11:28:16.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/man/impacts.sarlm.Rd 2018-11-20 20:19:41.000000000 +0000 @@ -6,6 +6,8 @@ \alias{impacts.gmsar} \alias{impacts.SLX} \alias{impacts.MCMC_sar_g} +\alias{impacts.MCMC_sem_g} +\alias{impacts.MCMC_sac_g} \alias{impacts.lagmess} \alias{plot.lagImpact} \alias{print.lagImpact} @@ -29,6 +31,8 @@ empirical=FALSE) \method{impacts}{SLX}(obj, ...) \method{impacts}{MCMC_sar_g}(obj, ..., tr=NULL, listw=NULL, evalues=NULL, Q=NULL) +\method{impacts}{MCMC_sem_g}(obj, ..., tr=NULL, listw=NULL, evalues=NULL, Q=NULL) +\method{impacts}{MCMC_sac_g}(obj, ..., tr=NULL, listw=NULL, evalues=NULL, Q=NULL) \method{plot}{lagImpact}(x, ..., choice="direct", trace=FALSE, density=TRUE) \method{print}{lagImpact}(x, ..., reportQ=NULL) \method{summary}{lagImpact}(object, ..., zstats=FALSE, short=FALSE, reportQ=NULL) diff -Nru r-cran-spdep-0.7-9+dfsg/man/lagsarlm.Rd r-cran-spdep-0.8-1+dfsg/man/lagsarlm.Rd --- r-cran-spdep-0.7-9+dfsg/man/lagsarlm.Rd 2018-09-27 12:07:03.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/man/lagsarlm.Rd 2018-10-26 09:18:32.000000000 +0000 @@ -98,6 +98,7 @@ \item{detval}{default \code{NULL}; not yet in use, precomputed matrix of log determinants} \item{prior}{a list with the following components: \describe{ + \item{rhoMH}{default FALSE; use Metropolis or griddy Gibbs} \item{Tbeta}{default \code{NULL}; values of the betas variance-covariance matrix, set to \code{diag(k)*1e+12} if \code{NULL}} \item{c_beta}{default \code{NULL}; values of the betas set to 0 if \code{NULL}} \item{rho}{default \code{0.5}; value of the autoregressive coefficient} diff -Nru r-cran-spdep-0.7-9+dfsg/man/sacsarlm.Rd r-cran-spdep-0.8-1+dfsg/man/sacsarlm.Rd --- r-cran-spdep-0.7-9+dfsg/man/sacsarlm.Rd 2018-09-21 08:23:05.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/man/sacsarlm.Rd 2018-11-20 20:19:41.000000000 +0000 @@ -1,6 +1,7 @@ % Copyright 2010 by Roger S. Bivand \name{sacsarlm} \alias{sacsarlm} +\alias{spBreg_sac} %- Also NEED an '\alias' for EACH other topic documented here. \title{Spatial simultaneous autoregressive SAC model estimation} \description{ @@ -15,6 +16,8 @@ method = "eigen", quiet = NULL, zero.policy = NULL, tol.solve = 1e-10, llprof=NULL, interval1=NULL, interval2=NULL, trs1=NULL, trs2=NULL, control = list()) +spBreg_sac(formula, data = list(), listw, listw2=NULL, na.action, + Durbin, type, zero.policy=NULL, control=list()) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -58,6 +61,31 @@ \item{pre_eig1, pre_eig2}{default NULL; may be used to pass pre-computed vectors of eigenvalues} }} +\section{Extra Bayesian control arguments}{ +\describe{ + \item{ldet_method}{default \dQuote{SE_classic}; equivalent to the \code{method} argument in \code{lagsarlm}} + \item{interval1}{default \code{c(-1, 1)}; used unmodified or set internally by \code{jacobianSetup}} + \item{interval2}{default \code{c(-1, 1)}; used unmodified or set internally by \code{jacobianSetup}} + \item{ndraw}{default \code{2500L}; integer total number of draws} + \item{nomit}{default \code{500L}; integer total number of omitted burn-in draws} + \item{thin}{default \code{1L}; integer thinning proportion} + \item{detval1}{default \code{NULL}; not yet in use, precomputed matrix of log determinants} + \item{detval2}{default \code{NULL}; not yet in use, precomputed matrix of log determinants} + \item{prior}{a list with the following components: + \describe{ + \item{Tbeta}{default \code{NULL}; values of the betas variance-covariance matrix, set to \code{diag(k)*1e+12} if \code{NULL}} + \item{c_beta}{default \code{NULL}; values of the betas set to 0 if \code{NULL}} + \item{lambda}{default \code{0.5}; value of the autoregressive coefficient} + \item{rho}{default \code{0.5}; value of the autoregressive coefficient} + \item{sige}{default \code{1}; value of the residual variance} + \item{nu}{default \code{0}; informative Gamma(nu,d0) prior on sige} + \item{d0}{default \code{0}; informative Gamma(nu,d0) prior on sige} + \item{a1}{default \code{1.01}; parameter for beta(a1,a2) prior on rho} + \item{a2}{default \code{1.01}; parameter for beta(a1,a2) prior on rho} + \item{cc1}{default \code{0.2}; initial tuning parameter for M-H sampling} + \item{cc2}{default \code{0.2}; initial tuning parameter for M-H sampling} + }} +}} \value{ A list object of class \code{sarlm} @@ -121,7 +149,21 @@ summary(COL.sacW.eig) W <- as(listw, "CsparseMatrix") trMatc <- trW(W, type="mult") +set.seed(1) summary(impacts(COL.sacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) +\dontrun{ +if (require(coda, quietly=TRUE)) { +set.seed(1) +COL.sacW.B0 <- spBreg_sac(CRIME ~ INC + HOVAL, data=COL.OLD, listw, + Durbin=FALSE, control=list(ndraw=1500L, nomit=500L)) +print(summary(COL.sacW.B0)) +print(summary(impacts(COL.sacW.B0, tr=trMatc), zstats=TRUE, short=TRUE)) +set.seed(1) +COL.sacW.B1 <- spBreg_sac(CRIME ~ INC + HOVAL, data=COL.OLD, listw, + Durbin=TRUE, control=list(ndraw=1500L, nomit=500L)) +print(summary(COL.sacW.B1)) +print(summary(impacts(COL.sacW.B1, tr=trMatc), zstats=TRUE, short=TRUE)) +}} COL.msacW.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, type="sacmixed", control=list(pre_eig1=ev, pre_eig2=ev)) summary(COL.msacW.eig) diff -Nru r-cran-spdep-0.7-9+dfsg/MD5 r-cran-spdep-0.8-1+dfsg/MD5 --- r-cran-spdep-0.7-9+dfsg/MD5 2018-09-28 11:40:04.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/MD5 2018-11-21 11:30:03.000000000 +0000 @@ -1,7 +1,7 @@ 2428ff9c8fc1c5390035b6f09bcc9b13 *ChangeLog -1c571035c97af4a0dcecb04d74a28eb1 *DESCRIPTION -66845965ba6a0335e7f5047a2e6335b0 *NAMESPACE -ca3af94aff175d2e462ffb33bd9ac568 *R/AAA.R +cbbbd754df9ad22ad5a3dadf116a2cb2 *DESCRIPTION +09a2ab4c31094ce5610e513919eefed6 *NAMESPACE +b4e50cc6e10db793dad9495b12ec0b47 *R/AAA.R b0e1b90e65dd7b62bdd3948361cb5237 *R/EBI.R 8010b6a2e0a46f013ad55ec541812c8e *R/LOSH.R 94be548ba7c37f33ef8076b74533092f *R/LOSH.cs.R @@ -17,7 +17,7 @@ de45513cf2ca3a506f98a26f3327dc52 *R/autocov.R a2d046225a488a388bac40c86e791658 *R/bptest.sarlm.R fe7171c29d8d6c25abf04f41b6c41623 *R/card.R -f852d0ee75f87db10d90f57a0b9664ca *R/cellneighbours.R +29a5fb402eadaae2d252438f73fb60b4 *R/cellneighbours.R 4007c384b0abdb7e46000b6552b540a3 *R/choynowski.R 848467dae136cea24da9c34e544954c2 *R/components.R e64d4fcc81be06647777c3c71fc125b3 *R/cyclical.R @@ -32,7 +32,7 @@ 5c7fc104be2f03783a4193fbf2351b19 *R/geary.R 494162cf871ebd36494283ccde5dc1df *R/globalG.R d6320c56e705e0f9e81e46a851a74fa0 *R/graph2nb.R -ba11cc59fe221a88ec891917557e9c69 *R/impacts.R +f0e03ce89a6d35746eb64787fddd5e51 *R/impacts.R d37880c837d66c5167434da6858c1c40 *R/jacobian.R f842c1af5dcab63b5f1eb0bd3099341a *R/jacobian_setup.R e9042b5d26783b2224f601afca856bbf *R/jc.R @@ -83,13 +83,13 @@ 64633a828f14988e6d3aa97d6ce8ea2d *R/relneigh.R 40534663fb14a40fc659768038487304 *R/rotation.R b8ae91f8e8af0ad2ddb651840f879a3d *R/s2sls.R -7e345218e6398e5541bdb6795b7542ed *R/sacsarlm.R +c8830dd2e3bc8b2be50a8d5935331d7c *R/sacsarlm.R f37ded41d4bc7b2dbb9eb97797234328 *R/sarlm.R 18ef3fcc6a1af79afbf767716ad411c0 *R/skater.R dfd3c7d74224e215cc90db24b16d7948 *R/soi.R e73d855ad4dc75409201fbf6c36b43f9 *R/sp.correlogram.R f6e0b551534ecbd40f91583a3fee8e0f *R/sp.mantel.R -3dfced1679d79034afe0e8dabb86305e *R/spBreg.R +f83123988526a32f0cc314e069a7cba3 *R/spBreg.R 93c531e3c39a5ccd1506fe92a95fa9c0 *R/spChkOption.R 05f368670f558994e9886590eb057c61 *R/spautolm.R 73870459b0081b9975e674cd64bebcc9 *R/subset.nb.R @@ -97,7 +97,7 @@ c54d106f445c7133868cb332cc0334c6 *R/summary.spsarlm.R 700f61092fae8a0b1fce15adf29df47a *R/tolerance.nb.R 79530fa1378ad438ea0695b44a3b134f *R/tri2nb.R -e4c7a7d4e336723499fe98a9e13fd67f *R/utils.R +10aacfda12b849f39c221bc5207bfda9 *R/utils.R 6f1fcd3c6b9caaf4d4e78e0ad8028046 *R/weights-utils.R 3b6aba89dced1b9ea9d9ae9bac66e772 *build/vignette.rds b818650dc1a3e0335a7767fe78f721f5 *data/columbus.rda @@ -108,24 +108,24 @@ 2428ff9c8fc1c5390035b6f09bcc9b13 *inst/ChangeLog 99c4d523114b0e2769464b65a9cf1150 *inst/doc/CO69.R 709cb9502be2c6176e03747da963bc15 *inst/doc/CO69.Rnw -0f24c04dd8b59dcd04bc527241b0cb20 *inst/doc/CO69.pdf +6979458dfc43e04a3963b1cf208f29e7 *inst/doc/CO69.pdf 27fb91633d64bff1470c754abdc09ac1 *inst/doc/SpatialFiltering.R 00a3b786362437b0d02c578910457a5e *inst/doc/SpatialFiltering.Rnw -0d0662ca2e60dffec16fe94ca0494832 *inst/doc/SpatialFiltering.pdf +51509ae4e67f9c56b82bdf2aa527b836 *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 -5bdc643639aa64546687b037d2700f57 *inst/doc/nb.pdf +e9a07622ceb2de0600def52ce3e4765a *inst/doc/nb.pdf 52619007cbaa4f7cb7e64423347d720c *inst/doc/nb_igraph.R 0a4d2c611361f8e9242b91d09a7d7354 *inst/doc/nb_igraph.Rmd -874e81a99d52041f5118675d3d1a9cec *inst/doc/nb_igraph.html +4b906dce62a4ea71617365b2d2423b59 *inst/doc/nb_igraph.html 522b2484fbd34290eaa545c19089dd77 *inst/doc/nb_sf.R 941aa1849341d87ce66e36c168b8ad7f *inst/doc/nb_sf.Rmd -5c29dfe4b8be1fed33517153a4a37176 *inst/doc/nb_sf.html +67416923a86a0cdba072e2bfdb83f6fd *inst/doc/nb_sf.html 98ba2b15ccd8b39b51dfd25ed6a25d15 *inst/doc/sids.R af6bf80e9a912760fa14ff55f82d2d00 *inst/doc/sids.Rnw -9297523cbff884d9d5bbedb86d48733c *inst/doc/sids.pdf +2fc2ce634f0ec95aaec3a58a8235d693 *inst/doc/sids.pdf 53178ac43dd150bbb835b63e018f0cd4 *inst/etc/misc/geary_eire.txt e3fb11232e6cb27770d15aaf4f632aa4 *inst/etc/misc/nc_xtra.dbf 2591ac95bd9191d0aee07450f7caf9bb *inst/etc/misc/nyadjwts.dbf @@ -169,7 +169,7 @@ 81ff6bdff475ba95c31840696eb785b1 *man/bhicv.Rd f8e83f726753614ab62c204ee2806e1c *man/bptest.sarlm.Rd 8dd642b512505da779f1e2ec058d3b85 *man/card.Rd -cdf723571f5373da8abfdd7b23b9be2e *man/cell2nb.Rd +779f1768dc606712d7189da4c07efcf3 *man/cell2nb.Rd 897c57c4006f615c5976ffe5cc378d1c *man/choynowski.Rd 4ba5ab6de2b94cb5eecbc6a6c4dd8c6b *man/columbus.Rd d4954f89422d430aa1a10601925127dd *man/compon.Rd @@ -180,14 +180,14 @@ e557ef79c38a6b08b541b17b9d3ad190 *man/edit.nb.Rd 10cfcfcc90c5486a46f3f532a58e3d2d *man/eigenw.Rd c28e1263be5dbd60f5da17b124efc534 *man/eire.Rd -93f8bf62dbe3680450a7f0e588889f0d *man/errorsarlm.Rd +7d0f3653af77add34dc971bfa409ab11 *man/errorsarlm.Rd 25e1a53cd6d9cbbb7bea0c1fea3e1ef0 *man/geary.Rd e5e5f9840a79014d86c7758de589b0f9 *man/geary.mc.Rd 7958dfcf21f532bbf6b75486556fce52 *man/geary.test.Rd dbc6c20ca81659ee9531af9074299f5d *man/globalG.test.Rd e876fa32337f7b7e9dbe701e84dc33c3 *man/graphneigh.Rd fe5102da5d3290cdbb47395805511ec2 *man/gstsls.Rd -6c77656fbeefa246777a5c490338406f *man/impacts.sarlm.Rd +bcce3ae95d41a9155c9e2873b4e1c072 *man/impacts.sarlm.Rd 3ebb36e6b2d2b8aea7d122dcfdddb527 *man/include.self.Rd 9a65838851fd6b5e78cc853b3eab8dd1 *man/invIrM.Rd 3a5aa076bb3c12d7c608a5d57eae9fb7 *man/joincount.mc.Rd @@ -197,7 +197,7 @@ 3666112819e0fb6937f4684147755588 *man/knn2nb.Rd 530a137b1bad26cbbc658a52e2f48562 *man/lag.listw.Rd 3527e506449803f2057451cd4366509a *man/lagmess.Rd -9160fa3d78ee90f1e21abe9d0785b71f *man/lagsarlm.Rd +c1e15075d4bfe68dad1cb2e15e38c664 *man/lagsarlm.Rd c4859675fde07731b661ffd7452f996b *man/lee.Rd 50357b77a1351040444263892230e680 *man/lee.mc.Rd 92b00a952642c858fcf2fd9115381d89 *man/lee.test.Rd @@ -240,7 +240,7 @@ fe86bf830f291a16d60102c52bfdf01e *man/read.gwt2nb.Rd 22587da32f258cb04ac771791ff2a5f4 *man/residuals.sarlm.Rd e3f28cd38c69d847efdf6143af31637d *man/rotation.Rd -77299c2be09abb48984cb02ee5dc93e0 *man/sacsarlm.Rd +9704b9a2833c0e18bdfa6dad671d562f *man/sacsarlm.Rd d4ec4879830f577619f8194ac46d2876 *man/set.mcOption.Rd 1ca46f773740c2b0f9ae03b6f0a81388 *man/set.spChkOption.Rd 923c93d9b47db0ba93a8d3156e3fcff0 *man/similar.listw.Rd diff -Nru r-cran-spdep-0.7-9+dfsg/NAMESPACE r-cran-spdep-0.8-1+dfsg/NAMESPACE --- r-cran-spdep-0.7-9+dfsg/NAMESPACE 2018-03-08 12:11:45.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/NAMESPACE 2018-11-20 20:19:41.000000000 +0000 @@ -3,6 +3,7 @@ import(spData) importFrom(stats, model.matrix, model.response, aggregate, optimize, optim, nlminb, influence.measures, lag, punif, na.fail, as.formula, terms, lm, model.extract, coefficients, residuals, var, deviance, integrate, summary.lm, coef, pchisq, logLik, gaussian, pnorm, model.weights, model.offset, glm.fit, sd, predict, formula, lm.fit, anova, ppois, model.frame, AIC, qnorm, napredict, nlm, optimHess, resid, fitted, weighted.residuals, weights, naresid, rnorm, optimise, uniroot, mahalanobis, dist, p.adjust, density, fitted.values, delete.response, printCoefmat, quantile, na.omit, cor, predict.lm) importFrom("stats", "dbeta", "rchisq", "runif") +importFrom("stats", "spline") importFrom(nlme, fdHess) importFrom(deldir, deldir) importFrom(boot, boot) @@ -112,7 +113,7 @@ export(lee.mc, lee, lee.test) -export(spBreg_lag) +export(spBreg_lag, spBreg_err, spBreg_sac) export(LOSH, LOSH.cs, LOSH.mc) @@ -172,6 +173,8 @@ S3method(MCMCsamp, sarlm) S3method(impacts, MCMC_sar_g) +S3method(impacts, MCMC_sem_g) +S3method(impacts, MCMC_sac_g) S3method(impacts, SLX) S3method(predict, SLX) diff -Nru r-cran-spdep-0.7-9+dfsg/R/AAA.R r-cran-spdep-0.8-1+dfsg/R/AAA.R --- r-cran-spdep-0.7-9+dfsg/R/AAA.R 2017-11-01 16:24:24.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/AAA.R 2018-10-12 09:04:32.000000000 +0000 @@ -12,6 +12,8 @@ assign("rlecuyerSeed", rep(12345, 6), envir = .spdepOptions) assign("listw_is_CsparseMatrix", FALSE, envir = .spdepOptions) +setOldClass(c("listw")) + #.conflicts.OK <- TRUE #.onLoad <- function(lib, pkg) { diff -Nru r-cran-spdep-0.7-9+dfsg/R/cellneighbours.R r-cran-spdep-0.8-1+dfsg/R/cellneighbours.R --- r-cran-spdep-0.7-9+dfsg/R/cellneighbours.R 2017-11-01 16:24:24.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/cellneighbours.R 2018-11-16 08:49:19.000000000 +0000 @@ -71,11 +71,11 @@ res } -cell2nb <- function(nrow, ncol, type="rook", torus=FALSE) { +cell2nb <- function(nrow, ncol, type="rook", torus=FALSE, legacy=FALSE) { nrow <- as.integer(nrow) if (nrow < 1) stop("nrow nonpositive") ncol <- as.integer(ncol) - if (ncol < 1) stop("nrow nonpositive") + if (ncol < 1) stop("ncol nonpositive") xcell <- NULL if (type == "rook") xcell <- rookcell if (type == "queen") xcell <- queencell @@ -85,11 +85,19 @@ if (n < 0) stop("non-positive number of cells") res <- vector(mode="list", length=n) rownames <- character(n) - for (i in 1:n) { + if (legacy) { + for (i in 1:n) { res[[i]] <- sort(mrc2vi(xcell(vi2mrc(i, nrow, ncol), nrow, ncol, torus), nrow, ncol)) rownames[i] <- paste(vi2mrc(i, nrow, ncol), collapse=":") - } + } + } else { + for (i in 1:n) { + res[[i]] <- sort(mrc2vi(xcell(vi2mrc(i, ncol, nrow), + ncol, nrow, torus), ncol, nrow)) + rownames[i] <- paste(vi2mrc(i, ncol, nrow), collapse=":") + } + } class(res) <- "nb" attr(res, "call") <- match.call() attr(res, "region.id") <- rownames diff -Nru r-cran-spdep-0.7-9+dfsg/R/impacts.R r-cran-spdep-0.8-1+dfsg/R/impacts.R --- r-cran-spdep-0.7-9+dfsg/R/impacts.R 2018-09-27 11:28:16.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/impacts.R 2018-10-31 12:58:54.000000000 +0000 @@ -128,16 +128,16 @@ stopifnot(!is.null(attr(obj, "mixedImps"))) n <- nrow(obj$model) k <- obj$qr$rank - impactsWX(attr(obj, "mixedImps"), n, k, type="SLX") + impactsWX(attr(obj, "mixedImps"), n, k, type="SLX", method="estimable") } impactSDEM <- function(obj) { n <- nrow(obj$tarX) k <- ncol(obj$tarX) - impactsWX(obj$emixedImps, n, k, type="SDEM") + impactsWX(obj$emixedImps, n, k, type="SDEM", method="estimable") } -impactsWX <- function(obj, n, k, type="SLX") { +impactsWX <- function(obj, n, k, type="SLX", method="estimable") { imps <- lapply(obj, function(x) x[, 1]) names(imps) <- c("direct", "indirect", "total") attr(imps, "bnames") <- rownames(obj[[1]]) @@ -148,7 +148,7 @@ attr(res, "n") <- n attr(res, "k") <- k attr(res, "type") <- type - attr(res, "method") <- "estimable" + attr(res, "method") <- method attr(res, "bnames") <- rownames(obj[[1]]) class(res) <- "WXImpact" res @@ -822,15 +822,20 @@ if (zstats) { # 100928 Eelke Folmer if (length(attr(object, "bnames")) == 1L) { + semat <- sapply(lres, function(x) x$statistics[2]) + semat <- matrix(semat, ncol=3) + colnames(semat) <- c("Direct", "Indirect", "Total") zmat <- sapply(lres, function(x) x$statistics[1]/x$statistics[2]) zmat <- matrix(zmat, ncol=3) colnames(zmat) <- c("Direct", "Indirect", "Total") } else { + semat <- sapply(lres, function(x) x$statistics[,2]) + colnames(semat) <- c("Direct", "Indirect", "Total") zmat <- sapply(lres, function(x) x$statistics[,1]/x$statistics[,2]) colnames(zmat) <- c("Direct", "Indirect", "Total") } pzmat <- 2*(1-pnorm(abs(zmat))) - res <- c(res, list(zmat=zmat, pzmat=pzmat)) + res <- c(res, list(semat=semat, zmat=zmat, pzmat=pzmat)) if (!is.null(Qmcmc) && !is.null(reportQ) && reportQ) { Qzmats <- lapply(Qmcmc, function(x) { Qm <- matrix(x$statistics[,1]/x$statistics[,2], @@ -928,7 +933,11 @@ } if (!is.null(x$zmat)) { cat("========================================================\n") - cat("Simulated z-values:\n") + cat("Simulated standard errors\n") + mat <- x$semat + rownames(mat) <- attr(x, "bnames") + print(mat) + cat("\nSimulated z-values:\n") mat <- x$zmat rownames(mat) <- attr(x, "bnames") print(mat) diff -Nru r-cran-spdep-0.7-9+dfsg/R/sacsarlm.R r-cran-spdep-0.8-1+dfsg/R/sacsarlm.R --- r-cran-spdep-0.7-9+dfsg/R/sacsarlm.R 2018-09-27 11:28:16.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/sacsarlm.R 2018-11-19 08:34:16.000000000 +0000 @@ -124,8 +124,15 @@ cn <- names(aliased) names(aliased) <- substr(cn, 2, nchar(cn)) if (any(aliased)) { - nacoef <- which(aliased) + if (is.formula(Durbin)) { + stop("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + } else { + warning("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + nacoef <- which(aliased) x <- x[,-nacoef] + } } LL_null_lm <- NULL if ("(Intercept)" %in% colnames(x)) LL_null_lm <- logLik(lm(y ~ 1)) diff -Nru r-cran-spdep-0.7-9+dfsg/R/spBreg.R r-cran-spdep-0.8-1+dfsg/R/spBreg.R --- r-cran-spdep-0.7-9+dfsg/R/spBreg.R 2018-09-27 11:28:16.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/spBreg.R 2018-11-21 09:48:38.000000000 +0000 @@ -12,8 +12,9 @@ in_coef=0.1, type="MC", correct=TRUE, trunc=TRUE, SE_method="LU", nrho=200, interpn=2000, SElndet=NULL, LU_order=FALSE, pre_eig=NULL, interval=c(-1, 1), ndraw=2500L, nomit=500L, thin=1L, - verbose=FALSE, detval=NULL, prior=list(Tbeta=NULL, c_beta=NULL, - rho=0.5, sige=1, nu=0, d0=0, a1 = 1.01, a2 = 1.01)) + verbose=FALSE, detval=NULL, prior=list(rhoMH=FALSE, Tbeta=NULL, + c_beta=NULL, rho=0.5, sige=1, nu=0, d0=0, a1 = 1.01, a2 = 1.01, + cc = 0.2, c=NULL, T=NULL)) priors <- con$prior nmsP <- names(priors) priors[(namp <- names(control$prior))] <- control$prior @@ -181,6 +182,8 @@ psave <- numeric(con$ndraw) ssave <- numeric(con$ndraw) lsave <- numeric(con$ndraw) + acc_rate <- NULL + if (priors$rhoMH) acc_rate <- numeric(con$ndraw) #% ====== initializations #% compute this stuff once to save time @@ -194,6 +197,15 @@ nu1 = n + 2*priors$nu nrho = length(detval1) rho_out = 0 + gsize = detval1[2] - detval1[1] + acc = 0 + noninf <- TRUE + if (!is.null(priors$c) && !is.null(priors$T)) { + if (length(priors$c) == 1L && is.numeric(priors$c) && + length(priors$T) == 1L && is.numeric(priors$T)) noninf <- FALSE + } + nmk = (n-k)/2 + cc <- priors$cc # nano_1 = 0 # nano_2 = 0 # nano_3 = 0 @@ -222,53 +234,107 @@ chi = rchisq(1, nu1) #chi = chis_rnd(1,nu1); sige = as.numeric(d1/chi) # see eq 5.30, p. 141 ssave[iter] = as.vector(sige) -# nano_2 <- nano_2 + microbenchmark::get_nanotime() - nano_p - + + if (priors$rhoMH) { + +## % metropolis step to get rho update + i1 = max(which(detval1 <= (rho + gsize))) + i2 = max(which(detval1 <= (rho - gsize))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 + detm = detval2[index] + if (noninf) { + epe = (crossprod(e))/(2*sige) + } else { + epe = (crossprod(e))/(2*sige) + 0.5*(((rho-priors$c)^2)/(priors$T*sige)) + } + rhox = detm - epe + accept = 0L; + rho2 = rho + cc*rnorm(1); + while(accept == 0L) { + if((rho2 > con$interval[1]) & (rho2 < con$interval[2])) { + accept=1 + } else { + rho2 = rho + cc*rnorm(1) + } + } + i1 = max(which(detval1 <= (rho2 + gsize))) + i2 = max(which(detval1 <= (rho2 - gsize))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 + detm = detval2[index] + yss = y - rho2*wy + e = yss - x %*% bhat + if (noninf) { + epe = (crossprod(e))/(2*sige) + } else { + epe = (crossprod(e))/(2*sige) + + 0.5*(((rho-priors$c)^2)/(priors$T*sige)) + } + rhoy = detm - epe + if ((rhoy - rhox) > exp(1)) { + p = 1 + } else { + ratio = exp(rhoy-rhox) + p = min(1, ratio) + } + ru = runif(1) + if(ru < p) { + rho = rho2 + acc = acc + 1 + } + acc_rate[iter] = acc/iter + if(acc_rate[iter] < 0.4) cc=cc/1.1 + if(acc_rate[iter] > 0.6) cc=cc*1.1 + + + AI = solve((xpx + sige*TI)) + b0 = AI %*% (xpy + sige*TIc) + bd = AI %*% (xpWy + sige*TIc) + e0 = y - x%*%b0 + ed = wy - x%*%bd + epe0 = as.vector(crossprod(e0)) + eped = as.vector(crossprod(ed)) + epe0d = as.vector(crossprod(ed, e0)) + } else { ###% update rho using griddy Gibbs -# nano_p <- microbenchmark::get_nanotime() - AI = solve((xpx + sige*TI)) - b0 = AI %*% (xpy + sige*TIc) - bd = AI %*% (xpWy + sige*TIc) - e0 = y - x%*%b0 - ed = wy - x%*%bd - epe0 = as.vector(crossprod(e0)) - eped = as.vector(crossprod(ed)) - epe0d = as.vector(crossprod(ed, e0)) - nmk = (n-k)/2 - z = epe0 - 2*detval1*epe0d + detval1*detval1*eped - den = detval2 - nmk*log(z) - den = den + bprior - - nd = length(den) - adj = max(den) - den = den - adj - xd = exp(den) - - ## trapezoid rule - isum = sum((detval1[2:nd] + detval1[1:(nd-1)])*(xd[2:nd] - - xd[1:(nd-1)])/2)#VIRGILIO:FIXED - zd = abs(xd/isum) - den = cumsum(zd) - - rnd = runif(1)*sum(zd) - cond <- den <= rnd - if (any(cond)) { + AI = solve((xpx + sige*TI)) + b0 = AI %*% (xpy + sige*TIc) + bd = AI %*% (xpWy + sige*TIc) + e0 = y - x%*%b0 + ed = wy - x%*%bd + epe0 = as.vector(crossprod(e0)) + eped = as.vector(crossprod(ed)) + epe0d = as.vector(crossprod(ed, e0)) + nmk = (n-k)/2 + z = epe0 - 2*detval1*epe0d + detval1*detval1*eped + den = detval2 - nmk*log(z) + den = den + bprior + + nd = length(den) + adj = max(den) + den = den - adj + xd = exp(den) + + ## trapezoid rule + isum = sum((detval1[2:nd] + detval1[1:(nd-1)])*(xd[2:nd] + - xd[1:(nd-1)])/2)#VIRGILIO:FIXED + zd = abs(xd/isum) + den = cumsum(zd) + + rnd = runif(1)*sum(zd) + cond <- den <= rnd + if (any(cond)) { # ind = which(den <= rnd) - idraw = which.min(cond) - 1 #max(ind) + idraw = which.min(cond) - 1 #max(ind) # if (idraw > 0 & idraw < nrho) - rho = detval1[idraw]#FIXME: This sometimes fail... - } else { - rho_out = rho_out+1 + rho = detval1[idraw]#FIXME: This sometimes fail... + } else { + rho_out = rho_out+1 + } } - z = epe0 - 2*rho*epe0d + rho*rho*eped - if (idraw > 0 & idraw < nrho) - ldet = detval2[idraw] - s2 <- z/n - ll_iter <- (ldet - ((n/2)*log(2*pi)) - (n/2)*log(s2) - - (1/(2*s2))*z) psave[iter] = as.vector(rho) - lsave[iter] <- as.vector(ll_iter) # nano_3 <- nano_3 + microbenchmark::get_nanotime() - nano_p } @@ -290,8 +356,6 @@ ldet <- do_ldet(rho, env) ll_mean <- (ldet - ((n/2)*log(2*pi)) - (n/2)*log(s2) - (1/(2*s2))*sse) - lsave <- as.vector(coda::mcmc(matrix(lsave, ncol=1), start=con$nomit+1, - end=con$ndraw, thin=con$thin)[,1]) timings[["finalise"]] <- proc.time() - .ptime_start attr(res, "timings") <- do.call("rbind", timings) @@ -300,10 +364,11 @@ attr(res, "type") <- type attr(res, "rho_out") <- rho_out attr(res, "listw_style") <- listw$style - attr(res, "lsave") <- lsave attr(res, "ll_mean") <- as.vector(ll_mean) attr(res, "aliased") <- aliased + attr(res, "acc_rate") <- acc_rate attr(res, "dvars") <- dvars + attr(res, "MH") <- priors$rhoMH class(res) <- c("MCMC_sar_g", class(res)) res @@ -401,3 +466,890 @@ } + + +# translated from Matlab code sem_g.m in the Spatial Econometrics toolbox by +# James LeSage and R. Kelley Pace (http://www.spatial-econometrics.com/). +# GSoc 2011 project by Abhirup Mallik mentored by Virgilio Gómez-Rubio + +spBreg_err <- function(formula, data = list(), listw, na.action, Durbin, etype, + zero.policy=NULL, control=list()) { + timings <- list() + .ptime_start <- proc.time() +#control + con <- list(tol.opt=.Machine$double.eps^0.5, ldet_method="SE_classic", + Imult=2, cheb_q=5, MC_p=16L, MC_m=30L, super=NULL, spamPivot="MMD", + in_coef=0.1, type="MC", correct=TRUE, trunc=TRUE, LAPACK=FALSE, + SE_method="LU", nrho=200, interpn=2000, SElndet=NULL, LU_order=FALSE, + pre_eig=NULL, interval=c(-1, 1), ndraw=2500L, nomit=500L, thin=1L, + verbose=FALSE, detval=NULL, prior=list(lambdaMH=FALSE, Tbeta=NULL, + c_beta=NULL, lambda=0.5, sige=1, nu=0, d0=0, a1 = 1.01, a2 = 1.01, + cc = 0.2, gG_sige=TRUE)) + priors <- con$prior + nmsP <- names(priors) + priors[(namp <- names(control$prior))] <- control$prior + if (length(noNms <- namp[!namp %in% nmsP])) + warning("unknown names in control$prior: ", + paste(noNms, collapse = ", ")) + control$prior <- NULL + con$prior <- NULL + nmsC <- names(con) + con[(namc <- names(control))] <- control + if (length(noNms <- namc[!namc %in% nmsC])) + warning("unknown names in control: ", paste(noNms, collapse = ", ")) + stopifnot(is.logical(con$verbose)) + stopifnot(is.integer(con$ndraw)) + stopifnot(is.integer(con$nomit)) + stopifnot(is.integer(con$thin)) + stopifnot(is.logical(con$LAPACK)) + + if (is.null(zero.policy)) + zero.policy <- get.ZeroPolicyOption() + stopifnot(is.logical(zero.policy)) + if (class(formula) != "formula") formula <- as.formula(formula) + mt <- terms(formula, data = data) + mf <- lm(formula, data, na.action=na.action, method="model.frame") + na.act <- attr(mf, "na.action") + if (!inherits(listw, "listw")) stop("No neighbourhood list") + can.sim <- FALSE + if (listw$style %in% c("W", "S")) can.sim <- can.be.simmed(listw) + if (!is.null(na.act)) { + subset <- !(1:length(listw$neighbours) %in% na.act) + listw <- subset(listw, subset, zero.policy=zero.policy) + } + y <- model.extract(mf, "response") +#MatrixModels::model.Matrix() +# x <- Matrix::sparse.model.matrix(mt, mf) + x <- model.matrix(mt, mf) + n <- nrow(x) + if (n != length(listw$neighbours)) + stop("Input data and weights have different dimensions") + xcolnames <- colnames(x) + K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) + wy <- lag.listw(listw, y, zero.policy=zero.policy) + if (anyNA(wy)) stop("NAs in lagged dependent variable") +#create_WX +# check for dgCMatrix + if (missing(etype)) etype <- "error" + if (missing(Durbin)) Durbin <- ifelse(etype == "error", FALSE, TRUE) + if (listw$style != "W" && is.formula(Durbin)) { + Durbin <- TRUE + warning("formula Durbin requires row-standardised weights; set TRUE") + } + if (is.logical(Durbin) && isTRUE(Durbin)) etype <- "emixed" + if (is.formula(Durbin)) etype <- "emixed" + if (is.logical(Durbin) && !isTRUE(Durbin)) etype <- "error" + + m <- ncol(x) + xxcolnames <- colnames(x) + dvars <- c(NCOL(x), 0L) + + if (is.formula(Durbin) || isTRUE(Durbin)) { + prefix <- "lag" + if (isTRUE(Durbin)) { + WX <- create_WX(x, listw, zero.policy=zero.policy, + prefix=prefix) + } else { + dmf <- lm(Durbin, data, na.action=na.action, + method="model.frame") + fx <- try(model.matrix(Durbin, dmf), silent=TRUE) + if (class(fx) == "try-error") + stop("Durbin variable mis-match") + WX <- create_WX(fx, listw, zero.policy=zero.policy, + prefix=prefix) + inds <- match(substring(colnames(WX), 5, + nchar(colnames(WX))), colnames(x)) + if (anyNA(inds)) stop("WX variables not in X: ", + paste(substring(colnames(WX), 5, + nchar(colnames(WX)))[is.na(inds)], collapse=" ")) + icept <- grep("(Intercept)", colnames(x)) + iicept <- length(icept) > 0L + if (iicept) { + xn <- colnames(x)[-1] + } else { + xn <- colnames(x) + } + wxn <- substring(colnames(WX), nchar(prefix)+2, + nchar(colnames(WX))) + zero_fill <- NULL + if (length((which(!(xn %in% wxn)))) > 0L) + zero_fill <- length(xn) + (which(!(xn %in% wxn))) + } + dvars <- c(NCOL(x), NCOL(WX)) + if (is.formula(Durbin)) { + attr(dvars, "f") <- Durbin + attr(dvars, "inds") <- inds + attr(dvars, "zero_fill") <- zero_fill + } + x <- cbind(x, WX) + xcolnames <- colnames(x) + m <- NCOL(x) + rm(WX) + } +# x <- cbind(x, WX) +# rm(WX) +# } else if (type != "lag") stop("No such type:", type) + lm.base <- lm(y ~ x - 1) # doesn't like dgCMatrix + aliased <- is.na(coefficients(lm.base)) + cn <- names(aliased) + names(aliased) <- substr(cn, 2, nchar(cn)) + if (any(aliased)) { + if (is.formula(Durbin)) { + stop("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + } else { + warning("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + nacoef <- which(aliased) + x <- x[,-nacoef] + } + } + m <- ncol(x) + timings[["set_up"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + env <- new.env() + assign("can.sim", can.sim, envir=env) + assign("listw", listw, envir=env) + assign("similar", FALSE, envir=env) + assign("n", n, envir=env) + assign("verbose", con$verbose, envir=env) + assign("family", "SAR", envir=env) + assign("method", con$ldet_method, envir=env) + assign("LAPACK", con$LAPACK, envir=env) + W <- as(listw, "CsparseMatrix") + assign("W", W, envir=env) + + con$interval <- jacobianSetup(con$ldet_method, env, con, + pre_eig=con$pre_eig, interval=con$interval) + detval1 <- get("detval1", envir=env)[,1] + detval2 <- get("detval1", envir=env)[,2] + bprior <- dbeta(detval1, priors$a1, priors$a2) + + nm <- paste(con$ldet_method, "set_up", sep="_") + timings[[nm]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + k <- m + + if (is.null(priors$c_beta)) + priors$c_beta <- rep(0, k) + else + stopifnot(length(priors$c_beta) == k) + + if (is.null(priors$Tbeta)) + priors$Tbeta <- diag(k)*1e+12 + else + stopifnot(nrow(priors$Tbeta) == k && ncol(priors$Tbeta) == k) + + sige <- priors$sige + lambda <- priors$lambda + cc <- priors$cc + +#% storage for draws + bsave <- matrix(0, nrow=con$ndraw, ncol=k) + psave <- numeric(con$ndraw) + ssave <- numeric(con$ndraw) + acc_rate <- NULL + if (priors$lambdaMH) acc_rate <- numeric(con$ndraw) +#% ====== initializations +#% compute this stuff once to save time + + TI = solve(priors$Tbeta); # see eq 5.29, Lesage & Pace (2009) p. 140 + TIc = TI%*%priors$c_beta; + + if (m > 1 || (m == 1 && K == 1)) { + WX <- matrix(nrow=n,ncol=(m-(K-1))) + for (k in K:m) { + wx <- lag.listw(listw, x[,k], zero.policy=zero.policy) + if (any(is.na(wx))) + stop("NAs in lagged independent variable") + WX[,(k-(K-1))] <- wx + } + } + if (K == 2) { +# modified to meet other styles, email from Rein Halbersma + wx1 <- as.double(rep(1, n)) + wx <- lag.listw(listw, wx1, zero.policy=zero.policy) + if (m > 1) WX <- cbind(wx, WX) + else WX <- matrix(wx, nrow=n, ncol=1) + } + colnames(WX) <- xcolnames + rm(wx) + + +# xpx = crossprod(x) +# xpy = crossprod(x, y) +# xpWy = crossprod(x, wy) + nu1 = n + 2*priors$nu + nrho = length(detval1) + gsize = detval1[2] - detval1[1] + acc = 0 + nmk = (n-k)/2 + +# RSB Lifted out of draw_rho_sem for homoscedastic model + if (!priors$lambdaMH) { + rgrid = seq(con$interval[1]+0.01, con$interval[2]-0.01, 0.01) + ng = length(rgrid) + epet = numeric(ng) + detxt = numeric(ng) + for(i in 1:ng) { + xs = x - rgrid[i]*WX + ys = y - rgrid[i]*wy + xsxs <- crossprod(xs) + AI = solve(xsxs) + bs = AI %*% crossprod(xs, ys) + e = ys - xs%*%bs + epet[i] = crossprod(e) + detxt[i] = det(xsxs) + } + EPE = exp((spline(x=rgrid, y=log(epet), xout=detval1))$y) + DETX = exp(spline(x=rgrid, y=log(detxt), xout=detval1)$y) + } + + timings[["complete_setup"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + for (iter in 1:con$ndraw) { #% start sampling; + + #% update beta + xss = x - lambda*WX + AI = solve(crossprod(xss) + sige*TI) + yss = y - lambda*wy + b = crossprod(xss, yss) + sige*TIc + b0 = AI %*% b + bhat = MASS::mvrnorm(1, b0, sige*AI) + bsave[iter, 1:k] = as.vector(bhat) + + #update sige: + e = yss - xss %*% bhat + ed = e - lambda*lag.listw(listw, e, zero.policy=zero.policy) + d1 = 2*priors$d0 + crossprod(ed) + chi = rchisq(1, nu1) + sige = as.numeric(d1/chi) + ssave[iter] = as.vector(sige) + + if (priors$lambdaMH) { + #update lambda using M-H + i1 = max(which(detval1 <= (lambda + gsize))) + i2 = max(which(detval1 <= (lambda - gsize))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 #Fixed this + detm = detval2[index] + epe = (crossprod(e))/(2*sige) + detx = 0.5*log(det(crossprod(xss))) + rhox = detm - detx - epe + accept = 0L; + lambda2 = lambda + cc*rnorm(1); + while(accept == 0L) { + if((lambda2 > con$interval[1]) & (lambda2 < con$interval[2])) { + accept=1 + } else { + lambda2 = lambda + cc*rnorm(1) + } + } + i1 = max(which(detval1 <= (lambda2 + gsize))) + i2 = max(which(detval1 <= (lambda2 - gsize))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 #Fixed this + detm = detval2[index] + xss = x - lambda2*WX + yss = y - lambda2*wy + detx = 0.5*log(det(crossprod(xss))) + e = yss - xss %*% bhat + epe = (crossprod(e))/(2*sige) + rhoy = detm - detx - epe + if ((rhoy - rhox) > exp(1)) { + p = 1 + } else { + ratio = exp(rhoy-rhox) + p = min(1, ratio) + } + ru = runif(1) + if(ru < p) { + lambda = lambda2 + acc = acc + 1 + } + acc_rate[iter] = acc/iter + if(acc_rate[iter] < 0.4) cc=cc/1.1 + if(acc_rate[iter] > 0.6) cc=cc*1.1 + } else { +# RSB updated sige added + if (priors$gG_sige) { + den = detval2 - 0.5*log(DETX) - nmk*log(EPE/(2*sige)) + } else { + den = detval2 - 0.5*log(DETX) - nmk*log(EPE) + } + adj = max(den) + den = den - adj + den = exp(den) + nden = length(den) + + isum = sum((detval1[2:nden] + detval1[1:(nden-1)]) * (den[2:nden] - + den[1:(nden-1)])/2) + z = abs(den/isum) + den = cumsum(z) + rnd = runif(1) * sum(z) + ind = which(den <= rnd) + idraw = max(ind) + if((idraw > 0) && (idraw < nrho)) { + lambda = detval1[idraw] + } + } + psave[iter] = as.vector(lambda) + } +### % end of sampling loop + timings[["sampling"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + mat <- cbind(bsave, psave, ssave) + colnames(mat) <- c(colnames(x), "lambda", "sige") + res <- coda::mcmc(mat, start=con$nomit+1, end=con$ndraw, thin=con$thin) + + emixedImps <- NULL + if (etype == "emixed") { + sum_stats <- summary(res)$statistics + if (isTRUE(Durbin)) { + odd <- (m%/%2) > 0 + if (odd) { + m2 <- (m-1)/2 + } else { + m2 <- m/2 + } + if (K == 1 && odd) { + warning("model configuration issue: no total impacts") + } else { + cm <- matrix(0, ncol=m, nrow=m2) + if (K == 2) { + if (odd) { + rownames(cm) <- xxcolnames[2:(m2+1)] + } else { + rownames(cm) <- xxcolnames[1:m2] + } + for (i in 1:m2) cm[i, c(i+1, i+(m2+1))] <- 1 +# drop bug fix 2016-09-21 Philipp Hunziker + dirImps <- sum_stats[2:(m2+1), 1:2, drop=FALSE] + rownames(dirImps) <- rownames(cm) + indirImps <- sum_stats[(m2+2):m, 1:2, drop=FALSE] + rownames(indirImps) <- rownames(cm) + tI <- res[, 2:(m2+1), drop=FALSE] + + res[, (m2+2):m, drop=FALSE] + totImps <- t(apply(tI, 2, function(x) c(mean(x), sd(x)))) + rownames(totImps) <- rownames(cm) + } else { + rownames(cm) <- xxcolnames[1:m2] + for (i in 1:m2) cm[i, c(i, i+m2)] <- 1 + dirImps <- sum_stats[1:m2, 1:2, drop=FALSE] + rownames(dirImps) <- rownames(cm) + indirImps <- sum_stats[(m2+1):m, 1:2, drop=FALSE] + rownames(indirImps) <- rownames(cm) + tI <- res[, 1:m2, drop=FALSE] + + res[, (m2+1):m, drop=FALSE] + totImps <- t(apply(tI, 2, function(x) c(mean(x), sd(x)))) + rownames(totImps) <- rownames(cm) + colnames(totImps) <- colnames(dirImps) + } + } + } else if (is.formula(Durbin)) { +#FIXME + m <- sum(dvars) + m2 <- dvars[2] + cm <- matrix(0, ncol=m, nrow=m2) + for (i in 1:m2) { + cm[i, c(inds[i], i+dvars[1])] <- 1 + } + rownames(cm) <- wxn + dirImps <- sum_stats[2:dvars[1], 1:2, drop=FALSE] + rownames(dirImps) <- xn + indirImps <- sum_stats[(dvars[1]+1):m, 1:2, drop=FALSE] + if (!is.null(zero_fill)) { + if (length(zero_fill) > 0L) { + lres <- vector(mode="list", length=2L) + for (j in 1:2) { + jindirImps <- rep(as.numeric(NA), (dvars[1]-1)) + for (i in seq(along=inds)) { + jindirImps[(inds[i]-1)] <- indirImps[i, j] + } + lres[[j]] <- jindirImps + } + indirImps <- do.call("cbind", lres) + } + } + rownames(indirImps) <- xn + tI <- res[, inds, drop=FALSE] + + res[, (dvars[1]+1):m, drop=FALSE] + totImps <- t(apply(tI, 2, function(x) c(mean(x), sd(x)))) + if (!is.null(zero_fill)) { + if (length(zero_fill) > 0L) { + lres <- vector(mode="list", length=2L) + for (j in 1:2) { + jtotImps <- dirImps[, j] + for (i in seq(along=inds)) { + jtotImps[(inds[i]-1)] <- totImps[i, j] + } + lres[[j]] <- jtotImps + } + totImps <- do.call("cbind", lres) + } + } + rownames(totImps) <- xn + colnames(totImps) <- colnames(dirImps) + } else stop("undefined emixed state") + emixedImps <- list(dirImps=dirImps, indirImps=indirImps, + totImps=totImps) + } + + + timings[["finalise"]] <- proc.time() - .ptime_start + attr(res, "timings") <- do.call("rbind", timings) +# attr(res, "nano") <- c(nano_1, nano_2, nano_3) + attr(res, "control") <- con + attr(res, "etype") <- etype + attr(res, "listw_style") <- listw$style +# attr(res, "ll_mean") <- as.vector(ll_mean) + attr(res, "aliased") <- aliased + attr(res, "dvars") <- dvars + attr(res, "emixedImps") <- emixedImps + attr(res, "acc_rate") <- acc_rate + attr(res, "cc") <- cc + attr(res, "n") <- n + attr(res, "k") <- k + attr(res, "MH") <- priors$lambdaMH + class(res) <- c("MCMC_sem_g", class(res)) + res + +#output mcmc object + +} + +impacts.MCMC_sem_g <- function(obj, ..., tr=NULL, listw=NULL, evalues=NULL, + Q=NULL) { + emixedImps <- attr(obj, "emixedImps") + if (is.null(emixedImps)) { + stop("No indirect impacts, use summary()") + } + n <- attr(obj, "n") + k <- attr(obj, "k") + impactsWX(emixedImps, n, k, type="SDEM", method="MCMC") +} + + +# translated from Matlab code sac_g.m in the Spatial Econometrics toolbox by +# James LeSage and R. Kelley Pace (http://www.spatial-econometrics.com/). +# GSoc 2011 project by Abhirup Mallik mentored by Virgilio Gómez-Rubio + +spBreg_sac <- function(formula, data = list(), listw, listw2=NULL, na.action, + Durbin, type, zero.policy=NULL, control=list()) { + timings <- list() + .ptime_start <- proc.time() +#control + con <- list(tol.opt=.Machine$double.eps^0.5, ldet_method="SE_classic", + Imult=2, cheb_q=5, MC_p=16L, MC_m=30L, super=NULL, spamPivot="MMD", + in_coef=0.1, type="MC", correct=TRUE, trunc=TRUE, + SE_method="LU", nrho=200, interpn=2000, SElndet=NULL, LU_order=FALSE, + pre_eig1=NULL, pre_eig2=NULL, interval1=c(-1, 1), interval2=c(-1, 1), + ndraw=2500L, nomit=500L, thin=1L, verbose=FALSE, detval1=NULL, + detval2=NULL, prior=list(Tbeta=NULL, c_beta=NULL, lambda=0.5, + rho=0.5, sige=1, nu=0, d0=0, a1 = 1.01, a2 = 1.01, cc1 = 0.2, + cc2 = 0.2)) + priors <- con$prior + nmsP <- names(priors) + priors[(namp <- names(control$prior))] <- control$prior + if (length(noNms <- namp[!namp %in% nmsP])) + warning("unknown names in control$prior: ", + paste(noNms, collapse = ", ")) + control$prior <- NULL + con$prior <- NULL + nmsC <- names(con) + con[(namc <- names(control))] <- control + if (length(noNms <- namc[!namc %in% nmsC])) + warning("unknown names in control: ", paste(noNms, collapse = ", ")) + stopifnot(is.logical(con$verbose)) + stopifnot(is.integer(con$ndraw)) + stopifnot(is.integer(con$nomit)) + stopifnot(is.integer(con$thin)) + + if (is.null(zero.policy)) + zero.policy <- get.ZeroPolicyOption() + stopifnot(is.logical(zero.policy)) + if (class(formula) != "formula") formula <- as.formula(formula) + mt <- terms(formula, data = data) + mf <- lm(formula, data, na.action=na.action, method="model.frame") + na.act <- attr(mf, "na.action") + if (!inherits(listw, "listw")) stop("No neighbourhood list") + if (is.null(listw2)) listw2 <- listw + if (!is.null(con$pre_eig1) && is.null(con$pre_eig2)) + con$pre_eig2 <- con$pre_eig1 + else if (!inherits(listw2, "listw")) stop("No 2nd neighbourhood list") + can.sim <- FALSE + if (listw$style %in% c("W", "S")) can.sim <- can.be.simmed(listw) + can.sim2 <- FALSE + if (listw2$style %in% c("W", "S")) can.sim2 <- can.be.simmed(listw2) + if (!is.null(na.act)) { + subset <- !(1:length(listw$neighbours) %in% na.act) + listw <- subset(listw, subset, zero.policy=zero.policy) + } + y <- model.extract(mf, "response") + x <- model.matrix(mt, mf) + n <- nrow(x) + if (n != length(listw$neighbours)) + stop("Input data and weights have different dimensions") + xcolnames <- colnames(x) + K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) + wy <- lag.listw(listw, y, zero.policy=zero.policy) + if (anyNA(wy)) stop("NAs in lagged dependent variable") +#create_WX +# check for dgCMatrix + if (missing(type)) type <- "sac" + if (type == "sacmixed") { + type <- "Durbin" + } + if (missing(Durbin)) Durbin <- ifelse(type == "sac", FALSE, TRUE) + if (listw$style != "W" && is.formula(Durbin)) { + Durbin <- TRUE + warning("formula Durbin requires row-standardised weights; set TRUE") + } + if (is.logical(Durbin) && isTRUE(Durbin)) type <- "Durbin" + if (is.formula(Durbin)) type <- "Durbin" + if (is.logical(Durbin) && !isTRUE(Durbin)) type <- "sac" + + m <- ncol(x) + dvars <- c(NCOL(x), 0L) + +# if (type == "Durbin") { +# WX <- create_WX(x, listw, zero.policy=zero.policy, prefix="lag") +#FIXME + if (is.formula(Durbin) || isTRUE(Durbin)) { + prefix <- "lag" + if (isTRUE(Durbin)) { + WX <- create_WX(x, listw, zero.policy=zero.policy, + prefix=prefix) + } else { + dmf <- lm(Durbin, data, na.action=na.action, + method="model.frame") + fx <- try(model.matrix(Durbin, dmf), silent=TRUE) + if (class(fx) == "try-error") + stop("Durbin variable mis-match") + WX <- create_WX(fx, listw, zero.policy=zero.policy, + prefix=prefix) + inds <- match(substring(colnames(WX), 5, + nchar(colnames(WX))), colnames(x)) + if (anyNA(inds)) stop("WX variables not in X: ", + paste(substring(colnames(WX), 5, + nchar(colnames(WX)))[is.na(inds)], collapse=" ")) + icept <- grep("(Intercept)", colnames(x)) + iicept <- length(icept) > 0L + if (iicept) { + xn <- colnames(x)[-1] + } else { + xn <- colnames(x) + } + wxn <- substring(colnames(WX), nchar(prefix)+2, + nchar(colnames(WX))) + zero_fill <- NULL + if (length((which(!(xn %in% wxn)))) > 0L) + zero_fill <- length(xn) + (which(!(xn %in% wxn))) + } + dvars <- c(NCOL(x), NCOL(WX)) + if (is.formula(Durbin)) { + attr(dvars, "f") <- Durbin + attr(dvars, "inds") <- inds + attr(dvars, "zero_fill") <- zero_fill + } + x <- cbind(x, WX) + m <- NCOL(x) + rm(WX) + } + if (NROW(x) != length(listw2$neighbours)) + stop("Input data and neighbourhood list2 have different dimensions") + w2y <- lag.listw(listw2, y, zero.policy=zero.policy) + w2wy <- lag.listw(listw2, wy, zero.policy=zero.policy) + lm.base <- lm(y ~ x - 1) # doesn't like dgCMatrix + aliased <- is.na(coefficients(lm.base)) + cn <- names(aliased) + names(aliased) <- substr(cn, 2, nchar(cn)) + if (any(aliased)) { + if (is.formula(Durbin)) { + stop("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + } else { + warning("Aliased variables found: ", + paste(names(aliased)[aliased], collapse=" ")) + nacoef <- which(aliased) + x <- x[,-nacoef] + } + } + m <- ncol(x) + xcolnames <- colnames(x) + K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) + if (any(is.na(wy))) + stop("NAs in lagged dependent variable") + if (m > 1 || (m == 1 && K == 1)) { + W2X <- matrix(nrow=n,ncol=(m-(K-1))) + for (k in K:m) { + wx <- lag.listw(listw2, x[,k], zero.policy=zero.policy) + if (any(is.na(wx))) + stop("NAs in lagged independent variable") + W2X[,(k-(K-1))] <- wx + } + } + if (K == 2) { + wx1 <- as.double(rep(1, n)) + wx <- lag.listw(listw2, wx1, zero.policy=zero.policy) + if (m > 1) W2X <- cbind(wx, W2X) + else W2X <- matrix(wx, nrow=n, ncol=1) + } + colnames(W2X) <- xcolnames + rm(wx) + + timings[["set_up"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + env <- new.env() + assign("can.sim", can.sim, envir=env) + assign("can.sim2", can.sim2, envir=env) + assign("listw", listw, envir=env) + assign("listw2", listw2, envir=env) + assign("similar", FALSE, envir=env) + assign("similar2", FALSE, envir=env) + assign("n", n, envir=env) + assign("verbose", con$verbose, envir=env) + assign("family", "SAR", envir=env) + assign("method", con$ldet_method, envir=env) + W <- as(listw, "CsparseMatrix") + assign("W", W, envir=env) + + con$interval1 <- jacobianSetup(con$ldet_method, env, con, + pre_eig=con$pre_eig1, interval=con$interval1, which=1) + detval11 <- get("detval1", envir=env)[,1] + detval12 <- get("detval1", envir=env)[,2] + con$interval2 <- jacobianSetup(con$ldet_method, env, con, + pre_eig=con$pre_eig2, interval=con$interval2, which=2) + detval21 <- get("detval2", envir=env)[,1] + detval22 <- get("detval2", envir=env)[,2] + + nm <- paste(con$ldet_method, "set_up", sep="_") + timings[[nm]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + k <- m + + if (is.null(priors$c_beta)) + priors$c_beta <- rep(0, k) + else + stopifnot(length(priors$c_beta) == k) + + if (is.null(priors$Tbeta)) + priors$Tbeta <- diag(k)*1e+12 + else + stopifnot(nrow(priors$Tbeta) == k && ncol(priors$Tbeta) == k) + + sige <- priors$sige + rho <- priors$rho + lambda <- priors$lambda + +#% storage for draws + bsave <- matrix(0, nrow=con$ndraw, ncol=k) + psave <- numeric(con$ndraw) + lsave <- numeric(con$ndraw) + ssave <- numeric(con$ndraw) + acc_rate1 <- numeric(con$ndraw) + acc_rate2 <- numeric(con$ndraw) + +#% ====== initializations +#% compute this stuff once to save time + + TI = solve(priors$Tbeta); # see eq 5.29, Lesage & Pace (2009) p. 140 + TIc = TI%*%priors$c_beta; + + xpx = crossprod(x) + xpy = crossprod(x, y) + xpWy = crossprod(x, wy) + nu1 = n + 2*priors$nu + nrho = length(detval11) + gsize1 = detval11[2] - detval11[1] + gsize2 = detval21[2] - detval21[1] + acc1 = 0 + acc2 = 0 + nmk = (n-k)/2 + cc1 <- priors$cc1 + cc2 <- priors$cc2 +# nano_1 = 0 +# nano_2 = 0 +# nano_3 = 0 + + timings[["complete_setup"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + + for (iter in 1:con$ndraw) { #% start sampling; + + ##% update beta + xss = x - lambda*W2X + AI = solve(crossprod(xss) + sige*TI) + yss = y - rho*wy -lambda*w2y + rho*lambda*w2wy + b = crossprod(xss, yss) + sige*TIc + b0 = AI %*% b + bhat = MASS::mvrnorm(1, b0, sige*AI) + bsave[iter, 1:k] = as.vector(bhat) + + ##% update sige +# FIXME - sources do not use bhat + lbhat <- solve(crossprod(xss)) %*% crossprod(xss, yss) + e = yss - xss %*% lbhat + ed = e - lambda*lag.listw(listw2, e, zero.policy=zero.policy) + d1 = 2*priors$d0 + crossprod(ed) + chi = rchisq(1, nu1) + sige = as.numeric(d1/chi) + ssave[iter] = as.vector(sige) + + #update lambda using M-H + i1 = max(which(detval21 <= (lambda + gsize2))) + i2 = max(which(detval21 <= (lambda - gsize2))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 #Fixed this + detm = detval22[index] + epe = (crossprod(e))/(2*sige) + rhox = detm - epe + accept = 0L; + lambda2 = lambda + cc2*rnorm(1); + while(accept == 0L) { + if((lambda2 > con$interval2[1]) & (lambda2 < con$interval2[2])) { + accept=1 + } else { + lambda2 = lambda + cc2*rnorm(1) + } + } + i1 = max(which(detval21 <= (lambda2 + gsize2))) + i2 = max(which(detval21 <= (lambda2 - gsize2))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 #Fixed this + detm = detval22[index] + xss = x - lambda2*W2X + yss = y - rho*wy - lambda2*w2y + rho*lambda2*w2wy + lbhat <- solve(crossprod(xss)) %*% crossprod(xss, yss) + e = yss - xss %*% lbhat + epe = (crossprod(e))/(2*sige) + rhoy = detm - epe + if ((rhoy - rhox) > exp(1)) { + p = 1 + } else { + ratio = exp(rhoy-rhox) + p = min(1, ratio) + } + ru = runif(1) + if(ru < p) { + lambda = lambda2 + acc2 = acc2 + 1 + } + acc_rate2[iter] = acc2/iter + if(acc_rate2[iter] < 0.4) cc2=cc2/1.1 + if(acc_rate2[iter] > 0.6) cc2=cc2*1.1 + lsave[iter] = as.vector(lambda) + +## % metropolis step to get rho update + i1 = max(which(detval11 <= (rho + gsize1))) + i2 = max(which(detval11 <= (rho - gsize1))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 + detm = detval12[index] + xss = x - lambda*W2X + yss = y - rho*wy - lambda*w2y + rho*lambda*w2wy + lbhat <- solve(crossprod(xss)) %*% crossprod(xss, yss) + e = yss - xss %*% lbhat + epe = (crossprod(e))/(2*sige) + rhox = detm - epe + accept = 0L + rho2 = rho + cc1*rnorm(1); + while(accept == 0L) { + if((rho2 > con$interval1[1]) & (rho2 < con$interval1[2])) { + accept=1 + } else { + rho2 = rho + cc1*rnorm(1) + } + } + i1 = max(which(detval11 <= (rho2 + gsize1))) + i2 = max(which(detval11 <= (rho2 - gsize1))) + index = round((i1+i2)/2) + if (!is.finite(index)) index = 1 + detm = detval12[index] + yss = y - rho2*wy - lambda*w2y + rho2*lambda*w2wy + lbhat <- solve(crossprod(xss)) %*% crossprod(xss, yss) + e = yss - xss %*% lbhat + epe = (crossprod(e))/(2*sige) + rhoy = detm - epe + if ((rhoy - rhox) > exp(1)) { + p = 1 + } else { + ratio = exp(rhoy-rhox) + p = min(1, ratio) + } + ru = runif(1) + if(ru < p) { + rho = rho2 + acc1 = acc1 + 1 + } + acc_rate1[iter] = acc1/iter + if(acc_rate1[iter] < 0.4) cc1=cc1/1.1 + if(acc_rate1[iter] > 0.6) cc1=cc1*1.1 + + psave[iter] = as.vector(rho) + + } +### % end of sampling loop + timings[["sampling"]] <- proc.time() - .ptime_start + .ptime_start <- proc.time() + + mat <- cbind(bsave, psave, lsave, ssave) + colnames(mat) <- c(colnames(x), "rho", "lambda", "sige") + res <- coda::mcmc(mat, start=con$nomit+1, end=con$ndraw, thin=con$thin) + means <- summary(res)$statistics[,1] + beta <- means[1:(length(means)-3)] + rho <- means[(length(means)-2)] + lambda <- means[(length(means)-1)] + xb = (x - lambda*W2X) %*% beta + ys = y - rho*wy - lambda*w2y + rho*lambda*w2wy + e = (ys - xb) + sse <- crossprod(e) + s2 <- sse/n + ldet <- do_ldet(rho, env) + ll_mean <- (ldet - ((n/2)*log(2*pi)) - (n/2)*log(s2) + - (1/(2*s2))*sse) + + timings[["finalise"]] <- proc.time() - .ptime_start + attr(res, "timings") <- do.call("rbind", timings) + attr(res, "control") <- con + attr(res, "type") <- type + attr(res, "listw_style") <- listw$style + attr(res, "ll_mean") <- as.vector(ll_mean) + attr(res, "aliased") <- aliased + attr(res, "acc_rate1") <- acc_rate1 + attr(res, "acc_rate2") <- acc_rate2 + attr(res, "dvars") <- dvars + class(res) <- c("MCMC_sac_g", class(res)) + res + +#output mcmc object +} + + +impacts.MCMC_sac_g <- function(obj, ..., tr=NULL, listw=NULL, evalues=NULL, + Q=NULL) { + obj_lag <- obj[, -(which(colnames(obj) == "lambda"))] + attributes(obj_lag) <- c(attributes(obj_lag), + attributes(obj)[5:13]) + if (attr(obj, "type") == "sac") { + attr(obj_lag, "type") <- "lag" + } else { + attr(obj_lag, "type") <- "Durbin" + } + class(obj_lag) <- c("MCMC_sar_g", class(obj_lag)) + res <- impacts.MCMC_sar_g(obj_lag, tr=tr, listw=listw, evalues=evalues, + Q=Q) + if (attr(obj, "type") == "sac") { + attr(res, "type") <- "sac" + } else { + attr(res, "type") <- "sacmixed" + } + res +} diff -Nru r-cran-spdep-0.7-9+dfsg/R/utils.R r-cran-spdep-0.8-1+dfsg/R/utils.R --- r-cran-spdep-0.7-9+dfsg/R/utils.R 2017-11-01 16:24:24.000000000 +0000 +++ r-cran-spdep-0.8-1+dfsg/R/utils.R 2018-10-12 09:17:43.000000000 +0000 @@ -79,6 +79,8 @@ if (!inherits(listw, "listw")) stop(paste(deparse(substitute(x)), "is not a listw object")) x <- var +# https://github.com/r-spatial/spdep/issues/19 + if (!is.vector(x) && !is.matrix(x)) x <- c(x) if (!is.vector(c(x)) && !is.matrix(x)) stop(paste(deparse(substitute(var)), "not a vector or matrix")) if (!is.numeric(x)) stop(paste(deparse(substitute(var)), @@ -93,7 +95,9 @@ as.logical(zero.policy), naok=NAOK, PACKAGE="spdep") } else { if (nrow(x) != n) stop("object lengths differ") - res <- matrix(0, nrow=nrow(x), ncol=ncol(x)) +# https://github.com/r-spatial/spdep/issues/19 + res <- matrix(0, nrow=nrow(x), + ncol=ifelse(is.na(ncol(x)), 1, ncol(x))) for (i in 1:ncol(x)) { res[,i] <- .Call("lagw", listw$neighbours, listw$weights, x[,i],