MCMC sample from fitted spatial regression
The MCMCsamp
method uses rwmetrop
, a random walk Metropolis algorithm, from LearnBayes to make MCMC samples from fitted maximum likelihood spatial regression models.
MCMCsamp(object, mcmc = 1L, verbose = NULL, ...) ## S3 method for class 'Spautolm' MCMCsamp(object, mcmc = 1L, verbose = NULL, ..., burnin = 0L, scale=1, listw, control = list()) ## S3 method for class 'Sarlm' MCMCsamp(object, mcmc = 1L, verbose = NULL, ..., burnin=0L, scale=1, listw, listw2=NULL, control=list())
object |
A spatial regression model object fitted by maximum likelihood with |
mcmc |
The number of MCMC iterations after burnin |
verbose |
default NULL, use global option value; if TRUE, reports progress |
... |
Arguments passed through |
burnin |
The number of burn-in iterations for the sampler |
scale |
a positive scale parameter |
listw, listw2 |
|
control |
list of extra control arguments - see |
An object of class “mcmc” suited to coda, with attributes: “accept” acceptance rate; “type” input ML fitted model type “SAR”, “CAR”, “SMA”, “lag”, “mixed”, “error”, “sac”, “sacmixed”; “timings” run times
If the acceptance rate is below 0.05, a warning will be issued; consider increasing mcmc.
Roger Bivand Roger.Bivand@nhh.no
Jim Albert (2007) Bayesian Computation with R, Springer, New York, pp. 104-105.
require("sf", quietly=TRUE) nydata <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData")[1], quiet=TRUE) suppressMessages(nyadjmat <- as.matrix(foreign::read.dbf(system.file( "misc/nyadjwts.dbf", package="spData")[1])[-1])) suppressMessages(ID <- as.character(names(foreign::read.dbf(system.file( "misc/nyadjwts.dbf", package="spData")[1]))[-1])) identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) #require("spdep", quietly=TRUE) nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") summary(esar1f) res <- MCMCsamp(esar1f, mcmc=1000, burnin=200, listw=listw_NY) summary(res) ## Not run: esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="eigen") summary(esar1fw) res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ecar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="CAR", method="eigen") summary(ecar1f) res <- MCMCsamp(ecar1f, mcmc=5000, burnin=500, listw=listw_NY) summary(res) esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="eigen") summary(esar1fw) res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ecar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="CAR", method="eigen") summary(ecar1fw) res <- MCMCsamp(ecar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ## End(Not run) esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(esar0) res <- MCMCsamp(esar0, mcmc=1000, burnin=200, listw=listw_NY) summary(res) ## Not run: esar0w <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8) summary(esar0) res <- MCMCsamp(esar0w, mcmc=5000, burnin=500, listw=listw_NY) summary(res) esar1 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, etype="emixed") summary(esar1) res <- MCMCsamp(esar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) lsar0 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(lsar0) res <- MCMCsamp(lsar0, mcmc=5000, burnin=500, listw=listw_NY) summary(res) lsar1 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, type="mixed") summary(lsar1) res <- MCMCsamp(lsar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ssar0 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(ssar0) res <- MCMCsamp(ssar0, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ssar1 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, type="sacmixed") summary(ssar1) res <- MCMCsamp(ssar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) ## End(Not run)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.