MCMCsamp.RdThe MCMCsamp method uses rwmetrop, a random walk Metropolis algorithm, from LearnBayes to make MCMC samples from fitted maximum likelihood spatial regression models.
A spatial regression model object fitted by maximum likelihood with spautolm
The number of MCMC iterations after burnin
default NULL, use global option value; if TRUE, reports progress
Arguments passed through
The number of burn-in iterations for the sampler
a positive scale parameter
listw objects created for example by nb2listw; should be the same object(s) used for fitting the model
list of extra control arguments - see spautolm
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.
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))
#> [1] TRUE
#require("spdep", quietly=TRUE)
nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY))
#> Warning: style is M (missing); style should be set to a valid value
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.017201
#>
#> 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.04611 0.01498 0.0004736 0.001756
#> (Intercept) -0.65331 0.19263 0.0060916 0.026540
#> PEXPOSURE 0.08968 0.04897 0.0015486 0.006607
#> PCTAGE65P 3.69656 0.52425 0.0165784 0.057839
#> PCTOWNHOME -0.39485 0.23949 0.0075732 0.037397
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda 0.015004 0.03557 0.04640 0.05531 0.07351
#> (Intercept) -1.089847 -0.77125 -0.63878 -0.51969 -0.34027
#> PEXPOSURE -0.001914 0.05288 0.09113 0.12367 0.18519
#> PCTAGE65P 2.672002 3.33982 3.67753 4.04290 4.74392
#> PCTOWNHOME -0.845780 -0.56384 -0.41154 -0.21860 0.07637
#>
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)
}