MCMC sample from fitted spatial regression
MCMCsamp.RdThe MCMCsamp method uses rwmetrop, a random walk Metropolis algorithm, from LearnBayes to make MCMC samples from fitted maximum likelihood spatial regression models.
Usage
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())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
listwobjects created for example bynb2listw; 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
Author
Roger Bivand Roger.Bivand@nhh.no
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)
listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), 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.017197
#>
#> 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.04424 0.01759 0.0005561 0.002587
#> (Intercept) -0.60387 0.17127 0.0054160 0.019939
#> PEXPOSURE 0.07818 0.03895 0.0012317 0.004971
#> PCTAGE65P 3.61380 0.61721 0.0195179 0.078471
#> PCTOWNHOME -0.42547 0.19215 0.0060763 0.022985
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda 0.010601 0.03057 0.04463 0.05636 0.07645
#> (Intercept) -0.909163 -0.73807 -0.58875 -0.48939 -0.25674
#> PEXPOSURE 0.003061 0.04789 0.07901 0.10917 0.15954
#> PCTAGE65P 2.376203 3.26135 3.65765 3.97016 4.82490
#> PCTOWNHOME -0.813304 -0.55636 -0.43884 -0.30986 -0.08202
#>
# \dontrun{
esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="SAR", method="eigen")
summary(esar1fw)
#>
#> Call:
#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, weights = POP8, family = "SAR", method = "eigen")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.48488 -0.26823 0.09489 0.46552 4.28343
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08
#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473
#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11
#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975
#>
#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764
#> Numerical Hessian standard error of lambda: 0.016356
#>
#> Log likelihood: -251.6017
#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
#> Number of observations: 281
#> Number of parameters estimated: 6
#> AIC: NA (not available for weighted model)
#>
res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.01332 0.01713 0.0002422 0.001040
#> (Intercept) -0.81092 0.14834 0.0020979 0.008623
#> PEXPOSURE 0.08370 0.03067 0.0004338 0.001764
#> PCTAGE65P 3.78962 0.58047 0.0082091 0.031503
#> PCTOWNHOME -0.37144 0.16490 0.0023320 0.009809
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda -0.02191 0.001749 0.01315 0.02592 0.04579
#> (Intercept) -1.09936 -0.906434 -0.80683 -0.70908 -0.52586
#> PEXPOSURE 0.02775 0.063093 0.08289 0.10397 0.14414
#> PCTAGE65P 2.61660 3.431381 3.79540 4.18104 4.91572
#> PCTOWNHOME -0.69168 -0.487400 -0.38088 -0.25202 -0.05274
#>
ecar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, family="CAR", method="eigen")
summary(ecar1f)
#>
#> Call:
#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, family = "CAR", method = "eigen")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.539732 -0.384311 -0.030646 0.335126 3.808848
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.648362 0.181129 -3.5796 0.0003442
#> PEXPOSURE 0.077899 0.043692 1.7829 0.0745986
#> PCTAGE65P 3.703830 0.627185 5.9055 3.516e-09
#> PCTOWNHOME -0.382789 0.195564 -1.9574 0.0503053
#>
#> Lambda: 0.084123 LR test value: 5.8009 p-value: 0.016018
#> Numerical Hessian standard error of lambda: 0.030853
#>
#> Log likelihood: -275.8283
#> ML residual variance (sigma squared): 0.40758, (sigma: 0.63842)
#> Number of observations: 281
#> Number of parameters estimated: 6
#> AIC: 563.66
#>
res <- MCMCsamp(ecar1f, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.08683 0.02838 0.0004014 0.001541
#> (Intercept) -0.67584 0.21215 0.0030003 0.012677
#> PEXPOSURE 0.08942 0.05299 0.0007493 0.003300
#> PCTAGE65P 3.64558 0.65339 0.0092403 0.038815
#> PCTOWNHOME -0.36821 0.22630 0.0032004 0.013385
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda 0.024844 0.06884 0.08851 0.1078 0.13505
#> (Intercept) -1.108206 -0.81500 -0.66971 -0.5255 -0.27879
#> PEXPOSURE -0.009761 0.05174 0.08631 0.1236 0.20173
#> PCTAGE65P 2.378676 3.19803 3.63371 4.1368 4.85338
#> PCTOWNHOME -0.808235 -0.52404 -0.36883 -0.2197 0.07086
#>
esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="SAR", method="eigen")
summary(esar1fw)
#>
#> Call:
#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, weights = POP8, family = "SAR", method = "eigen")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.48488 -0.26823 0.09489 0.46552 4.28343
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08
#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473
#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11
#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975
#>
#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764
#> Numerical Hessian standard error of lambda: 0.016356
#>
#> Log likelihood: -251.6017
#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
#> Number of observations: 281
#> Number of parameters estimated: 6
#> AIC: NA (not available for weighted model)
#>
res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.01324 0.01634 0.0002311 0.0009562
#> (Intercept) -0.79796 0.15741 0.0022261 0.0097176
#> PEXPOSURE 0.08088 0.03206 0.0004534 0.0018172
#> PCTAGE65P 3.78034 0.58893 0.0083288 0.0345906
#> PCTOWNHOME -0.37812 0.16071 0.0022727 0.0090656
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda -0.01751 0.001606 0.01336 0.02384 0.04498
#> (Intercept) -1.11299 -0.907786 -0.79649 -0.69600 -0.49941
#> PEXPOSURE 0.02125 0.059565 0.08099 0.10241 0.14586
#> PCTAGE65P 2.68306 3.369273 3.77554 4.21312 4.93407
#> PCTOWNHOME -0.70673 -0.481221 -0.37593 -0.27548 -0.06440
#>
ecar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="CAR", method="eigen")
summary(ecar1fw)
#>
#> Call:
#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, weights = POP8, family = "CAR", method = "eigen")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.491042 -0.270906 0.081435 0.451556 4.198134
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.790154 0.144862 -5.4545 4.910e-08
#> PEXPOSURE 0.081922 0.028593 2.8651 0.004169
#> PCTAGE65P 3.825858 0.577720 6.6223 3.536e-11
#> PCTOWNHOME -0.386820 0.157436 -2.4570 0.014010
#>
#> Lambda: 0.022419 LR test value: 0.38785 p-value: 0.53343
#> Numerical Hessian standard error of lambda: 0.038916
#>
#> Log likelihood: -251.5711
#> ML residual variance (sigma squared): 1102.9, (sigma: 33.21)
#> Number of observations: 281
#> Number of parameters estimated: 6
#> AIC: NA (not available for weighted model)
#>
res <- MCMCsamp(ecar1fw, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.03676 0.04170 0.0005897 0.002549
#> (Intercept) -0.81674 0.16846 0.0023824 0.010592
#> PEXPOSURE 0.08920 0.03669 0.0005189 0.002637
#> PCTAGE65P 3.74362 0.58834 0.0083203 0.033752
#> PCTOWNHOME -0.37703 0.17603 0.0024894 0.010758
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda -0.04118 0.009495 0.03863 0.0640 0.11787
#> (Intercept) -1.18055 -0.920940 -0.81408 -0.7065 -0.48666
#> PEXPOSURE 0.02424 0.063144 0.08627 0.1108 0.16979
#> PCTAGE65P 2.54117 3.383146 3.76084 4.1392 4.87221
#> PCTOWNHOME -0.72276 -0.504797 -0.37302 -0.2542 -0.04689
#>
# }
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.04431 0.01630 0.0005154 0.002628
#> (Intercept) -0.63273 0.18769 0.0059354 0.030330
#> PEXPOSURE 0.08199 0.04317 0.0013653 0.006289
#> PCTAGE65P 3.49773 0.50939 0.0161083 0.053385
#> PCTOWNHOME -0.38817 0.20745 0.0065601 0.031330
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda 0.012461 0.03083 0.04728 0.05706 0.07158
#> (Intercept) -1.073224 -0.74209 -0.64456 -0.50287 -0.28163
#> PEXPOSURE 0.004178 0.05825 0.07866 0.10630 0.17035
#> PCTAGE65P 2.588327 3.17304 3.47441 3.86058 4.50569
#> PCTOWNHOME -0.761138 -0.53907 -0.40222 -0.24125 0.03827
#>
# \dontrun{
esar0w <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8)
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(esar0w, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.01328 0.01740 0.0002460 0.001022
#> (Intercept) -0.79295 0.14970 0.0021171 0.008980
#> PEXPOSURE 0.08111 0.02939 0.0004157 0.001769
#> PCTAGE65P 3.75967 0.60116 0.0085017 0.036245
#> PCTOWNHOME -0.37938 0.16255 0.0022988 0.009875
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda -0.02138 0.001516 0.01327 0.02514 0.04744
#> (Intercept) -1.09797 -0.894409 -0.78552 -0.68430 -0.52764
#> PEXPOSURE 0.02601 0.061233 0.08006 0.10111 0.13853
#> PCTAGE65P 2.59602 3.323984 3.74509 4.19572 4.87748
#> PCTOWNHOME -0.67869 -0.497202 -0.38603 -0.26800 -0.05950
#>
esar1 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, etype="emixed")
summary(esar1)
#>
#> Call:errorsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
#> data = nydata, listw = listw_NY, etype = "emixed")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.81562 -0.37641 -0.02224 0.33638 4.00054
#>
#> Type: error
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.118019 0.247425 -4.5186 6.225e-06
#> PEXPOSURE 0.218279 0.079245 2.7545 0.005879
#> PCTAGE65P 3.416477 0.645587 5.2920 1.210e-07
#> PCTOWNHOME 0.036593 0.249835 0.1465 0.883551
#> lag.(Intercept) 0.121515 0.057636 2.1083 0.035003
#> lag.PEXPOSURE -0.035075 0.015943 -2.2000 0.027808
#> lag.PCTAGE65P 0.263096 0.220118 1.1953 0.231989
#> lag.PCTOWNHOME -0.155680 0.059213 -2.6291 0.008560
#>
#> Lambda: 0.022723, LR test value: 1.6846, p-value: 0.19432
#> Asymptotic standard error: 0.017169
#> z-value: 1.3235, p-value: 0.18567
#> Wald statistic: 1.7516, p-value: 0.18567
#>
#> Log likelihood: -269.5398 for error model
#> ML residual variance (sigma squared): 0.39759, (sigma: 0.63055)
#> Number of observations: 281
#> Number of parameters estimated: 10
#> AIC: 559.08, (AIC for lm: 558.76)
#>
res <- MCMCsamp(esar1, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> lambda 0.02675 0.01582 0.0002237 0.001178
#> (Intercept) -1.09743 0.25377 0.0035888 0.020706
#> PEXPOSURE 0.21730 0.08427 0.0011917 0.006770
#> PCTAGE65P 3.41168 0.64492 0.0091206 0.049641
#> PCTOWNHOME -0.01402 0.22948 0.0032453 0.016651
#> lag.(Intercept) 0.12159 0.06097 0.0008622 0.005116
#> lag.PEXPOSURE -0.03595 0.01712 0.0002422 0.001468
#> lag.PCTAGE65P 0.24807 0.21767 0.0030783 0.016869
#> lag.PCTOWNHOME -0.14533 0.06036 0.0008536 0.004839
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> lambda -0.004004 0.01637 0.02646 0.03706 0.058414
#> (Intercept) -1.596596 -1.26665 -1.09711 -0.92951 -0.569075
#> PEXPOSURE 0.060047 0.15865 0.21404 0.27377 0.378642
#> PCTAGE65P 2.125474 2.99396 3.38851 3.80196 4.759975
#> PCTOWNHOME -0.432544 -0.18179 -0.02074 0.14040 0.469260
#> lag.(Intercept) 0.003745 0.08293 0.12176 0.16425 0.246727
#> lag.PEXPOSURE -0.069168 -0.04666 -0.03371 -0.02598 -0.003682
#> lag.PCTAGE65P -0.204498 0.09091 0.26104 0.40579 0.654141
#> lag.PCTOWNHOME -0.270422 -0.18572 -0.14366 -0.10596 -0.030349
#>
lsar0 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY)
summary(lsar0)
#>
#> Call:lagsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.586752 -0.391580 -0.022469 0.338017 4.029430
#>
#> Type: lag
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.514495 0.156154 -3.2948 0.000985
#> PEXPOSURE 0.047627 0.034509 1.3801 0.167542
#> PCTAGE65P 3.648198 0.599046 6.0900 1.129e-09
#> PCTOWNHOME -0.414601 0.169554 -2.4453 0.014475
#>
#> Rho: 0.038893, LR test value: 6.9683, p-value: 0.0082967
#> Asymptotic standard error: 0.015053
#> z-value: 2.5837, p-value: 0.0097755
#> Wald statistic: 6.6754, p-value: 0.0097755
#>
#> Log likelihood: -275.2447 for lag model
#> ML residual variance (sigma squared): 0.41166, (sigma: 0.6416)
#> Number of observations: 281
#> Number of parameters estimated: 6
#> AIC: 562.49, (AIC for lm: 567.46)
#> LM test for residual autocorrelation
#> test value: 1.4633, p-value: 0.22641
#>
res <- MCMCsamp(lsar0, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> rho 0.03870 0.01553 0.0002197 0.0009692
#> (Intercept) -0.52594 0.15180 0.0021467 0.0083453
#> PEXPOSURE 0.04892 0.03234 0.0004573 0.0017550
#> PCTAGE65P 3.69568 0.60023 0.0084885 0.0326213
#> PCTOWNHOME -0.40816 0.17766 0.0025125 0.0101486
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> rho 0.009508 0.02819 0.03772 0.04933 0.07033
#> (Intercept) -0.824584 -0.63648 -0.51958 -0.42098 -0.23010
#> PEXPOSURE -0.014658 0.02721 0.04915 0.07180 0.10892
#> PCTAGE65P 2.618603 3.25730 3.68615 4.11619 4.93251
#> PCTOWNHOME -0.748801 -0.52479 -0.40866 -0.29381 -0.01644
#>
lsar1 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, type="mixed")
summary(lsar1)
#>
#> Call:lagsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, type = "mixed")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.799308 -0.390125 -0.021371 0.346128 3.965251
#>
#> Type: mixed
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.131233 0.249631 -4.5316 5.853e-06
#> PEXPOSURE 0.218364 0.079301 2.7536 0.005894
#> PCTAGE65P 3.361158 0.654123 5.1384 2.771e-07
#> PCTOWNHOME 0.071903 0.253967 0.2831 0.777085
#> lag.(Intercept) 0.132544 0.056175 2.3595 0.018300
#> lag.PEXPOSURE -0.035239 0.015536 -2.2681 0.023322
#> lag.PCTAGE65P 0.161685 0.223690 0.7228 0.469798
#> lag.PCTOWNHOME -0.140681 0.058529 -2.4036 0.016234
#>
#> Rho: 0.026981, LR test value: 2.558, p-value: 0.10974
#> Asymptotic standard error: 0.016766
#> z-value: 1.6093, p-value: 0.10755
#> Wald statistic: 2.5899, p-value: 0.10755
#>
#> Log likelihood: -269.1031 for mixed model
#> ML residual variance (sigma squared): 0.39587, (sigma: 0.62918)
#> Number of observations: 281
#> Number of parameters estimated: 10
#> AIC: 558.21, (AIC for lm: 558.76)
#> LM test for residual autocorrelation
#> test value: 4.908, p-value: 0.026732
#>
res <- MCMCsamp(lsar1, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> rho 0.02683 0.01652 0.0002336 0.001222
#> (Intercept) -1.12357 0.25891 0.0036615 0.020126
#> PEXPOSURE 0.21348 0.08108 0.0011467 0.006288
#> PCTAGE65P 3.36104 0.68392 0.0096721 0.052730
#> PCTOWNHOME 0.09177 0.26971 0.0038143 0.021975
#> lag.(Intercept) 0.14089 0.06108 0.0008638 0.005039
#> lag.PEXPOSURE -0.03549 0.01605 0.0002270 0.001229
#> lag.PCTAGE65P 0.13688 0.23196 0.0032805 0.018439
#> lag.PCTOWNHOME -0.15040 0.05940 0.0008400 0.004527
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> rho -0.005181 0.01561 0.02730 0.03756 0.059973
#> (Intercept) -1.618594 -1.29277 -1.12902 -0.95147 -0.583037
#> PEXPOSURE 0.060313 0.15941 0.21675 0.26835 0.378858
#> PCTAGE65P 2.091132 2.87639 3.37052 3.84396 4.618430
#> PCTOWNHOME -0.439832 -0.09099 0.09702 0.28782 0.615503
#> lag.(Intercept) 0.016612 0.10280 0.14452 0.18005 0.256387
#> lag.PEXPOSURE -0.066728 -0.04674 -0.03595 -0.02417 -0.006011
#> lag.PCTAGE65P -0.306106 -0.02674 0.12948 0.29879 0.606519
#> lag.PCTOWNHOME -0.270533 -0.18772 -0.14930 -0.11312 -0.041099
#>
ssar0 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY)
summary(ssar0)
#>
#> Call:sacsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.468382 -0.375687 -0.034996 0.314714 3.833950
#>
#> Type: sac
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.386572 0.123188 -3.1381 0.001701
#> PEXPOSURE 0.026684 0.024013 1.1112 0.266479
#> PCTAGE65P 3.089824 0.562851 5.4896 4.029e-08
#> PCTOWNHOME -0.323052 0.137449 -2.3503 0.018756
#>
#> Rho: 0.089451
#> Asymptotic standard error: 0.019427
#> z-value: 4.6046, p-value: 4.1325e-06
#> Lambda: -0.08192
#> Asymptotic standard error: 0.033201
#> z-value: -2.4674, p-value: 0.01361
#>
#> LR test value: 10.114, p-value: 0.0063661
#>
#> Log likelihood: -273.672 for sac model
#> ML residual variance (sigma squared): 0.3766, (sigma: 0.61368)
#> Number of observations: 281
#> Number of parameters estimated: 7
#> AIC: 561.34, (AIC for lm: 567.46)
#>
res <- MCMCsamp(ssar0, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> rho 0.003025 0.07120 0.001007 0.01775
#> lambda 0.023449 0.07639 0.001080 0.01584
#> (Intercept) -0.633518 0.27623 0.003906 0.03590
#> PEXPOSURE 0.078740 0.07188 0.001017 0.01103
#> PCTAGE65P 3.447683 0.64893 0.009177 0.04671
#> PCTOWNHOME -0.323160 0.20883 0.002953 0.01820
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> rho -0.12891 -0.06108 0.01500 0.06611 0.1038
#> lambda -0.10782 -0.04264 0.02638 0.09726 0.1322
#> (Intercept) -1.23395 -0.83597 -0.58227 -0.43465 -0.2052
#> PEXPOSURE -0.02448 0.02465 0.05975 0.12415 0.2415
#> PCTAGE65P 2.26394 3.00999 3.42170 3.90465 4.7714
#> PCTOWNHOME -0.72455 -0.46375 -0.32118 -0.19006 0.1020
#>
ssar1 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, type="sacmixed")
summary(ssar1)
#>
#> Call:sacsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
#> listw = listw_NY, type = "sacmixed")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.633958 -0.363826 -0.019927 0.348238 3.655509
#>
#> Type: sacmixed
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.133298 0.247495 -4.5791 4.670e-06
#> PEXPOSURE 0.206963 0.074480 2.7788 0.005456
#> PCTAGE65P 3.083983 0.671081 4.5955 4.316e-06
#> PCTOWNHOME 0.174800 0.256280 0.6821 0.495196
#> lag.(Intercept) 0.153427 0.050817 3.0192 0.002534
#> lag.PEXPOSURE -0.033400 0.013817 -2.4173 0.015634
#> lag.PCTAGE65P -0.079738 0.222144 -0.3589 0.719634
#> lag.PCTOWNHOME -0.102502 0.056760 -1.8059 0.070940
#>
#> Rho: 0.092495
#> Asymptotic standard error: 0.023829
#> z-value: 3.8817, p-value: 0.00010375
#> Lambda: -0.091069
#> Asymptotic standard error: 0.038431
#> z-value: -2.3697, p-value: 0.017804
#>
#> LR test value: 22.379, p-value: 0.0010335
#>
#> Log likelihood: -267.5392 for sacmixed model
#> ML residual variance (sigma squared): 0.35617, (sigma: 0.5968)
#> Number of observations: 281
#> Number of parameters estimated: 11
#> AIC: 557.08, (AIC for lm: 567.46)
#>
res <- MCMCsamp(ssar1, mcmc=5000, burnin=500, listw=listw_NY)
summary(res)
#>
#> Iterations = 1:5000
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> rho -0.01142 0.06745 0.0009538 0.014884
#> lambda 0.02714 0.06616 0.0009357 0.012108
#> (Intercept) -1.06996 0.24244 0.0034287 0.019867
#> PEXPOSURE 0.21671 0.07858 0.0011113 0.006317
#> PCTAGE65P 3.28772 0.66128 0.0093519 0.052019
#> PCTOWNHOME 0.01692 0.25217 0.0035663 0.020679
#> lag.(Intercept) 0.11616 0.06420 0.0009079 0.006229
#> lag.PEXPOSURE -0.03494 0.01674 0.0002367 0.001567
#> lag.PCTAGE65P 0.26774 0.29277 0.0041405 0.030786
#> lag.PCTOWNHOME -0.15841 0.07020 0.0009927 0.007588
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> rho -0.12817 -0.06749 -0.0096060 0.04658 0.104484
#> lambda -0.11365 -0.02404 0.0381424 0.08516 0.119884
#> (Intercept) -1.56301 -1.21034 -1.0493309 -0.90101 -0.620537
#> PEXPOSURE 0.05671 0.16419 0.2144963 0.27104 0.380724
#> PCTAGE65P 1.90152 2.83869 3.3276412 3.73813 4.673860
#> PCTOWNHOME -0.45412 -0.17097 0.0003399 0.18841 0.536266
#> lag.(Intercept) -0.02415 0.07789 0.1131960 0.16018 0.241134
#> lag.PEXPOSURE -0.06796 -0.04608 -0.0356251 -0.02399 -0.002681
#> lag.PCTAGE65P -0.25665 0.05732 0.2597481 0.44511 0.926273
#> lag.PCTOWNHOME -0.29581 -0.20585 -0.1603293 -0.10924 -0.024885
#>
# }