+
-
+
+
+
+
+
+Mortality rates and other health statistics for counties are, in many
+cases, highly unstable estimates that cannot be relied upon for public
+advisories or inference (due to small population sizes). Hence, we need
+models to make inferences from small area data.
+
The following code fits a spatial conditional autoregressive (CAR) model
to female county mortality data. These models are used for estimating
disease risk in small areas like counties, and for analyzing covariation
-of health outcomes with other area qualities. The R syntax for fitting
+of health outcomes with other area variables. The R syntax for fitting
the models is similar to using `lm` or `glm`. We provide the population
at risk (the denominator for mortality rates) as an offset term, using
-the log-transform. In this case, three of the observations are missing
+the log-transform.
+
+In our Georgia mortality data, three of the observations are missing
because they have been censored; per CDC criteria, this means that there
were 9 or fewer deaths in those counties. By using the `censor_point`
-argument and setting it to `censor_point = 9`, the model will account
-for the censoring process when providing estimates of the mortality
-rates:
+argument and setting it to `censor_point = 9`, we can easily obtain
+estimates for the censored counties (along with all the others) using
+models account for the censoring process:
``` r
-cars <- prep_car_data(A)
+# prepare a list of data for CAR models in Stan
+cars <- prep_car_data(C)
#> Range of permissible rho values: -1.661, 1
+
+# fit the model to female mortality rates
fit <- stan_car(deaths.female ~ offset(log(pop.at.risk.female)),
censor_point = 9,
data = georgia,
car_parts = cars,
family = poisson(),
- cores = 4, # for multi-core processing
- refresh = 0) # to silence some printing
+ iter = 1e3, # no. MCMC samples
+ quiet = TRUE) # to silence printing
#> 3 NA values identified in the outcome variable
#> Found in rows: 55, 126, 157
-#>
-#> *Setting prior parameters for intercept
-#> Distribution: normal
-#> location scale
-#> 1 -4.7 5
-#>
-#> *Setting prior for CAR scale parameter (car_scale)
-#> Distribution: student_t
-#> df location scale
-#> 1 10 0 3
-#>
-#> *Setting prior for CAR spatial autocorrelation parameter (car_rho)
-#> Distribution: uniform
-#> lower upper
-#> 1 -1.7 1
```
Passing a fitted model to the `sp_diag` function will return a set of
diagnostics for spatial models:
``` r
-sp_diag(fit, georgia, w = A)
-#> Using sp_diag(y, shape, rates = TRUE, ...). To examine data as (unstandardized) counts, use rates = FALSE.
-#> 3 NA values found in x will be dropped from data x and matrix w
-#> Warning: Removed 3 rows containing missing values or values
-#> outside the scale range (`geom_pointrange()`).
+sp_diag(fit, georgia)
+#> 3 NA values found in x will be dropped from data x and from matrix w (nb: this disrupts row-standardization of w)
+#> Warning: Removed 3 rows containing missing values or values outside the
+#> scale range (`geom_pointrange()`).
```
@@ -200,58 +205,149 @@ diagnostics from Stan (Monte Carlo standard errors of the mean
print(fit)
#> Spatial Model Results
#> Formula: deaths.female ~ offset(log(pop.at.risk.female))
-#> Spatial method (outcome): CAR
-#> Likelihood function: poisson
-#> Link function: log
-#> Residual Moran Coefficient: 0.0011525
-#> WAIC: 1227.47
+#> Likelihood: poisson
+#> Link: log
+#> Spatial method: CAR
+#> Residual Moran Coefficient: -0.0031375
#> Observations: 156
-#> Data models (ME): none
+#>
#> Inference for Stan model: foundation.
-#> 4 chains, each with iter=2000; warmup=1000; thin=1;
-#> post-warmup draws per chain=1000, total post-warmup draws=4000.
+#> 4 chains, each with iter=1000; warmup=500; thin=1;
+#> post-warmup draws per chain=500, total post-warmup draws=2000.
#>
#> mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
-#> intercept -4.674 0.002 0.089 -4.849 -4.730 -4.674 -4.621 -4.505 2362 1.000
-#> car_rho 0.923 0.001 0.058 0.778 0.879 0.937 0.973 0.995 3319 1.000
-#> car_scale 0.458 0.001 0.036 0.395 0.428 0.456 0.488 0.534 3618 0.999
+#> intercept -4.677 0.003 0.092 -4.871 -4.734 -4.677 -4.619 -4.513 861 1.004
+#> car_rho 0.924 0.001 0.058 0.784 0.883 0.936 0.973 0.996 1606 0.999
+#> car_scale 0.456 0.001 0.036 0.391 0.424 0.454 0.486 0.532 1899 1.002
#>
-#> Samples were drawn using NUTS(diag_e) at Tue Sep 17 16:44:56 2024.
+#> Samples were drawn using NUTS(diag_e) at Wed Nov 13 16:09:50 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
```
-Applying the `fitted` method to the fitted model will return the fitted
-values from the model - in this case, the fitted values are the
-estimates of the county mortality rates. Multiplying them by 10,000
+To extract estimates of the county mortality rates from this, we apply
+the `fitted` method - in this case, the fitted values from the model are
+the estimates of the county mortality rates. Multiplying them by 10,000
gives mortality rate per 10,000 at risk:
``` r
+# mortality rates per 10,000 at risk
mortality_est <- fitted(fit) * 10e3
+
+# display rates with county names
county_name <- georgia$NAME
head( cbind(county_name, mortality_est) )
#> county_name mean sd 2.5% 20% 50%
-#> fitted[1] Crisp 101.48785 9.604829 83.99009 93.31163 101.17610
-#> fitted[2] Candler 136.99885 15.905146 109.27395 123.11823 136.31355
-#> fitted[3] Barrow 94.25470 6.071597 82.80270 89.20105 94.16678
-#> fitted[4] DeKalb 59.76214 1.579194 56.72962 58.44624 59.75766
-#> fitted[5] Columbia 53.33958 3.257549 47.19615 50.56654 53.28387
-#> fitted[6] Cobb 54.12983 1.498260 51.24933 52.85101 54.10133
+#> fitted[1] Crisp 101.89147 9.784748 83.87184 93.37621 101.37227
+#> fitted[2] Candler 137.13841 15.643722 108.50033 123.57906 136.51359
+#> fitted[3] Barrow 94.35971 6.364805 82.62612 88.73920 94.14574
+#> fitted[4] DeKalb 59.74315 1.575741 56.77068 58.33325 59.71677
+#> fitted[5] Columbia 53.34581 3.207432 47.33439 50.61504 53.26339
+#> fitted[6] Cobb 54.12259 1.495041 51.24262 52.86109 54.12912
#> 80% 97.5%
-#> fitted[1] 109.30723 121.16598
-#> fitted[2] 150.17348 169.77611
-#> fitted[3] 99.19399 106.44508
-#> fitted[4] 61.07091 62.86805
-#> fitted[5] 56.08790 59.78086
-#> fitted[6] 55.42278 57.02966
+#> fitted[1] 110.47617 122.04006
+#> fitted[2] 150.26114 169.29720
+#> fitted[3] 99.74677 107.51136
+#> fitted[4] 61.13469 62.84502
+#> fitted[5] 56.05293 59.83064
+#> fitted[6] 55.38574 57.00559
```
The mortality estimates are stored in the column named “mean”, and the
limits of the 95% credible interval are found in the columns “2.5%” and
-“97.5%”.
+“97.5%”. Here we create a map of estimates (with some help from `sf`
+package):
+
+``` r
+library(sf)
+
+# put estimates into bins for map colors
+x <- mortality_est$mean
+brks <- quantile(x, probs = c(0, 0.2, 0.4, 0.6, 0.8, 1))
+est_cut <- cut(x, breaks = brks, include.lowest = TRUE)
+
+# assign colors to values
+rank <- as.numeric( est_cut )
+pal_fun <- colorRampPalette( c("#5D74A5FF", "gray90", "#A8554EFF") )
+pal <- pal_fun( max(rank) )
+colors <- pal[ rank ]
+
+# set plot margins
+og=par(mar=rep(1, 4))
+
+# get boundaries
+geom <- sf::st_geometry(georgia)
+
+# map estimates
+plot(geom,
+ lwd = 0.2,
+ col = colors)
+
+# legend
+legend("right",
+ fill = pal,
+ title = 'Mortality per 10,000',
+ legend = levels(est_cut),
+ bty = 'n'
+)
+
+mtext('County mortality rates per 10,000, Georgia women ages 55-64', side = 3, font = 2)
+```
+
+
+
+``` r
+# reset margins
+par(og)
+```
+
+Using the credible intervals, we can complement our map with a
+point-interval plot:
+
+``` r
+# order counties by mortality rate
+index <- order(mortality_est$mean, decreasing = TRUE)
+dat <- mortality_est[index, ]
+
+# gather estimate with credible interval (95%)
+est <- dat$mean
+lwr <- dat$`2.5%`
+upr <- dat$`97.5%`
+y <- seq_along(county_name)
+x_lim <- c(min(lwr), max(upr)) |>
+ round()
+
+og=par(mar = c(3, 0, 0, 0))
+
+# points
+plot(est,
+ y,
+ pch = 5,
+ col = 'gray50',
+ bty = 'L',
+ axes = FALSE,
+ xlim = x_lim,
+ ylab = NA,
+ xlab = NA)
+
+# intervals
+segments(x0 = lwr, x1 = upr,
+ y0 = y, y1 = y,
+ col = colors[ index ])
+
+# x axis
+axis(1, at = seq(x_lim[1], x_lim[2], by = 20))
+mtext('County mortality rates per 10,000, Georgia women ages 55-64', side = 1, line = 2)
+```
+
+
+
+``` r
+par(og)
+```
-Details and demonstrations can be found in the package [help
+More details and demonstrations can be found in the package [help
pages](https://connordonegan.github.io/geostan/reference/index.html) and
[vignettes](https://connordonegan.github.io/geostan/articles/index.html).
diff --git a/docs/index.html b/docs/index.html
index 3592e456..34cb49ad 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -85,7 +85,7 @@ - -Disease mapping and spatial regression Statistical models for data recorded across areal units like states, counties, or census tracts. +Spatial regression and disease mapping Statistical models for data recorded across areal units (states, counties, or census tracts) or networks, including spatial econometric models.
-
Spatial analysis tools Tools for visualizing and measuring spatial autocorrelation and map patterns, for exploratory analysis and model diagnostics.
@@ -96,7 +96,7 @@ - -Custom spatial models Tools for building custom spatial models in Stan. +Custom spatial models Tools for building custom spatial or network models in Stan.
geostan: Bayesian spatial analysis
The RStan ecosystem Interfaces easily with many high-quality R packages for Bayesian modeling.
For public health research, geostan complements the surveil R package for the study of time trends in disease incidence or mortality data.
@@ -142,10 +142,12 @@Usage -
This has county population and mortality data by sex for ages 55-64, and for the period 2014-2018. As is common for public access data, some of the observations missing because the CDC has censored them.
-The sp_diag function provides visual summaries of spatial data, including a histogram, Moran scatter plot, and map. Here is a visual summary of crude female mortality rates (as deaths per 10,000):
This has county population and mortality data by sex for ages 55-64, and for the period 2014-2018. As is common for public access data, some of the observations missing because the CDC has censored them to protect privacy.
+The sp_diag function provides visual summaries of spatial data, including a histogram, Moran scatter plot, and map. The Moran scatter plot displays the values against a summary of their neighboring values, so that the slope of the regression line gives a measure of their degree of autocorrelation.
Here is a quick visual summary of crude female mortality rates (as deaths per 10,000):
-A <- shape2mat(georgia, style = "B")
+# create adjacency matrix ("B" is for binary)
+C <- shape2mat(georgia, style = "B")
#> Contiguity condition: queen
#> Number of neighbors per unit, summary:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
@@ -154,94 +156,165 @@ Usage
#> Spatial weights, summary:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1 1 1 1 1 1
+
+# crude mortality rate per 10,000 at risk
mortality_rate <- georgia$rate.female * 10e3
-sp_diag(mortality_rate, georgia, w = A)
-#> 3 NA values found in x will be dropped from data x and matrix w
-#> Warning: Removed 3 rows containing non-finite outside the scale
-#> range (`stat_bin()`).
The following code fits a spatial conditional autoregressive (CAR) model to female county mortality data. These models are used for estimating disease risk in small areas like counties, and for analyzing covariation of health outcomes with other area qualities. The R syntax for fitting the models is similar to using lm or glm. We provide the population at risk (the denominator for mortality rates) as an offset term, using the log-transform. In this case, three of the observations are missing because they have been censored; per CDC criteria, this means that there were 9 or fewer deaths in those counties. By using the censor_point argument and setting it to censor_point = 9, the model will account for the censoring process when providing estimates of the mortality rates:
Mortality rates and other health statistics for counties are, in many cases, highly unstable estimates that cannot be relied upon for public advisories or inference (due to small population sizes). Hence, we need models to make inferences from small area data.
+The following code fits a spatial conditional autoregressive (CAR) model to female county mortality data. These models are used for estimating disease risk in small areas like counties, and for analyzing covariation of health outcomes with other area variables. The R syntax for fitting the models is similar to using lm or glm. We provide the population at risk (the denominator for mortality rates) as an offset term, using the log-transform.
In our Georgia mortality data, three of the observations are missing because they have been censored; per CDC criteria, this means that there were 9 or fewer deaths in those counties. By using the censor_point argument and setting it to censor_point = 9, we can easily obtain estimates for the censored counties (along with all the others) using models account for the censoring process:
-cars <- prep_car_data(A)
+# prepare a list of data for CAR models in Stan
+cars <- prep_car_data(C)
#> Range of permissible rho values: -1.661, 1
+
+# fit the model to female mortality rates
fit <- stan_car(deaths.female ~ offset(log(pop.at.risk.female)),
censor_point = 9,
data = georgia,
car_parts = cars,
family = poisson(),
- cores = 4, # for multi-core processing
- refresh = 0) # to silence some printing
+ iter = 1e3, # no. MCMC samples
+ quiet = TRUE) # to silence printing
#> 3 NA values identified in the outcome variable
-#> Found in rows: 55, 126, 157
-#>
-#> *Setting prior parameters for intercept
-#> Distribution: normal
-#> location scale
-#> 1 -4.7 5
-#>
-#> *Setting prior for CAR scale parameter (car_scale)
-#> Distribution: student_t
-#> df location scale
-#> 1 10 0 3
-#>
-#> *Setting prior for CAR spatial autocorrelation parameter (car_rho)
-#> Distribution: uniform
-#> lower upper
-#> 1 -1.7 1Passing a fitted model to the sp_diag function will return a set of diagnostics for spatial models:
-sp_diag(fit, georgia, w = A)
-#> Using sp_diag(y, shape, rates = TRUE, ...). To examine data as (unstandardized) counts, use rates = FALSE.
-#> 3 NA values found in x will be dropped from data x and matrix w
-#> Warning: Removed 3 rows containing missing values or values
-#> outside the scale range (`geom_pointrange()`).sp_diag(fit, georgia)
+#> 3 NA values found in x will be dropped from data x and from matrix w (nb: this disrupts row-standardization of w)
+#> Warning: Removed 3 rows containing missing values or values outside the
+#> scale range (`geom_pointrange()`).
The print method returns a summary of the probability distributions for model parameters, as well as Markov chain Monte Carlo (MCMC) diagnostics from Stan (Monte Carlo standard errors of the mean se_mean, effective sample size n_eff, and the R-hat statistic Rhat):
print(fit)
#> Spatial Model Results
#> Formula: deaths.female ~ offset(log(pop.at.risk.female))
-#> Spatial method (outcome): CAR
-#> Likelihood function: poisson
-#> Link function: log
-#> Residual Moran Coefficient: 0.0011525
-#> WAIC: 1227.47
+#> Likelihood: poisson
+#> Link: log
+#> Spatial method: CAR
+#> Residual Moran Coefficient: -0.0031375
#> Observations: 156
-#> Data models (ME): none
+#>
#> Inference for Stan model: foundation.
-#> 4 chains, each with iter=2000; warmup=1000; thin=1;
-#> post-warmup draws per chain=1000, total post-warmup draws=4000.
+#> 4 chains, each with iter=1000; warmup=500; thin=1;
+#> post-warmup draws per chain=500, total post-warmup draws=2000.
#>
#> mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
-#> intercept -4.674 0.002 0.089 -4.849 -4.730 -4.674 -4.621 -4.505 2362 1.000
-#> car_rho 0.923 0.001 0.058 0.778 0.879 0.937 0.973 0.995 3319 1.000
-#> car_scale 0.458 0.001 0.036 0.395 0.428 0.456 0.488 0.534 3618 0.999
+#> intercept -4.677 0.003 0.092 -4.871 -4.734 -4.677 -4.619 -4.513 861 1.004
+#> car_rho 0.924 0.001 0.058 0.784 0.883 0.936 0.973 0.996 1606 0.999
+#> car_scale 0.456 0.001 0.036 0.391 0.424 0.454 0.486 0.532 1899 1.002
#>
-#> Samples were drawn using NUTS(diag_e) at Tue Sep 17 16:44:56 2024.
+#> Samples were drawn using NUTS(diag_e) at Wed Nov 13 16:09:50 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).Applying the fitted method to the fitted model will return the fitted values from the model - in this case, the fitted values are the estimates of the county mortality rates. Multiplying them by 10,000 gives mortality rate per 10,000 at risk:
To extract estimates of the county mortality rates from this, we apply the fitted method - in this case, the fitted values from the model are the estimates of the county mortality rates. Multiplying them by 10,000 gives mortality rate per 10,000 at risk:
-mortality_est <- fitted(fit) * 10e3
+# mortality rates per 10,000 at risk
+mortality_est <- fitted(fit) * 10e3
+
+# display rates with county names
county_name <- georgia$NAME
head( cbind(county_name, mortality_est) )
#> county_name mean sd 2.5% 20% 50%
-#> fitted[1] Crisp 101.48785 9.604829 83.99009 93.31163 101.17610
-#> fitted[2] Candler 136.99885 15.905146 109.27395 123.11823 136.31355
-#> fitted[3] Barrow 94.25470 6.071597 82.80270 89.20105 94.16678
-#> fitted[4] DeKalb 59.76214 1.579194 56.72962 58.44624 59.75766
-#> fitted[5] Columbia 53.33958 3.257549 47.19615 50.56654 53.28387
-#> fitted[6] Cobb 54.12983 1.498260 51.24933 52.85101 54.10133
+#> fitted[1] Crisp 101.89147 9.784748 83.87184 93.37621 101.37227
+#> fitted[2] Candler 137.13841 15.643722 108.50033 123.57906 136.51359
+#> fitted[3] Barrow 94.35971 6.364805 82.62612 88.73920 94.14574
+#> fitted[4] DeKalb 59.74315 1.575741 56.77068 58.33325 59.71677
+#> fitted[5] Columbia 53.34581 3.207432 47.33439 50.61504 53.26339
+#> fitted[6] Cobb 54.12259 1.495041 51.24262 52.86109 54.12912
#> 80% 97.5%
-#> fitted[1] 109.30723 121.16598
-#> fitted[2] 150.17348 169.77611
-#> fitted[3] 99.19399 106.44508
-#> fitted[4] 61.07091 62.86805
-#> fitted[5] 56.08790 59.78086
-#> fitted[6] 55.42278 57.02966The mortality estimates are stored in the column named “mean”, and the limits of the 95% credible interval are found in the columns “2.5%” and “97.5%”.
-Details and demonstrations can be found in the package help pages and vignettes.
+#> fitted[1] 110.47617 122.04006 +#> fitted[2] 150.26114 169.29720 +#> fitted[3] 99.74677 107.51136 +#> fitted[4] 61.13469 62.84502 +#> fitted[5] 56.05293 59.83064 +#> fitted[6] 55.38574 57.00559 +The mortality estimates are stored in the column named “mean”, and the limits of the 95% credible interval are found in the columns “2.5%” and “97.5%”. Here we create a map of estimates (with some help from sf package):
+library(sf)
+
+# put estimates into bins for map colors
+x <- mortality_est$mean
+brks <- quantile(x, probs = c(0, 0.2, 0.4, 0.6, 0.8, 1))
+est_cut <- cut(x, breaks = brks, include.lowest = TRUE)
+
+# assign colors to values
+rank <- as.numeric( est_cut )
+pal_fun <- colorRampPalette( c("#5D74A5FF", "gray90", "#A8554EFF") )
+pal <- pal_fun( max(rank) )
+colors <- pal[ rank ]
+
+# set plot margins
+og=par(mar=rep(1, 4))
+
+# get boundaries
+geom <- sf::st_geometry(georgia)
+
+# map estimates
+plot(geom,
+ lwd = 0.2,
+ col = colors)
+
+# legend
+legend("right",
+ fill = pal,
+ title = 'Mortality per 10,000',
+ legend = levels(est_cut),
+ bty = 'n'
+)
+
+mtext('County mortality rates per 10,000, Georgia women ages 55-64', side = 3, font = 2)
+# reset margins
+par(og)Using the credible intervals, we can complement our map with a point-interval plot:
+
+# order counties by mortality rate
+index <- order(mortality_est$mean, decreasing = TRUE)
+dat <- mortality_est[index, ]
+
+# gather estimate with credible interval (95%)
+est <- dat$mean
+lwr <- dat$`2.5%`
+upr <- dat$`97.5%`
+y <- seq_along(county_name)
+x_lim <- c(min(lwr), max(upr)) |>
+ round()
+
+og=par(mar = c(3, 0, 0, 0))
+
+# points
+plot(est,
+ y,
+ pch = 5,
+ col = 'gray50',
+ bty = 'L',
+ axes = FALSE,
+ xlim = x_lim,
+ ylab = NA,
+ xlab = NA)
+
+# intervals
+segments(x0 = lwr, x1 = upr,
+ y0 = y, y1 = y,
+ col = colors[ index ])
+
+# x axis
+axis(1, at = seq(x_lim[1], x_lim[2], by = 20))
+mtext('County mortality rates per 10,000, Georgia women ages 55-64', side = 1, line = 2)
+par(og)More details and demonstrations can be found in the package help pages and vignettes.
Citing geostan
@@ -249,16 +322,16 @@ Citing geostanIf you use geostan in published work, please include a citation.
Donegan, Connor (2022) “geostan: An R package for Bayesian spatial analysis” The Journal of Open Source Software. 7, no. 79: 4716. https://doi.org/10.21105/joss.04716.
-@Article{,
- title = {{geostan}: An {R} package for {B}ayesian spatial analysis},
- author = {Connor Donegan},
- journal = {The Journal of Open Source Software},
- year = {2022},
- volume = {7},
- number = {79},
- pages = {4716},
- doi = {10.21105/joss.04716},
-}
+@Article{,
+ title = {{geostan}: An {R} package for {B}ayesian spatial analysis},
+ author = {Connor Donegan},
+ journal = {The Journal of Open Source Software},
+ year = {2022},
+ volume = {7},
+ number = {79},
+ pages = {4716},
+ doi = {10.21105/joss.04716},
+}
Donegan, Connor (2022) “geostan: An R package for Bayesian spatial analysis” The Journal of Open Source Software. 7, no. 79: 4716. https://doi.org/10.21105/joss.04716.
@Article{,
- title = {{geostan}: An {R} package for {B}ayesian spatial analysis},
- author = {Connor Donegan},
- journal = {The Journal of Open Source Software},
- year = {2022},
- volume = {7},
- number = {79},
- pages = {4716},
- doi = {10.21105/joss.04716},
-}@Article{,
+ title = {{geostan}: An {R} package for {B}ayesian spatial analysis},
+ author = {Connor Donegan},
+ journal = {The Journal of Open Source Software},
+ year = {2022},
+ volume = {7},
+ number = {79},
+ pages = {4716},
+ doi = {10.21105/joss.04716},
+}Spillover/diffusion effects for spatial lag models
spill(beta, gamma = 0, rho, W, method = c("quick", "proper"), K = 20)
-impacts(object, method = c("quick", "proper"), K = 20)Arguments
A fitted spatial lag model (from stan_sar)
An object of class 'impacts_slm', as returned by geostan::impacts
Round results to this many digits
Additional arguments will be passed to base::print
Details
diff --git a/docs/reference/index.html b/docs/reference/index.html index b0187c1e..71b8cb24 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -210,7 +210,7 @@Methods and diagnostics spill() impacts()
+
spill() impacts() print(<impacts_slm>)
Spillover/diffusion effects for spatial lag models
Predict method for geostan_fit models
summary = TRUE,
type = c("link", "response"),
add_slx = FALSE,
+ approx = FALSE,
+ K = 15,
...
)
@@ -105,6 +107,14 @@ Arguments
Logical. If add_slx = TRUE, any spatially-lagged covariates that were specified through the 'slx' argument (of the model fitting function, e.g., stan_glm) will be added to the linear predictor. The spatial lag terms will be calculated internally using object$C, the spatial weights matrix used to fit the model. Hence, newdata must have N = object$N rows. Predictions from spatial lag models (SAR models of type 'SLM' and 'SDLM') always include the SLX terms (i.e., any value passed to add_slx will be overwritten with TRUE).
For SAR models of type 'SLM' or 'SDLM' only; use an approximation for matrix inversion? See details below.
Number of matrix powers to use with approx.
Not used
Value
Details
The primary purpose of the predict method is to explore marginal effects of covariates.
-The model formula will be taken from object$formula, and then a model matrix will be created by passing newdata to the model.frame function (as in: model.frame(newdata, object$formula). Parameters are taken from as.matrix(object, pars = c("intercept", "beta")). If add_slx = TRUE, SLX coefficients will be taken from as.matrix(object, pars = "gamma").
Spatially-lagged covariates added via the slx argument ('spillover effects') will be included if add_slx = TRUE or if a spatial lag model is provided (a SAR model of type 'SLM' or 'SDLM'). In either of those cases, newdata must have the same number of rows as were used to fit the original data. For details on these 'spillover effects', see LeSage and Pace (2009) and LeSage (2014).
The model formula will be taken from object$formula, and then a model matrix will be created by passing newdata to the model.frame function (as in: model.frame(object$formula, newdata). Parameters are taken from as.matrix(object, pars = c("intercept", "beta")).
Spatial lag of X
+ + +Spatially-lagged covariates which were included via the slx argument will, by default, not be included. They will be be included in predictions if add_slx = TRUE or if the fitted model is a SAR model of type 'SLM' or 'SDLM'. In either of those cases, newdata must have the same number of rows as were used to fit the original data.
Spatial lag of Y
+ + +The typical 'marginal effect' interpretation of the regression coefficients does not hold for the SAR models of type 'SLM' or 'SDLM'. For details on these 'spillover effects', see LeSage and Pace (2009), LeSage (2014), and geostan::impacts.
Predictions for the spatial lag model (SAR models of type 'SLM') are equal to: $$ (I - \rho W)^{-1} X \beta $$ -where \(X \beta\) contains the intercept and covariates. (For intercept-only models, the above term is equal to the constant intercept.) Predictions for the spatial Durbin lag model (SAR models of type 'SDLM') are equal to: +where \(X \beta\) contains the intercept and covariates. Predictions for the spatial Durbin lag model (SAR models of type 'SDLM') are equal to: $$ (I - \rho W)^{-1} (X \beta + WX \gamma) $$ -where \(WX \gamma\) are spatially lagged covariates multiplied by their coefficients. The SLM and SDLM differ from all other model types in that the spatial component of the model cannot be separated from the linear predictor and is, therefore, automatically incorporated into the predictions.
-In generalized linear models (such as Poisson and Binomial models) marginal effects plots on the response scale may be sensitive to the level of other covariates in the model and to location. If the model includes a spatial autocorrelation component (for example, you used a spatial CAR, SAR (error model), or ESF model), by default these terms will be fixed at zero for the purposes of calculating marginal effects. If you want to change this, you can introduce spatial trend values by specifying a varying intercept using the alpha argument.
The inverse of the matrix \((I - \rho W)\) can be time consuming to compute (especially when iterating over MCMC samples). You can use approx = TRUE to approximate the inverse using a series of matrix powers. The argument \(K\) controls how many powers to use for the approximation. As a rule, higher values of \(\rho\) require larger \(K\). Notice that \(\rho^K\) should be close to zero for the approximation to hold. For example, for \(\rho = .5\) a value of \(K=8\) may suffice (eqn0.5^8 = 0.004), but larger values of \(\rho\) require higher values of \(K\).
Generalized linear models
+ + +In generalized linear models (such as Poisson and Binomial models) marginal effects plots on the response scale may be sensitive to the level of other covariates in the model and to geographic location. If the model includes a spatial autocorrelation component (for example, you used a spatial CAR, SAR, or ESF model), by default these terms will be fixed at zero for the purposes of calculating marginal effects. If you want to change this, you can introduce spatial trend values by specifying a varying intercept using the alpha argument.