diff --git a/DESCRIPTION b/DESCRIPTION index b6099d7e..6a1b11d8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: geostan Title: Bayesian Spatial Analysis -Version: 0.8.0 -Date: 2024-11-15 +Version: 0.8.1 +Date: 2024-12-04 URL: https://connordonegan.github.io/geostan/ BugReports: https://github.com/ConnorDonegan/geostan/issues Authors@R: c( diff --git a/NEWS.md b/NEWS.md index 2e383bbc..05645923 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,8 @@ # geostan 0.8.1 -This version changes one line in the 'Spatial analysis with geostan' vignette that caused the package to fail on installation for old versions of windows and mac os. +This is a small patch to ensure that the package can install on older versions of R. + +The patch changes one line in the 'Spatial analysis with geostan' vignette that caused the package to fail on installation for old versions of windows and mac os. # geostan 0.8.0 diff --git a/R/stan_glm.R b/R/stan_glm.R index 85315cd3..6a874cdf 100644 --- a/R/stan_glm.R +++ b/R/stan_glm.R @@ -205,6 +205,7 @@ #' # example prior for two covariates #' pl <- list(beta = normal(c(0, 0), #' c(1, 1)) +#' ) #' #' ## #' ## Poisson model for count data diff --git a/cran-comments.md~ b/cran-comments.md~ deleted file mode 100644 index 06c56051..00000000 --- a/cran-comments.md~ +++ /dev/null @@ -1 +0,0 @@ -Passed devtools::check and win-builder \ No newline at end of file diff --git a/docs/404.html b/docs/404.html index 9adc22b2..f99c8fd5 100644 --- a/docs/404.html +++ b/docs/404.html @@ -39,7 +39,7 @@
diff --git a/docs/LICENSE-text.html b/docs/LICENSE-text.html index ef459082..32f4286e 100644 --- a/docs/LICENSE-text.html +++ b/docs/LICENSE-text.html @@ -17,7 +17,7 @@ diff --git a/docs/LICENSE.html b/docs/LICENSE.html index 110a0745..9d54e0b6 100644 --- a/docs/LICENSE.html +++ b/docs/LICENSE.html @@ -17,7 +17,7 @@ diff --git a/docs/articles/custom-spatial-models.html b/docs/articles/custom-spatial-models.html index ddefca07..cf441241 100644 --- a/docs/articles/custom-spatial-models.html +++ b/docs/articles/custom-spatial-models.html @@ -40,7 +40,7 @@ diff --git a/docs/articles/index.html b/docs/articles/index.html index 17104171..4e4946a4 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -17,7 +17,7 @@ diff --git a/docs/articles/measuring-sa.html b/docs/articles/measuring-sa.html index a9fcb883..a8642796 100644 --- a/docs/articles/measuring-sa.html +++ b/docs/articles/measuring-sa.html @@ -40,7 +40,7 @@ @@ -299,7 +299,10 @@For a summary of model results:
print(fit)
@@ -309,18 +312,18 @@ Model diagnostics#> Likelihood: poisson
#> Link: log
#> Spatial method: Exchangeable
-#> Residual Moran Coefficient: 0.023644
+#> Residual Moran Coefficient: 0.02472075
#> Observations: 159
#>
#> 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.
#>
-#> mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
-#> intercept -4.180 0.001 0.021 -4.221 -4.197 -4.180 -4.163 -4.140 444 1.005
-#> alpha_tau 0.247 0.000 0.016 0.219 0.234 0.246 0.261 0.281 4397 1.000
+#> mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
+#> intercept -4.181 0.001 0.022 -4.225 -4.200 -4.182 -4.162 -4.137 382 1.01
+#> alpha_tau 0.247 0.000 0.016 0.219 0.234 0.247 0.260 0.281 4579 1.00
#>
-#> Samples were drawn using NUTS(diag_e) at Sat Nov 16 10:37:03 2024.
+#> Samples were drawn using NUTS(diag_e) at Wed Dec 4 13:07:41 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).
fit <- stan_sar(y ~ z, data = grid, C = W)
-The stan_sar
function will take the spatial weights matrix W
and pass it through a function called prep_sar_data
which will calculate the eigenvalues of the spatial weights matrix using Matrix::Schur
. This step can be prohibitive for large data sets (e.g., \(N = 100,000\)).
The stan_sar
function will take the spatial weights matrix W
and pass it through a function called prep_sar_data
which will calculate the eigenvalues of the spatial weights matrix using Matrix::Schur
. This step can be prohibitive for large data sets (e.g., \(N = 100,000\)).
The following code would normally be used to fit a conditional autoregressive (CAR) model:
C <- shape2mat(grid, style = "B", queen = FALSE)
car_list <- prep_car_data(C, "WCAR")
fit <- stan_car(y ~ z, data = grid, car_parts = car_list)
Here, the prep_car_data
function calculates the eigenvalues of the spatial weights matrix using Matrix::Schur
, which is not feasible for large N.
Here, the prep_car_data
function calculates the eigenvalues of the spatial weights matrix using Matrix::Schur
, which is not feasible for large N.
The prep_sar_data2
and prep_car_data2
functions are designed for large raster layers. As input, they require the dimensions of the grid (number of rows and number of columns). The eigenvalues are produced very quickly using Equation 5 from Griffith (2000). The methods have some restrictions:
print(fit)
The approximate method uses matrix powers, and for high values of SA (especially values of rho greater than 0.9) the approximation is sensitive to the number of powers taken. In sim_sar
this is controlled by the argument K
. For rho = 0.9
, you may need to raise the default value of K for the approximation to hold. Notice that \(0.9^{45} = 0.009\), which is near zero and should provide sound results, whereas \(0.9^25 = 0.07\) which is not quite there, and \(0.9^{15} = 0.2\) will be completely inadequate. In geostan 0.8.0, the default is \(K = 20\).
fdf <- fitted(fit_lm)
head(fdf)
## mean sd 2.5% 20% 50% 80% 97.5%
-## fitted[1] 70.22671 0.3767095 69.49725 69.91060 70.21986 70.53322 71.00355
-## fitted[2] 63.48810 0.5714337 62.40300 62.99766 63.48573 63.94980 64.64769
-## fitted[3] 79.29656 0.6182122 78.11493 78.77971 79.28295 79.80699 80.47506
-## fitted[4] 80.31900 0.6666701 79.06651 79.74855 80.30555 80.86816 81.57068
-## fitted[5] 75.99801 0.4794856 75.09004 75.60299 75.99605 76.41061 76.94099
-## fitted[6] 67.89460 0.4132484 67.12661 67.55064 67.88613 68.23257 68.73600
+## fitted[1] 70.21184 0.3828826 69.45109 69.90147 70.20678 70.53523 70.97700
+## fitted[2] 63.46098 0.5897879 62.31106 62.95473 63.45715 63.95377 64.62738
+## fitted[3] 79.29817 0.5814509 78.17766 78.79123 79.29411 79.80983 80.42024
+## fitted[4] 80.32248 0.6268546 79.11845 79.78267 80.32545 80.86897 81.51848
+## fitted[5] 75.99363 0.4550348 75.13656 75.58926 76.00222 76.38572 76.87329
+## fitted[6] 67.87548 0.4279134 67.02802 67.51338 67.86547 68.23872 68.71756
The resid
method behaves similarly. Examining the Moran scatter plot using the residuals shows a moderate degree of positive SA as well as some skewness:
rdf <- resid(fit_lm)
@@ -385,7 +385,7 @@ Spatial regression## Likelihood: auto_gaussian
## Link: identity
## Spatial method: CAR
-## Residual Moran Coefficient: -0.3563787
+## Residual Moran Coefficient: -0.3565875
## Observations: 162
##
## Inference for Stan model: foundation.
@@ -393,11 +393,11 @@ Spatial regression## post-warmup draws per chain=400, total post-warmup draws=1600.
##
## mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
-## intercept 70.237 0.084 2.695 64.792 68.153 70.205 72.365 75.621 1027 1.000
-## car_rho 0.981 0.000 0.011 0.953 0.974 0.983 0.990 0.995 805 1.006
-## car_scale 7.827 0.014 0.458 7.017 7.433 7.803 8.215 8.779 1008 1.007
+## intercept 70.183 0.091 2.872 64.690 67.874 70.230 72.427 75.709 1005 1.001
+## car_rho 0.981 0.000 0.010 0.956 0.974 0.983 0.990 0.996 1111 0.999
+## car_scale 7.818 0.014 0.446 6.994 7.431 7.799 8.162 8.789 996 1.002
##
-## Samples were drawn using NUTS(diag_e) at Sat Nov 16 10:10:33 2024.
+## Samples were drawn using NUTS(diag_e) at Wed Dec 4 13:08:24 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).
@@ -441,13 +441,13 @@ The adjusted estimate of .59 is considerably different from the naive estimate of .80 and is outside the naive confidence intervals. (The adjusted estimate is .62 if we use SAR models.)
head(preds)
## gdpPercap mean sd 2.5% 20% 50% 80% 97.5%
-## 1 597.1352 60.14550 1.545635 57.40796 58.82435 60.12483 61.40456 63.43402
-## 2 1201.4715 62.94946 1.334939 60.58126 61.79015 62.91636 64.02657 65.75630
-## 3 1805.8079 64.58357 1.229620 62.38827 63.51346 64.54929 65.57171 67.11554
-## 4 2410.1442 65.74131 1.165277 63.67829 64.74130 65.70433 66.66620 68.15017
-## 5 3014.4805 66.63862 1.122363 64.62079 65.68522 66.58734 67.52480 68.97823
-## 6 3618.8169 67.37142 1.092381 65.43513 66.45664 67.33406 68.22549 69.66932
+## 1 597.1352 60.03988 1.538276 57.18041 58.76481 59.99106 61.31258 63.35620
+## 2 1201.4715 62.84811 1.332449 60.36180 61.72960 62.82398 63.91972 65.64284
+## 3 1805.8079 64.48471 1.229632 62.23882 63.45565 64.46347 65.46745 67.08958
+## 4 2410.1442 65.64422 1.166803 63.49242 64.66848 65.63881 66.56481 68.06477
+## 5 3014.4805 66.54289 1.124860 64.46880 65.62236 66.52958 67.42559 68.87927
+## 6 3618.8169 67.27680 1.095510 65.21997 66.37565 67.27810 68.11602 69.49700
These ‘predicted’ values represent the expectation of the outcome variable at the given level of the covariates. So we would expect actual observations to form a cloud of points around the ‘predicted’ values. To calculate these predicted values, the predict
function only includes covariates and the intercept, it does not include any spatial autocorrelation components. Its purpose is mainly to examine implications of the coefficient estimates on recognizable scales of variation, not to predict values for particular places. (The log-linear model can also be interpreted in terms of percent changes in the covariate, such as ’a 10% increase in GDP per capita, e.g., from 10,000 to 11,000, is associated with around 4 * log(11/10) = 0.38 additional years of life expectancy on average.)
# scale GDP
preds <- transform(preds, gdpPercap = gdpPercap / 1e3)
-yrange <- c(min(preds$`2.5%`), max(preds$`97.5%`))
+## yrange <- c(min(preds$`2.5%`), max(preds$`97.5%`))
+yrange <- c(57, 85)
+
plot(preds$gdpPercap, preds$mean,
t = 'l',
ylim = yrange,
@@ -541,7 +543,10 @@ Predicted values
# add credible intervals
lines(preds$gdpPercap, preds$`2.5%`, lty = 3)
-lines(preds$gdpPercap, preds$`97.5%`, lty = 3)
Per this dataset, about 50% of the world population lives in countries with GDP per capita below $12,300.
diff --git a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-21-1.png b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-21-1.png index 0348dd66..fbee081f 100644 Binary files a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-21-1.png and b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-21-1.png differ diff --git a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-23-1.png index 39cae0ec..2f43d60f 100644 Binary files a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-27-1.png b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-27-1.png index c69e2495..675a4e34 100644 Binary files a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-27-1.png and b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-27-1.png differ diff --git a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-32-1.png b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-32-1.png index b961248f..7fbbcfad 100644 Binary files a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-32-1.png and b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-32-1.png differ diff --git a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-37-1.png b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-37-1.png index a814b262..b6f8c4d0 100644 Binary files a/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-37-1.png and b/docs/articles/spatial-analysis_files/figure-html/unnamed-chunk-37-1.png differ diff --git a/docs/articles/spatial-me-models.html b/docs/articles/spatial-me-models.html index 1279cc92..ca1251ca 100644 --- a/docs/articles/spatial-me-models.html +++ b/docs/articles/spatial-me-models.html @@ -40,7 +40,7 @@ @@ -218,23 +218,17 @@## Warning: The largest R-hat is 1.06, indicating chains have not mixed.
-## Running the chains for more iterations may help. See
-## https://mc-stan.org/misc/warnings.html#r-hat
-## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
-## Running the chains for more iterations may help. See
-## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
-+print(fit)
## Spatial Model Results ## Formula: deaths.male ~ offset(log(pop.at.risk.male)) + insurance ## Likelihood: poisson ## Link: log ## Spatial method: CAR -## Residual Moran Coefficient: 0.001679231 +## Residual Moran Coefficient: 0.0008961538 ## Observations: 159 ## Data models (ME): insurance ## Data model (ME prior): CAR (auto-normal) @@ -243,18 +237,18 @@
Spatial ME model## post-warmup draws per chain=325, total post-warmup draws=1300. ## ## mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat -## intercept -4.134 0.029 0.177 -4.276 -4.202 -4.162 -4.119 -3.574 38 1.113 -## insurance -2.736 0.024 0.781 -4.280 -3.377 -2.726 -2.090 -1.226 1023 1.000 -## car_rho 0.871 0.005 0.093 0.640 0.806 0.888 0.948 0.998 312 1.016 -## car_scale 0.447 0.001 0.034 0.384 0.418 0.445 0.476 0.517 817 1.006 +## intercept -4.167 0.003 0.063 -4.283 -4.205 -4.164 -4.127 -4.033 367 1.006 +## insurance -2.762 0.025 0.732 -4.170 -3.404 -2.759 -2.137 -1.344 835 1.001 +## car_rho 0.865 0.003 0.091 0.639 0.793 0.883 0.943 0.990 938 1.004 +## car_scale 0.448 0.001 0.033 0.388 0.419 0.446 0.476 0.516 726 1.009 ## -## Samples were drawn using NUTS(diag_e) at Sat Nov 16 10:47:47 2024. +## Samples were drawn using NUTS(diag_e) at Wed Dec 4 13:09:02 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).
The CAR models in geostan can be highly efficient in terms of MCMC sampling, but the required number of iterations will depend on many characteristics of the model and data. Often the default
iter = 2000
is more than sufficient (withcores = 4
). To speed up sampling with multi-core processing, usecores = 4
(to sample 4 chains in parallel).In the following section, methods for examining MCMC samples of the modeled covariate values \(\boldsymbol x\) will be illustrated. Note that the process of storing MCMC samples for \(\boldsymbol x\) can become computationally burdensome if you have multiple covariates and moderately large N. If you do not need the samples of \(\boldsymbol x\) you can use the
-drop
argument, as in:+fit2 <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + insurance, centerx = TRUE, ME = ME_list, @@ -272,7 +266,7 @@
Visual diagnostics\(z_i\) an the modeled values \(x_i\): \[\delta_i = z_i - x_i.\] We want to look for any strong spatial pattern in these \(\delta_i\) values, because that would be an indication of a bias. However, the magnitude of the \(\delta_i\) value is important to consider—there may be a pattern, but if the amount of shrinkage is very small, that pattern may not matter. (The model returns \(M\) samples from the posterior distribution of each \(x_i\); or, indexing by MCMC sample \(x_i^m\) (\(m\) is an index, not an exponent). The reported \(\delta_i\) values is the MCMC mean \(\delta_i = \frac{1}{M} \sum_m z_i - x_i^m\).)
Two figures are provided to evaluate patterns in \(\delta_i\): first, a Moran scatter plot and, second, a map.
-@@ -132,7 +132,7 @@+me_diag(fit, 'insurance', georgia)
In this case, the results do not look too concerning insofar as there are no conspicuous patterns. However, a number of the credible intervals on the modeled values are large, which indicates low data reliability. The fact that some of the \(\delta_i\) values are substantial also points to low data quality for those estimates.
@@ -282,54 +276,54 @@Working with MCMC samples from
geostan consists of pre-compiled Stan models, and users can always access the Markov chain Monte Carlo (MCMC) samples returned by Stan. When extracted as a matrix of samples (as below), each row of the matrix represents a draw from the joint probability distribution for all model parameters, and each column consists of samples from the marginal distribution of each parameter.
The ME models return samples for \(x_i\) as well as the model parameters \(\mu\) (“mu_x_true”), \(\rho\) (“car_rho_x_true”), and \(\tau\) (“sigma_x_true”). We can access these using
-as.matrix
(alsoas.array
andas.data.frame
).diff --git a/docs/reference/index.html b/docs/reference/index.html index 71b8cb24..0b68e88a 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -17,7 +17,7 @@+-## [1] 1300 1
diff --git a/docs/reference/impacts-1.png b/docs/reference/impacts-1.png index 50645e60..0bd11987 100644 Binary files a/docs/reference/impacts-1.png and b/docs/reference/impacts-1.png differ diff --git a/docs/reference/impacts.html b/docs/reference/impacts.html index d4d74120..77e52918 100644 --- a/docs/reference/impacts.html +++ b/docs/reference/impacts.html @@ -17,7 +17,7 @@+head(mu.x)
-## parameters ## iterations mu_x_true[1] -## [1,] 1.878425 -## [2,] 1.764369 -## [3,] 1.811935 -## [4,] 1.787079 -## [5,] 1.812591 -## [6,] 1.767798
diff --git a/docs/reference/gr.html b/docs/reference/gr.html index ee0a60b1..092b213a 100644 --- a/docs/reference/gr.html +++ b/docs/reference/gr.html @@ -17,7 +17,7 @@+## [1,] 1.710997 +## [2,] 1.823856 +## [3,] 1.721450 +## [4,] 1.844658 +## [5,] 1.719539 +## [6,] 1.923031+-mean(mu.x)
+## [1] 1.784293
## [1] 1.786588
We can visualize these using
-plot
or print a summary:diff --git a/docs/reference/get_shp.html b/docs/reference/get_shp.html index 7302f31b..b40e534c 100644 --- a/docs/reference/get_shp.html +++ b/docs/reference/get_shp.html @@ -17,7 +17,7 @@+## Inference for Stan model: foundation. ## 4 chains, each with iter=650; warmup=325; thin=1; ## post-warmup draws per chain=325, total post-warmup draws=1300. ## ## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat -## mu_x_true[1] 1.78 0 0.08 1.62 1.74 1.79 1.83 1.95 521 1.01 -## car_rho_x_true[1] 0.91 0 0.06 0.76 0.88 0.93 0.96 0.99 1039 1.00 -## sigma_x_true[1] 0.48 0 0.03 0.42 0.46 0.48 0.50 0.55 1406 1.00 +## mu_x_true[1] 1.79 0 0.08 1.64 1.74 1.79 1.83 1.94 576 1 +## car_rho_x_true[1] 0.91 0 0.07 0.73 0.87 0.92 0.96 0.99 816 1 +## sigma_x_true[1] 0.48 0 0.04 0.42 0.46 0.48 0.51 0.56 933 1 ## -## Samples were drawn using NUTS(diag_e) at Sat Nov 16 10:47:47 2024. +## Samples were drawn using NUTS(diag_e) at Wed Dec 4 13:09:02 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).
To extract samples from the joint probability distribution for \(\boldsymbol x\), use the generic parameter name “x_true”:
-diff --git a/docs/reference/geostan-package.html b/docs/reference/geostan-package.html index 083c1274..50c75b9b 100644 --- a/docs/reference/geostan-package.html +++ b/docs/reference/geostan-package.html @@ -17,7 +17,7 @@+## [1] 1300 159
If we wanted to calculate the mean of each of these marginal distributions (one for every \(x_i\)), we could use
-apply
withMARGIN = 2
to summarize by column:diff --git a/docs/reference/georgia-1.png b/docs/reference/georgia-1.png new file mode 100644 index 00000000..0ff911ee Binary files /dev/null and b/docs/reference/georgia-1.png differ diff --git a/docs/reference/georgia.html b/docs/reference/georgia.html index 3f2e293c..caafa3c3 100644 --- a/docs/reference/georgia.html +++ b/docs/reference/georgia.html @@ -17,7 +17,7 @@++## 0.8702400## x_insurance[1] x_insurance[2] x_insurance[3] x_insurance[4] x_insurance[5] -## 0.8367041 0.8001162 0.8516995 0.8520539 0.9179434 +## 0.8368129 0.8006499 0.8518878 0.8520715 0.9177817 ## x_insurance[6] -## 0.8702477
The vector
-x.mu
contains estimates (posterior means) for \(x_i\). We might want to use these to plot the residuals or fitted values against the predictor:diff --git a/docs/reference/expected_mc.html b/docs/reference/expected_mc.html index 7a46b3cc..048e6a2b 100644 --- a/docs/reference/expected_mc.html +++ b/docs/reference/expected_mc.html @@ -17,7 +17,7 @@+diff --git a/docs/reference/eigen_grid.html b/docs/reference/eigen_grid.html index cc0c0b80..fc725877 100644 --- a/docs/reference/eigen_grid.html +++ b/docs/reference/eigen_grid.html @@ -17,7 +17,7 @@res <- resid(fit)$mean plot(x.mu, res, xlab = 'Insurance rate', @@ -347,7 +341,7 @@
Working with MCMC samples from
Non-spatial ME models
If the
-ME
list doesn’t have a slot withcar_parts
, geostan will automatically use a non-spatial Student’s t model instead of the CAR model: \[ p( x | \mathcal{M}) = Student(\boldsymbol x | \nu, \mu, \sigma), \] with degrees of freedom \(\nu\), mean \(\mu\), and scale \(\sigma\).diff --git a/docs/reference/edges-1.png b/docs/reference/edges-1.png new file mode 100644 index 00000000..24e29143 Binary files /dev/null and b/docs/reference/edges-1.png differ diff --git a/docs/reference/edges.html b/docs/reference/edges.html index 836cb6aa..1c5ef262 100644 --- a/docs/reference/edges.html +++ b/docs/reference/edges.html @@ -17,7 +17,7 @@+@@ -52,7 +52,8 @@ME_nsp <- prep_me_data( se = data.frame(insurance = georgia$insurance.se), logit = TRUE @@ -358,7 +352,7 @@
Non-spatial ME modelsMultiple covariates
To model multiple covariates, simply add them to the data frame of standard errors:
-diff --git a/docs/news/index.html b/docs/news/index.html index 10a40af9..9bfec058 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -17,7 +17,7 @@+diff --git a/docs/index.html b/docs/index.html index 34cb49ad..477a99f5 100644 --- a/docs/index.html +++ b/docs/index.html @@ -40,7 +40,7 @@georgia$college <- georgia$college / 100 georgia$college.se <- georgia$college.se / 100 @@ -381,7 +375,7 @@
Log transforms
Income and other magnitudes are often log transformed. In that case, the survey standard errors also need to be transformed. The
se_log
function provide approximate standard errors for this purpose.When using
-se_log
, pass in the variable itself (untransformed) and its standard errors. Here is a full workflow for a spatial mortality model with covariate measurement error:diff --git a/docs/authors.html b/docs/authors.html index 4b1e7983..1cefb671 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -17,7 +17,7 @@+@@ -97,7 +97,7 @@data(georgia) # income in $1,000s diff --git a/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-14-1.png index 55ff32c1..20f9df08 100644 Binary files a/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-14-1.png and b/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-9-1.png index 75df6bbb..1d1c61fe 100644 Binary files a/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-9-1.png and b/docs/articles/spatial-me-models_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/spatial-weights-matrix.html b/docs/articles/spatial-weights-matrix.html index 428c08ff..ed9732ec 100644 --- a/docs/articles/spatial-weights-matrix.html +++ b/docs/articles/spatial-weights-matrix.html @@ -40,7 +40,7 @@
Getting startedTo get started, load the geostan and sf packages:
Changelog
geostan 0.8.1
-This version changes one line in the ‘Spatial analysis with geostan’ vignette that caused the package to fail on installation for old versions of windows and mac os.
+This is a small patch to ensure that the package can install on older versions of R.
+The patch changes one line in the ‘Spatial analysis with geostan’ vignette that caused the package to fail on installation for old versions of windows and mac os.
diff --git a/docs/reference/auto_gaussian.html b/docs/reference/auto_gaussian.html index b3012929..f573ce00 100644 --- a/docs/reference/auto_gaussian.html +++ b/docs/reference/auto_gaussian.html @@ -17,7 +17,7 @@geostan 0.8.02024-11-16
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 82c345ee..464d57e5 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -8,7 +8,7 @@ articles: spatial-analysis: spatial-analysis.html spatial-me-models: spatial-me-models.html spatial-weights-matrix: spatial-weights-matrix.html -last_built: 2024-11-16T16:47Z +last_built: 2024-12-04T21:46Z urls: reference: https://connordonegan.github.io/geostan/reference article: https://connordonegan.github.io/geostan/articles diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png index 03f31164..8ee18a05 100644 Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index cdd59c35..c7873260 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/Rplot003.png b/docs/reference/Rplot003.png index ba9c0f30..68a64799 100644 Binary files a/docs/reference/Rplot003.png and b/docs/reference/Rplot003.png differ diff --git a/docs/reference/Rplot004.png b/docs/reference/Rplot004.png index 32fbe136..be01a383 100644 Binary files a/docs/reference/Rplot004.png and b/docs/reference/Rplot004.png differ diff --git a/docs/reference/Rplot005.png b/docs/reference/Rplot005.png index f3b1a70b..a9071896 100644 Binary files a/docs/reference/Rplot005.png and b/docs/reference/Rplot005.png differ diff --git a/docs/reference/Rplot006.png b/docs/reference/Rplot006.png new file mode 100644 index 00000000..0afe040a Binary files /dev/null and b/docs/reference/Rplot006.png differ diff --git a/docs/reference/aple.html b/docs/reference/aple.html index 9ea8ca63..17fe07e6 100644 --- a/docs/reference/aple.html +++ b/docs/reference/aple.html @@ -17,7 +17,7 @@Model prep
- + Extract eigenfunctions of a connectivity matrix for spatial filtering
Prepare data for spatial filtering
diff --git a/docs/reference/lg-1.png b/docs/reference/lg-1.png new file mode 100644 index 00000000..e842b509 Binary files /dev/null and b/docs/reference/lg-1.png differ diff --git a/docs/reference/lg-2.png b/docs/reference/lg-2.png new file mode 100644 index 00000000..5faa9527 Binary files /dev/null and b/docs/reference/lg-2.png differ diff --git a/docs/reference/lg.html b/docs/reference/lg.html index edb17874..9d5942be 100644 --- a/docs/reference/lg.html +++ b/docs/reference/lg.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/lisa-1.png b/docs/reference/lisa-1.png new file mode 100644 index 00000000..7ce8272e Binary files /dev/null and b/docs/reference/lisa-1.png differ diff --git a/docs/reference/lisa.html b/docs/reference/lisa.html index 19b622bf..ff8759e5 100644 --- a/docs/reference/lisa.html +++ b/docs/reference/lisa.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/log_lik.html b/docs/reference/log_lik.html index 90ea8989..d4433dfc 100644 --- a/docs/reference/log_lik.html +++ b/docs/reference/log_lik.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/make_EV-1.png b/docs/reference/make_EV-1.png new file mode 100644 index 00000000..9e03ae6f Binary files /dev/null and b/docs/reference/make_EV-1.png differ diff --git a/docs/reference/make_EV.html b/docs/reference/make_EV.html index 4bfe3091..84752db7 100644 --- a/docs/reference/make_EV.html +++ b/docs/reference/make_EV.html @@ -1,5 +1,5 @@ - Extract eigenfunctions of a connectivity matrix for spatial filtering — make_EV • geostan Prepare data for spatial filtering — make_EV • geostan @@ -17,7 +17,7 @@ @@ -46,13 +46,13 @@diff --git a/docs/reference/n_eff-1.png b/docs/reference/n_eff-1.png new file mode 100644 index 00000000..671a40fe Binary files /dev/null and b/docs/reference/n_eff-1.png differ diff --git a/docs/reference/n_eff.html b/docs/reference/n_eff.html index d679b34d..3a1f5e07 100644 --- a/docs/reference/n_eff.html +++ b/docs/reference/n_eff.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/n_nbs.html b/docs/reference/n_nbs.html index d53ea1cd..6def4775 100644 --- a/docs/reference/n_nbs.html +++ b/docs/reference/n_nbs.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/posterior_predict-1.png b/docs/reference/posterior_predict-1.png new file mode 100644 index 00000000..398c5c88 Binary files /dev/null and b/docs/reference/posterior_predict-1.png differ diff --git a/docs/reference/posterior_predict.html b/docs/reference/posterior_predict.html index 38e24847..df690c49 100644 --- a/docs/reference/posterior_predict.html +++ b/docs/reference/posterior_predict.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/predict.geostan_fit-1.png b/docs/reference/predict.geostan_fit-1.png new file mode 100644 index 00000000..dd15f8a8 Binary files /dev/null and b/docs/reference/predict.geostan_fit-1.png differ diff --git a/docs/reference/predict.geostan_fit-2.png b/docs/reference/predict.geostan_fit-2.png new file mode 100644 index 00000000..0b81a0ee Binary files /dev/null and b/docs/reference/predict.geostan_fit-2.png differ diff --git a/docs/reference/predict.geostan_fit.html b/docs/reference/predict.geostan_fit.html index 6d842996..037cc84a 100644 --- a/docs/reference/predict.geostan_fit.html +++ b/docs/reference/predict.geostan_fit.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_car_data.html b/docs/reference/prep_car_data.html index 612315ea..564b2b9a 100644 --- a/docs/reference/prep_car_data.html +++ b/docs/reference/prep_car_data.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_car_data2.html b/docs/reference/prep_car_data2.html index babd9f89..a492bc21 100644 --- a/docs/reference/prep_car_data2.html +++ b/docs/reference/prep_car_data2.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_icar_data.html b/docs/reference/prep_icar_data.html index 7a7b905b..e54357a4 100644 --- a/docs/reference/prep_icar_data.html +++ b/docs/reference/prep_icar_data.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_me_data.html b/docs/reference/prep_me_data.html index 7c16d052..a626ff92 100644 --- a/docs/reference/prep_me_data.html +++ b/docs/reference/prep_me_data.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_sar_data.html b/docs/reference/prep_sar_data.html index 91c2771a..a437c9f2 100644 --- a/docs/reference/prep_sar_data.html +++ b/docs/reference/prep_sar_data.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/prep_sar_data2.html b/docs/reference/prep_sar_data2.html index 82199c8a..ef80c108 100644 --- a/docs/reference/prep_sar_data2.html +++ b/docs/reference/prep_sar_data2.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/print_geostan_fit-1.png b/docs/reference/print_geostan_fit-1.png new file mode 100644 index 00000000..57460a49 Binary files /dev/null and b/docs/reference/print_geostan_fit-1.png differ diff --git a/docs/reference/print_geostan_fit.html b/docs/reference/print_geostan_fit.html index 313a82c3..399ba685 100644 --- a/docs/reference/print_geostan_fit.html +++ b/docs/reference/print_geostan_fit.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/priors-1.png b/docs/reference/priors-1.png new file mode 100644 index 00000000..e0763d25 Binary files /dev/null and b/docs/reference/priors-1.png differ diff --git a/docs/reference/priors.html b/docs/reference/priors.html index dbd64bdc..debf1f58 100644 --- a/docs/reference/priors.html +++ b/docs/reference/priors.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/resid_geostan_fit-1.png b/docs/reference/resid_geostan_fit-1.png new file mode 100644 index 00000000..12995dcb Binary files /dev/null and b/docs/reference/resid_geostan_fit-1.png differ diff --git a/docs/reference/resid_geostan_fit.html b/docs/reference/resid_geostan_fit.html index c2dceae1..610e540d 100644 --- a/docs/reference/resid_geostan_fit.html +++ b/docs/reference/resid_geostan_fit.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/row_standardize.html b/docs/reference/row_standardize.html index bf406c6b..cca73ef5 100644 --- a/docs/reference/row_standardize.html +++ b/docs/reference/row_standardize.html @@ -17,7 +17,7 @@ @@ -84,11 +84,11 @@diff --git a/docs/reference/moran_plot-1.png b/docs/reference/moran_plot-1.png new file mode 100644 index 00000000..7203c4dc Binary files /dev/null and b/docs/reference/moran_plot-1.png differ diff --git a/docs/reference/moran_plot.html b/docs/reference/moran_plot.html index 5025dd69..6f03862b 100644 --- a/docs/reference/moran_plot.html +++ b/docs/reference/moran_plot.html @@ -17,7 +17,7 @@-Extract eigenfunctions of a connectivity matrix for spatial filtering
+Prepare data for spatial filtering
Source:R/convenience-functions.R
make_EV.Rd
-Extract eigenfunctions of a connectivity matrix for spatial filtering
+Prepare data for spatial filtering
diff --git a/docs/reference/mc.html b/docs/reference/mc.html index 8a26f891..4f743582 100644 --- a/docs/reference/mc.html +++ b/docs/reference/mc.html @@ -17,7 +17,7 @@diff --git a/docs/reference/me_diag-1.png b/docs/reference/me_diag-1.png new file mode 100644 index 00000000..c9f70b72 Binary files /dev/null and b/docs/reference/me_diag-1.png differ diff --git a/docs/reference/me_diag-2.png b/docs/reference/me_diag-2.png new file mode 100644 index 00000000..9330d48a Binary files /dev/null and b/docs/reference/me_diag-2.png differ diff --git a/docs/reference/me_diag.html b/docs/reference/me_diag.html index 5db16b3f..1b558577 100644 --- a/docs/reference/me_diag.html +++ b/docs/reference/me_diag.html @@ -17,7 +17,7 @@Value
Examples
diff --git a/docs/reference/samples_geostan_fit.html b/docs/reference/samples_geostan_fit.html index 32392e16..3de50fdc 100644 --- a/docs/reference/samples_geostan_fit.html +++ b/docs/reference/samples_geostan_fit.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/se_log-1.png b/docs/reference/se_log-1.png new file mode 100644 index 00000000..97b8f584 Binary files /dev/null and b/docs/reference/se_log-1.png differ diff --git a/docs/reference/se_log.html b/docs/reference/se_log.html index 9c0d49ba..7d65e0f9 100644 --- a/docs/reference/se_log.html +++ b/docs/reference/se_log.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/sentencing.html b/docs/reference/sentencing.html index ed1bd965..bb5e22af 100644 --- a/docs/reference/sentencing.html +++ b/docs/reference/sentencing.html @@ -18,7 +18,7 @@ diff --git a/docs/reference/shape2mat.html b/docs/reference/shape2mat.html index 7bbb3d69..bc97e8d8 100644 --- a/docs/reference/shape2mat.html +++ b/docs/reference/shape2mat.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/sim_sar-1.png b/docs/reference/sim_sar-1.png new file mode 100644 index 00000000..47868570 Binary files /dev/null and b/docs/reference/sim_sar-1.png differ diff --git a/docs/reference/sim_sar-2.png b/docs/reference/sim_sar-2.png new file mode 100644 index 00000000..1f49f57a Binary files /dev/null and b/docs/reference/sim_sar-2.png differ diff --git a/docs/reference/sim_sar.html b/docs/reference/sim_sar.html index eafad13f..84871692 100644 --- a/docs/reference/sim_sar.html +++ b/docs/reference/sim_sar.html @@ -17,7 +17,7 @@ @@ -64,7 +64,7 @@A <- shape2mat(georgia) head(Matrix::summary(A)) -Matrix::rowSums(A) +Matrix::rowSums(A) W <- row_standardize(A) head(Matrix::summary(W)) -Matrix::rowSums(W) +Matrix::rowSums(W)
Simulate spatially autocorrelated data
w, type = c("SEM", "SLM"), approx = FALSE, - K = 15, + K = 20, ... ) diff --git a/docs/reference/sp_diag-1.png b/docs/reference/sp_diag-1.png new file mode 100644 index 00000000..fcf10b42 Binary files /dev/null and b/docs/reference/sp_diag-1.png differ diff --git a/docs/reference/sp_diag-2.png b/docs/reference/sp_diag-2.png new file mode 100644 index 00000000..d95b82b2 Binary files /dev/null and b/docs/reference/sp_diag-2.png differ diff --git a/docs/reference/sp_diag-3.png b/docs/reference/sp_diag-3.png new file mode 100644 index 00000000..656caf06 Binary files /dev/null and b/docs/reference/sp_diag-3.png differ diff --git a/docs/reference/sp_diag.html b/docs/reference/sp_diag.html index 8667c39f..d164bdec 100644 --- a/docs/reference/sp_diag.html +++ b/docs/reference/sp_diag.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/stan_car-1.png b/docs/reference/stan_car-1.png index 8067b752..497a7b97 100644 Binary files a/docs/reference/stan_car-1.png and b/docs/reference/stan_car-1.png differ diff --git a/docs/reference/stan_car-2.png b/docs/reference/stan_car-2.png index 8f42b067..4f179a26 100644 Binary files a/docs/reference/stan_car-2.png and b/docs/reference/stan_car-2.png differ diff --git a/docs/reference/stan_car-3.png b/docs/reference/stan_car-3.png index 5a2519a4..44f02b6c 100644 Binary files a/docs/reference/stan_car-3.png and b/docs/reference/stan_car-3.png differ diff --git a/docs/reference/stan_car.html b/docs/reference/stan_car.html index f03fd495..6cea3e6a 100644 --- a/docs/reference/stan_car.html +++ b/docs/reference/stan_car.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/stan_esf-1.png b/docs/reference/stan_esf-1.png new file mode 100644 index 00000000..bf787099 Binary files /dev/null and b/docs/reference/stan_esf-1.png differ diff --git a/docs/reference/stan_esf-2.png b/docs/reference/stan_esf-2.png new file mode 100644 index 00000000..0ceaf158 Binary files /dev/null and b/docs/reference/stan_esf-2.png differ diff --git a/docs/reference/stan_esf-3.png b/docs/reference/stan_esf-3.png new file mode 100644 index 00000000..34f3fdc2 Binary files /dev/null and b/docs/reference/stan_esf-3.png differ diff --git a/docs/reference/stan_esf.html b/docs/reference/stan_esf.html index b7e472ca..6152417a 100644 --- a/docs/reference/stan_esf.html +++ b/docs/reference/stan_esf.html @@ -17,7 +17,7 @@ @@ -289,25 +289,19 @@Value
diff --git a/docs/reference/stan_glm-1.png b/docs/reference/stan_glm-1.png index 94cc578b..cea05628 100644 Binary files a/docs/reference/stan_glm-1.png and b/docs/reference/stan_glm-1.png differ diff --git a/docs/reference/stan_glm-2.png b/docs/reference/stan_glm-2.png index ccf1d05f..40821e30 100644 Binary files a/docs/reference/stan_glm-2.png and b/docs/reference/stan_glm-2.png differ diff --git a/docs/reference/stan_glm-3.png b/docs/reference/stan_glm-3.png index 3a3778cb..dafcbbac 100644 Binary files a/docs/reference/stan_glm-3.png and b/docs/reference/stan_glm-3.png differ diff --git a/docs/reference/stan_glm-4.png b/docs/reference/stan_glm-4.png index 8da1b5e2..777cf23b 100644 Binary files a/docs/reference/stan_glm-4.png and b/docs/reference/stan_glm-4.png differ diff --git a/docs/reference/stan_glm-5.png b/docs/reference/stan_glm-5.png index f280d62a..383c2aba 100644 Binary files a/docs/reference/stan_glm-5.png and b/docs/reference/stan_glm-5.png differ diff --git a/docs/reference/stan_glm-6.png b/docs/reference/stan_glm-6.png new file mode 100644 index 00000000..66f52bde Binary files /dev/null and b/docs/reference/stan_glm-6.png differ diff --git a/docs/reference/stan_glm.html b/docs/reference/stan_glm.html index 3efe0ed5..0d999655 100644 --- a/docs/reference/stan_glm.html +++ b/docs/reference/stan_glm.html @@ -17,7 +17,7 @@ @@ -254,24 +254,6 @@Details
Eigenvector spatial filtering (ESF) is a method for spatial regression analysis. ESF is extensively covered in Griffith et al. (2019). This function implements the methodology introduced in Donegan et al. (2020), which uses Piironen and Vehtari's (2017) regularized horseshoe prior.
-ESF decomposes spatial autocorrelation into a linear combination of various patterns, typically at different scales (such as local, regional, and global trends). By adding a spatial filter to a regression model, these spatial autocorrelation patterns are shifted from the residuals to the spatial filter. ESF models take the spectral decomposition of a transformed spatial connectivity matrix, \(C\). The resulting eigenvectors, \(E\), are mutually orthogonal and uncorrelated map patterns. The spatial filter equals \(E \beta_{E}\) where \(\beta_{E}\) is a vector of coefficients.
+By adding a spatial filter to a regression model, spatial autocorrelation patterns are shifted from the residuals to the spatial filter. ESF models take the spectral decomposition of a transformed spatial connectivity matrix, \(C\). The resulting eigenvectors, \(E\), are mutually orthogonal and uncorrelated map patterns (at various scales, 'local' to 'regional' to 'global'). The spatial filter equals \(E \beta_{E}\) where \(\beta_{E}\) is a vector of coefficients.
ESF decomposes the data into a global mean, \(\alpha\), global patterns contributed by covariates \(X \beta\), spatial trends \(E \beta_{E}\), and residual variation. Thus, for
family=gaussian()
,$$ y \sim Gauss(\alpha + X * \beta + E \beta_{E}, \sigma). $$
-An ESF component can be incorporated into the linear predictor of any generalized linear model. For example, a spatial Poisson model for rare disease incidence may be specified as follows:
-$$y \sim Poisson(e^{O + \mu})$$ -$$\mu = \alpha + E \beta_{E} + A $$ -$$ A \sim Guass(0, \tau) $$ -$$ \tau \sim student(20, 0, 2) $$ -$$ \beta_{E} \sim horseshoe(.) $$
-The form of this model is similar to the BYM model (see stan_icar), in the sense that it contains a spatially structured trend term (\(E \beta_{E}\)) and an unstructured ('random effects') term (\(A\)).
+An ESF component can be incorporated into the linear predictor of any generalized linear model. For example, using
stan_esf
withfamily = poisson()
and adding a 'random effects' term for each spatial unit (via there
argument) will produce a model that resembles the BYM model (combining spatially structured and spatially-unstructured components).The spatial.geostan_fit method will return \(E \beta_{E}\).
The model can also be extended to the space-time domain; see shape2mat to specify a space-time connectivity matrix.
The coefficients \(\beta_{E}\) are assigned the regularized horseshoe prior (Piironen and Vehtari, 2017), resulting in a relatively sparse model specification. In addition, numerous eigenvectors are automatically dropped because they represent trace amounts of spatial autocorrelation (this is controlled by the
threshold
argument). By default,stan_esf
will drop all eigenvectors representing negative spatial autocorrelation patterns. You can change this behavior using thensa
argument.Additional functionality
-The CAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for stan_glm.
+The ESF models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for stan_glm.
Value
@@ -396,30 +396,6 @@Details
Fit a generalized linear model using the R formula interface. Default prior distributions are designed to be weakly informative relative to the data. Much of the functionality intended for spatial models, such as the ability to add spatially lagged covariates and observational error models, are also available in
stan_glm
. All ofgeostan
's spatial models build on top of the same Stan code used instan_glm
.-- -Poisson models and disease mapping
- - -In spatial statistics, Poisson models are often used to calculate incidence rates (mortality rates, or disease incidence rates) for administrative areas like counties or census tracts. If \(y\) are counts of cases, and \(P\) are populations at risk, then the crude rates are \(y/P\). The purpose is to model risk \(\eta\) for which crude rates are a (noisy) indicator. Our analysis should also respect the fact that the amount of information contained in the observations \(y/P\) increases with \(P\). Hierarchical Poisson models are often used to incorporate all of this information.
-For the Poisson model, \(y\) is specified as the outcome and the log of the population at risk
-log(P)
needs to be provided as an offset term. For such a case, disease incidence across the collection of areas could be modeled as:$$y \sim Poisson(e^{log(P) + \eta})$$ -$$ \eta = \alpha + A$$ -$$ A \sim Gauss(0, \tau)$$ -$$\tau \sim Student(20, 0, 2)$$ -where \(\alpha\) is the mean log-risk (incidence rate) and \(A\) is a vector of (so-called) random effects, which enable partial pooling of information across observations. Covariates can be added to the model for the log-rates, such that \(\eta = \alpha + X \beta + A\).
-Note that the denominator for the rates is specified as a log-offset to provide a consistent, formula-line interface to the model. Using the log-offest (as above) is equivalent to the following: -$$ -y \sim Poisson(P * e^{\eta}) -$$ -where \(P\) is still the population at risk and it is multiplied by \(e^{\eta}\), the incidence rate (risk).
-@@ -363,7 +363,7 @@Spatially lagged covariates (SLX)
@@ -381,6 +363,11 @@Examples
print(fit2) +# example prior for two covariates +pl <- list(beta = normal(c(0, 0), + c(1, 1)) + ) + ## ## Poisson model for count data ## with county 'random effects' @@ -409,15 +396,18 @@Examples
# summary of results with MCMC diagnostics print(fit.pois) +# \donttest{ # MCMC diagnostics plot: Rhat values should all by very near 1 rstan::stan_rhat(fit.pois$stanfit) + # effective sample size for all parameters and generated quantities # (including residuals, predicted values, etc.) rstan::stan_ess(fit.pois$stanfit) # or for a particular parameter rstan::stan_ess(fit.pois$stanfit, "alpha_re") +# } ## ## Visualize the posterior predictive distribution diff --git a/docs/reference/stan_icar-1.png b/docs/reference/stan_icar-1.png new file mode 100644 index 00000000..9f9d3300 Binary files /dev/null and b/docs/reference/stan_icar-1.png differ diff --git a/docs/reference/stan_icar-2.png b/docs/reference/stan_icar-2.png new file mode 100644 index 00000000..5a061b18 Binary files /dev/null and b/docs/reference/stan_icar-2.png differ diff --git a/docs/reference/stan_icar-3.png b/docs/reference/stan_icar-3.png new file mode 100644 index 00000000..6b86428e Binary files /dev/null and b/docs/reference/stan_icar-3.png differ diff --git a/docs/reference/stan_icar.html b/docs/reference/stan_icar.html index cb81eb4c..8f4e96a5 100644 --- a/docs/reference/stan_icar.html +++ b/docs/reference/stan_icar.html @@ -17,7 +17,7 @@BYM2
Additional functionality
-The CAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for stan_glm.
+The ICAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for stan_glm.
Examples
library(rstan) rstan::stan_ess(fit.bym$stanfit) rstan::stan_rhat(fit.bym$stanfit) - -# calculate log-standardized incidence ratios -# (observed/exected case counts) -library(ggplot2) -library(sf) - -f <- fitted(fit.bym, rates = FALSE)$mean -SSR <- f / sentencing$expected_sents -log.SSR <- log( SSR, base = 2) - -ggplot(sentencing) + - geom_sf(aes(fill = log.SSR)) + - scale_fill_gradient2( - low = "navy", - high = "darkred" - ) + - labs(title = "Log-standardized incidence ratios", - subtitle = "log( Fitted/Expected), base 2") + - theme_void() + - theme( - legend.position = "bottom", - legend.key.height = unit(0.35, "cm"), - legend.key.width = unit(1.5, "cm") - ) # } diff --git a/docs/reference/stan_sar-1.png b/docs/reference/stan_sar-1.png index fc9ab771..7e797922 100644 Binary files a/docs/reference/stan_sar-1.png and b/docs/reference/stan_sar-1.png differ diff --git a/docs/reference/stan_sar-2.png b/docs/reference/stan_sar-2.png index abbca300..702280d3 100644 Binary files a/docs/reference/stan_sar-2.png and b/docs/reference/stan_sar-2.png differ diff --git a/docs/reference/stan_sar.html b/docs/reference/stan_sar.html index 74ea5c74..1a721622 100644 --- a/docs/reference/stan_sar.html +++ b/docs/reference/stan_sar.html @@ -17,7 +17,7 @@ diff --git a/docs/reference/waic.html b/docs/reference/waic.html index 8a5d512d..bc16c6e5 100644 --- a/docs/reference/waic.html +++ b/docs/reference/waic.html @@ -17,7 +17,7 @@ diff --git a/man/make_EV.Rd b/man/make_EV.Rd index 832d71b7..7739249b 100644 --- a/man/make_EV.Rd +++ b/man/make_EV.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/convenience-functions.R \name{make_EV} \alias{make_EV} -\title{Extract eigenfunctions of a connectivity matrix for spatial filtering} +\title{Prepare data for spatial filtering} \source{ Daniel Griffith and Yongwan Chun. 2014. "Spatial Autocorrelation and Spatial Filtering." in M. M. Fischer and P. Nijkamp (eds.), \emph{Handbook of Regional Science.} Springer. } @@ -23,7 +23,7 @@ Setting \code{nsa = TRUE} will result in a candidate set of EVs that contains ei A \code{data.frame} of eigenvectors for spatial filtering. If \code{values=TRUE} then a named list is returned with elements \code{eigenvectors} and \code{eigenvalues}. } \description{ -Extract eigenfunctions of a connectivity matrix for spatial filtering +Prepare data for spatial filtering } \details{ Returns a set of eigenvectors related to the Moran coefficient (MC), limited to those eigenvectors with |MC| > \code{threshold} if \code{nsa = TRUE} or MC > \code{threshold} if \code{nsa = FALSE}, optionally with corresponding eigenvalues. diff --git a/man/stan_esf.Rd b/man/stan_esf.Rd index 30582e04..9f28a14e 100644 --- a/man/stan_esf.Rd +++ b/man/stan_esf.Rd @@ -147,7 +147,7 @@ Fit a spatial regression model using eigenvector spatial filtering (ESF). \details{ Eigenvector spatial filtering (ESF) is a method for spatial regression analysis. ESF is extensively covered in Griffith et al. (2019). This function implements the methodology introduced in Donegan et al. (2020), which uses Piironen and Vehtari's (2017) regularized horseshoe prior. -ESF decomposes spatial autocorrelation into a linear combination of various patterns, typically at different scales (such as local, regional, and global trends). By adding a spatial filter to a regression model, these spatial autocorrelation patterns are shifted from the residuals to the spatial filter. ESF models take the spectral decomposition of a transformed spatial connectivity matrix, \eqn{C}. The resulting eigenvectors, \eqn{E}, are mutually orthogonal and uncorrelated map patterns. The spatial filter equals \eqn{E \beta_{E}} where \eqn{\beta_{E}} is a vector of coefficients. +By adding a spatial filter to a regression model, spatial autocorrelation patterns are shifted from the residuals to the spatial filter. ESF models take the spectral decomposition of a transformed spatial connectivity matrix, \eqn{C}. The resulting eigenvectors, \eqn{E}, are mutually orthogonal and uncorrelated map patterns (at various scales, 'local' to 'regional' to 'global'). The spatial filter equals \eqn{E \beta_{E}} where \eqn{\beta_{E}} is a vector of coefficients. ESF decomposes the data into a global mean, \eqn{\alpha}, global patterns contributed by covariates \eqn{X \beta}, spatial trends \eqn{E \beta_{E}}, and residual variation. Thus, for \code{family=gaussian()}, @@ -155,15 +155,7 @@ ESF decomposes the data into a global mean, \eqn{\alpha}, global patterns contri y \sim Gauss(\alpha + X * \beta + E \beta_{E}, \sigma). } -An ESF component can be incorporated into the linear predictor of any generalized linear model. For example, a spatial Poisson model for rare disease incidence may be specified as follows: - -\deqn{y \sim Poisson(e^{O + \mu})} -\deqn{\mu = \alpha + E \beta_{E} + A } -\deqn{ A \sim Guass(0, \tau) } -\deqn{ \tau \sim student(20, 0, 2) } -\deqn{ \beta_{E} \sim horseshoe(.) } - -The form of this model is similar to the BYM model (see \link[geostan]{stan_icar}), in the sense that it contains a spatially structured trend term (\eqn{E \beta_{E}}) and an unstructured ('random effects') term (\eqn{A}). +An ESF component can be incorporated into the linear predictor of any generalized linear model. For example, using \code{stan_esf} with \code{family = poisson()} and adding a 'random effects' term for each spatial unit (via the \code{re} argument) will produce a model that resembles the BYM model (combining spatially structured and spatially-unstructured components). The \link[geostan]{spatial.geostan_fit} method will return \eqn{E \beta_{E}}. @@ -172,7 +164,7 @@ The model can also be extended to the space-time domain; see \link[geostan]{shap The coefficients \eqn{\beta_{E}} are assigned the regularized horseshoe prior (Piironen and Vehtari, 2017), resulting in a relatively sparse model specification. In addition, numerous eigenvectors are automatically dropped because they represent trace amounts of spatial autocorrelation (this is controlled by the \code{threshold} argument). By default, \code{stan_esf} will drop all eigenvectors representing negative spatial autocorrelation patterns. You can change this behavior using the \code{nsa} argument. \subsection{Additional functionality}{ -The CAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for \link[geostan]{stan_glm}. +The ESF models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for \link[geostan]{stan_glm}. } } \examples{ diff --git a/man/stan_icar.Rd b/man/stan_icar.Rd index 90d0d7b4..4e9db302 100644 --- a/man/stan_icar.Rd +++ b/man/stan_icar.Rd @@ -242,7 +242,7 @@ scale_c <- function(C) \{ \subsection{Additional functionality}{ -The CAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for \link[geostan]{stan_glm}. +The ICAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), missing outcome data, and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for \link[geostan]{stan_glm}. } } \examples{ @@ -264,30 +264,6 @@ sp_diag(fit.bym, sentencing) library(rstan) rstan::stan_ess(fit.bym$stanfit) rstan::stan_rhat(fit.bym$stanfit) - -# calculate log-standardized incidence ratios -# (observed/exected case counts) -library(ggplot2) -library(sf) - -f <- fitted(fit.bym, rates = FALSE)$mean -SSR <- f / sentencing$expected_sents -log.SSR <- log( SSR, base = 2) - -ggplot(sentencing) + - geom_sf(aes(fill = log.SSR)) + - scale_fill_gradient2( - low = "navy", - high = "darkred" - ) + - labs(title = "Log-standardized incidence ratios", - subtitle = "log( Fitted/Expected), base 2") + - theme_void() + - theme( - legend.position = "bottom", - legend.key.height = unit(0.35, "cm"), - legend.key.width = unit(1.5, "cm") - ) } } \seealso{ diff --git a/tests/check-ME-estimates-monte-carlo.R b/tests/check-ME-estimates-monte-carlo.R index d572dbe0..9b1634a4 100644 --- a/tests/check-ME-estimates-monte-carlo.R +++ b/tests/check-ME-estimates-monte-carlo.R @@ -5,8 +5,11 @@ library(geostan) +# flag for running Monte Carlo test or just code once +full_test <- FALSE + ## no. iterations -M = 15 +M <- ifelse(full_test, 20, 1) ## a regular grid ncol = 20 diff --git a/tests/check-SEM-CAR-GLM-estimates-monte-carlo.R b/tests/check-SEM-CAR-GLM-estimates-monte-carlo.R index c914b23b..69f70941 100644 --- a/tests/check-SEM-CAR-GLM-estimates-monte-carlo.R +++ b/tests/check-SEM-CAR-GLM-estimates-monte-carlo.R @@ -11,8 +11,11 @@ #devtools::load_all("~/dev/geostan") library(geostan) +# flag for running Monte Carlo test or just code once +full_test <- FALSE + ## no. iterations -M = 20 +M <- ifelse(full_test, 20, 1) # generate y using SLM or SEM #sar_type <- commandArgs(trailingOnly = TRUE) diff --git a/tests/check-SLM-estimates-monte-carlo-output.csv b/tests/check-SLM-estimates-monte-carlo-output.csv index 78c7d124..bc01387b 100644 --- a/tests/check-SLM-estimates-monte-carlo-output.csv +++ b/tests/check-SLM-estimates-monte-carlo-output.csv @@ -3,7 +3,7 @@ "beta",-0.5,-0.5,0 "gamma",0.1,0.1,0 "rho",0.7,0.7,0 -"sigma",0.1,0.1,0 +"sigma",0.25,0.23,-0.02 "direct",-0.57,-0.57,0 -"indirect",-0.76,-0.76,0 -"total",-1.33,-1.33,0 +"indirect",-0.76,-0.78,-0.01 +"total",-1.33,-1.34,-0.01 diff --git a/tests/check-SLM-estimates-monte-carlo.R b/tests/check-SLM-estimates-monte-carlo.R index 3d0f0dc4..32f1f679 100644 --- a/tests/check-SLM-estimates-monte-carlo.R +++ b/tests/check-SLM-estimates-monte-carlo.R @@ -8,12 +8,15 @@ #library(spatialreg) library(geostan) -# no. iterations [s > 90 required for impact est. to converge - they are high variance; others converge by 30] -S <- 30 +# flag for running Monte Carlo test or just code once +full_test <- FALSE # use measurement error in x has_me <- FALSE +# no. iterations +M <- ifelse(full_test, 20, 1) + # regular lattice parts <- prep_sar_data2(row = 10, col = 30, quiet = TRUE) cp = prep_car_data2(10, 30, quiet=T) @@ -28,16 +31,16 @@ N <- nrow(W) B <- -0.5 G <- 0.1 R <- 0.7 -Sig = 0.1 -sigma_me <- 0.1 +Sig = 0.25 +sigma_me <- 0.2 pars <- c(const = 0, beta = B, gamma = G, rho = R, sigma = Sig) pars <- c(pars, geostan::spill(beta = B, gamma = G, rho = R, W = W, approx=FALSE)) # results for S iterations ## res <- res2 <- matrix(NA, nrow = S, ncol = 8) -res <- matrix(NA, nrow = S, ncol = 8) +res <- matrix(NA, nrow = M, ncol = 8) -for (s in 1:S) { +for (m in 1:M) { x <- sim_sar(w=W, rho=R) if (has_me) x = x + rnorm(N, sd = sigma_me) Wx <- (W %*% x)[,1] @@ -70,7 +73,7 @@ for (s in 1:S) { suppressWarnings() } - res[s, ] <- c( + res[m, ] <- c( fit$summary[c('intercept', 'x', 'w.x', 'sar_rho', 'sar_scale'), 'mean'], geostan::impacts(fit, approx = FALSE)$summary$x[,'mean'] ) @@ -89,7 +92,7 @@ for (s in 1:S) { # spatialreg reports variance parameter, geostan reports scale param. ## res2[, 5] <- sqrt(res2[,5]) -cat("\n**\nSDLM Monte Carlo results\n", S, "iterations\nAverage values\n**\n") +cat("\n**\nSDLM Monte Carlo results\n", M, "iterations\nAverage values\n**\n") est <- apply(res, 2, mean)