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 Spautolm
MCMCsamp(object, mcmc = 1L, verbose = NULL, ...,
 burnin = 0L, scale=1, listw, control = list())
# S3 method for Sarlm
MCMCsamp(object, mcmc = 1L, verbose = NULL, ...,
    burnin=0L, scale=1, listw, listw2=NULL, control=list())

Arguments

object

A spatial regression model object fitted by maximum likelihood with spautolm

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

listw objects created for example by nb2listw; should be the same object(s) used for fitting the model

control

list of extra control arguments - see spautolm

Value

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

Note

If the acceptance rate is below 0.05, a warning will be issued; consider increasing mcmc.

References

Jim Albert (2007) Bayesian Computation with R, Springer, New York, pp. 104-105.

Author

Roger Bivand Roger.Bivand@nhh.no

See also

Examples

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))
#> [1] TRUE
#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)
#> 
#> Call: 
#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, 
#>     listw = listw_NY, family = "SAR", method = "eigen")
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.56754 -0.38239 -0.02643  0.33109  4.01219 
#> 
#> Coefficients: 
#>              Estimate Std. Error z value  Pr(>|z|)
#> (Intercept) -0.618193   0.176784 -3.4969 0.0004707
#> PEXPOSURE    0.071014   0.042051  1.6888 0.0912635
#> PCTAGE65P    3.754200   0.624722  6.0094 1.862e-09
#> PCTOWNHOME  -0.419890   0.191329 -2.1946 0.0281930
#> 
#> Lambda: 0.040487 LR test value: 5.2438 p-value: 0.022026 
#> Numerical Hessian standard error of lambda: 0.017209 
#> 
#> Log likelihood: -276.1069 
#> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333)
#> Number of observations: 281 
#> Number of parameters estimated: 6 
#> AIC: 564.21
#> 
res <- MCMCsamp(esar1f, mcmc=1000, burnin=200, listw=listw_NY)
summary(res)
#> 
#> Iterations = 1:1000
#> Thinning interval = 1 
#> Number of chains = 1 
#> Sample size per chain = 1000 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>                 Mean      SD  Naive SE Time-series SE
#> lambda       0.04312 0.01738 0.0005496       0.002193
#> (Intercept) -0.64883 0.18254 0.0057724       0.023256
#> PEXPOSURE    0.08483 0.04913 0.0015537       0.006320
#> PCTAGE65P    3.78232 0.58084 0.0183678       0.066337
#> PCTOWNHOME  -0.39612 0.22118 0.0069943       0.030304
#> 
#> 2. Quantiles for each variable:
#> 
#>                  2.5%      25%      50%      75%     97.5%
#> lambda       0.008390  0.03207  0.04215  0.05501  0.076133
#> (Intercept) -1.039676 -0.77587 -0.64419 -0.51195 -0.326557
#> PEXPOSURE   -0.004963  0.05359  0.08142  0.11521  0.202348
#> PCTAGE65P    2.627258  3.39980  3.73455  4.18334  4.962310
#> PCTOWNHOME  -0.819651 -0.56246 -0.39274 -0.22634  0.004321
#> 
if (FALSE) {
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)
}
esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
 listw=listw_NY)
summary(esar0)
#> 
#> Call:errorsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, 
#>     data = nydata, listw = listw_NY)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.56754 -0.38239 -0.02643  0.33109  4.01219 
#> 
#> Type: error 
#> Coefficients: (asymptotic standard errors) 
#>              Estimate Std. Error z value  Pr(>|z|)
#> (Intercept) -0.618193   0.176784 -3.4969 0.0004707
#> PEXPOSURE    0.071014   0.042051  1.6888 0.0912635
#> PCTAGE65P    3.754200   0.624722  6.0094 1.862e-09
#> PCTOWNHOME  -0.419890   0.191329 -2.1946 0.0281930
#> 
#> Lambda: 0.040487, LR test value: 5.2438, p-value: 0.022026
#> Asymptotic standard error: 0.016214
#>     z-value: 2.4971, p-value: 0.01252
#> Wald statistic: 6.2356, p-value: 0.01252
#> 
#> Log likelihood: -276.1069 for error model
#> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333)
#> Number of observations: 281 
#> Number of parameters estimated: 6 
#> AIC: 564.21, (AIC for lm: 567.46)
#> 
res <- MCMCsamp(esar0, mcmc=1000, burnin=200, listw=listw_NY)
summary(res)
#> 
#> Iterations = 1:1000
#> Thinning interval = 1 
#> Number of chains = 1 
#> Sample size per chain = 1000 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>                 Mean      SD  Naive SE Time-series SE
#> lambda       0.04086 0.01436 0.0004541       0.001639
#> (Intercept) -0.62555 0.18275 0.0057791       0.024443
#> PEXPOSURE    0.06900 0.04034 0.0012758       0.004965
#> PCTAGE65P    3.68791 0.64881 0.0205172       0.077680
#> PCTOWNHOME  -0.38263 0.19769 0.0062513       0.025188
#> 
#> 2. Quantiles for each variable:
#> 
#>                  2.5%      25%      50%      75%     97.5%
#> lambda       0.014717  0.03100  0.04103  0.05272  0.065715
#> (Intercept) -0.972641 -0.75486 -0.63637 -0.50400 -0.198074
#> PEXPOSURE   -0.004903  0.03891  0.06820  0.09560  0.152374
#> PCTAGE65P    2.382993  3.31371  3.67406  4.12972  4.919868
#> PCTOWNHOME  -0.779951 -0.49064 -0.38080 -0.25234 -0.003732
#> 
if (FALSE) {
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)
}