# Case study on spatial models for areal data - Poisson CAR/IAR

**URL:** <https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481>\
**Category:** Modeling\
**Created:** [May 9, 2017, 8:45pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481 "2017-05-09T20:45:03Z")\
**Posts on this page:** 17\
**Page:** 6

<div class="post-metadata">

**Author:** ![anon75146577](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@anon75146577](https://discourse.mc-stan.org/u/anon75146577)\
**Post date:** [September 18, 2017, 1:34pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/101 "2017-09-18T13:34:26Z")

</div>

There’s nothing stopping us from coding that. It’s just that I don’t know how to automatically detect the number of connected components / what’s in each connected component. But I’m very bad with graphs!

---

<div class="post-metadata">

**Author:** ![fabio](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fabio/32/935_2.png) [@fabio](https://discourse.mc-stan.org/u/fabio)\
**Post date:** [September 18, 2017, 2:22pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/102 "2017-09-18T14:22:53Z")

</div>

I usually use the _spdep_ package. Or is it something more complicated?

---

<div class="post-metadata">

**Author:** ![anon75146577](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@anon75146577](https://discourse.mc-stan.org/u/anon75146577)\
**Post date:** [September 18, 2017, 3:38pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/103 "2017-09-18T15:38:38Z")

</div>

Which function do you use? (Sorry - I don’t usually do these things manually)

---

<div class="post-metadata">

**Author:** ![fabio](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fabio/32/935_2.png) [@fabio](https://discourse.mc-stan.org/u/fabio)\
**Post date:** [September 19, 2017, 11:45am UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/104 "2017-09-19T11:45:27Z")

</div>

@anon75146577, Here is a little example.  
I use the package _maptools_ to have the ‘Sudden Infant deaths’ shapefile available, then use the package _rgdal_ to read the shapefiles from my PC by using _readOGR_

```
library(maptools)
library(rgdal)
sids <- readOGR(system.file("shapes/", package="maptools"), "sids")

```

Then I use the package _spdep_ to build the adjacency

```
library(spdep)
temp <- poly2nb(sids)

```

if you want an higher order lag

```
temp2= nblag(temp,2)

```

To build the adjacency matrix just use the _nb2mat_ function

```
matrix <- nb2mat(temp, style="B")

```

_style=“B”_ is the basic binary coding, W is row standardised, C is globally standardised, while S is the variance-stabilizing coding scheme proposed by Tiefelsdorf et al. 1999. W£ith this matrix you can check who is connected with whom, and who is disconnected.

Anyway, for a more extensive description you can check [https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf](https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf)  
or the manual of the _spdep_ package (it’s huge)  
I hope this will help!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 19, 2017, 8:50pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/105 "2017-09-19T20:50:48Z")

</div>

Easy peasy. I missed the matrix classes, but took _a lot_ of algorithms classes.

> **[Tarjan's strongly connected components algorithm](https://en.wikipedia.org/wiki/Tarjan&#39;s_strongly_connected_components_algorithm)**
>
> Tarjan's algorithm is an algorithm in graph theory for finding the strongly connected components of a graph. It runs in linear time, matching the time bound for alternative methods including Kosaraju's algorithm and the path-based strong component algorithm. Tarjan's Algorithm is named for its inventor, Robert Tarjan. The algorithm takes a directed graph as input, and produces a partition of the graph's vertices into the graph's strongly connected components. Each vertex of the graph appears in...

I could code that as a Stan function, but it may not be worthwhile to do it in Stan directly unless you somehow need the components themselves.

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [September 27, 2017, 8:15pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/106 "2017-09-27T20:15:06Z")

</div>

> [@anon75146577](#):
>
> There’s nothing stopping us from coding that. It’s just that I don’t know how to automatically detect the number of connected components / what’s in each connected component. But I’m very bad with graphs!

R’s spdep package does this for you, function `n.comp.nb`

> **[spdep source: R/components.R](https://rdrr.io/rforge/spdep/src/R/components.R)**
>
> R/components.R defines the following functions: n.comp.nb

given an `nb` object `nb_nyc_mm`, here a map of New York City census tracts:

```no-highlight
> components <- n.comp.nb(nb_nyc_mm)
> table(components$comp.id)

   1 2 3 4 5 6 7 
1983 2 1 1 22 1 108 

```

there are 7 components, 3 of which are islands.

---

<div class="post-metadata">

**Author:** ![anon75146577](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@anon75146577](https://discourse.mc-stan.org/u/anon75146577)\
**Post date:** [September 27, 2017, 8:17pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/107 "2017-09-27T20:17:43Z")

</div>

Cool! In that case there’s no reason not to code up Anna and Håvard’s stuff. @imadmali: you should consider this for rstanarm.

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [September 27, 2017, 8:19pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/108 "2017-09-27T20:19:56Z")

</div>

are islands a problem here?

---

<div class="post-metadata">

**Author:** ![anon75146577](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@anon75146577](https://discourse.mc-stan.org/u/anon75146577)\
**Post date:** [September 27, 2017, 8:23pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/109 "2017-09-27T20:23:13Z")

</div>

Kind of. The nullspace of the model with one connected component is just a constant.

When there are multiple connected components, you can add a constant to any of them independently. So there should be 7 constraints on the ICAR over the nb\_nyc\_mm graph.

The paper [linked above](https://arxiv.org/abs/1705.04854) gives good advice about what to do in that situation. (In particular, you should do something specific for the islands: basically take them out of the ICAR and give them a random effect)

---

<div class="post-metadata">

**Author:** ![mattwilliamson13](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattwilliamson13/32/773_2.png) [@mattwilliamson13](https://discourse.mc-stan.org/u/mattwilliamson13)\
**Post date:** [July 3, 2018, 10:04pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/110 "2018-07-03T22:04:31Z")

</div>

This case study is extremely helpful. I have one question/clarification (sorry to resurrect such an old thread). @mitzimorris notes in her description of the bym2 model results that:

> n\_eff numbers, while low for rho and sigma, are sufficient

I’m wondering how to gauge “sufficient” in this context. I ask because I’ve been running a similar model (albeit binomial rather than Poisson) and am getting pretty low (\<100) `n_eff` values for `rho` and a number of `phi`’s (despite running 8000 warmup iterations and 10000 iterations total). The `Rhat` values are all ok (~1), is that an indication that the n\_eff are sufficient?

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [July 3, 2018, 11:38pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/111 "2018-07-03T23:38:59Z")

</div>

are you running a single chain with 8000 warmup iterations or are you running multiple chains? how many warmup iterations per chain?

if you’re using R, the [Bayesplot package](http://mc-stan.org/users/interfaces/bayesplot) might be useful. Also see Michael Betancourt’s [RStan workflow case study](http://mc-stan.org/users/documentation/case-studies/rstan_workflow.html) - there’s a [utility script](https://github.com/betanalpha/knitr_case_studies/blob/master/rstan_workflow/stan_utility.R) which has a function to check n\_eff.

---

<div class="post-metadata">

**Author:** ![mattwilliamson13](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattwilliamson13/32/773_2.png) [@mattwilliamson13](https://discourse.mc-stan.org/u/mattwilliamson13)\
**Post date:** [July 3, 2018, 11:46pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/112 "2018-07-03T23:46:28Z")

</div>

I’m running 4 chains with `warmup = 8000` and `iter=10000` which I presume means each chain is running that number of warmup and iterations. I’ll give that utility script a look. Thanks for any additional info.

My model has a varying intercept in addition to the CAR term. I’ve tried running it with the varying intercept included in the `convolved_re` term and with it maintained in the actual model with an additional random effect (`theta`) included in the `convolved_re` term. Perhaps that is causing my problem?

Thanks for your help!

---

<div class="post-metadata">

**Author:** ![mattwilliamson13](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattwilliamson13/32/773_2.png) [@mattwilliamson13](https://discourse.mc-stan.org/u/mattwilliamson13)\
**Post date:** [July 4, 2018, 1:11am UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/113 "2018-07-04T01:11:45Z")

</div>

Just ran the utility script for `check_n_eff` and returned:

> [1] “n\_eff / iter looks reasonable for all parameters”

so I guess I’m good. I’ll need to read a bit more to understand why I’m good, but at least it seems safe to proceed. Thanks for your help!

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [July 4, 2018, 7:34am UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/114 "2018-07-04T07:34:01Z")

</div>

The `check_n_eff` function simply checks to see if you can trust the effective sample size _estimator_, not whether the effective sample size is sufficient in your particular application. In practice we have to estimate effective sample size from the empirical autocorrelation in the samples, and if the autocorrelations are too large the estimator can be biased high and strongly overestimate the true effective sample size.

Whether the effective sample size in your chains is sufficient depends on the details of your application. In ideal circumstances (where you don’t see any divergences or E-FMI warnings, see the other functions in `stan_utility.R`) the effective sample size controls how accurate your Markov chain Monte Carlo estimators are – more effective samples, more precise estimators. How much precision you need depends on the questions you want to answer.

The case study at [https://betanalpha.github.io/assets/case\_studies/rstan\_workflow.html](https://betanalpha.github.io/assets/case_studies/rstan_workflow.html) more carefully goes through these definitions. See also [https://arxiv.org/abs/1701.02434](https://arxiv.org/abs/1701.02434), Section 2 for a more elaborate discussion.

---

<div class="post-metadata">

**Author:** ![mattwilliamson13](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattwilliamson13/32/773_2.png) [@mattwilliamson13](https://discourse.mc-stan.org/u/mattwilliamson13)\
**Post date:** [July 4, 2018, 6:13pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/115 "2018-07-04T18:13:53Z")

</div>

Thanks @betanalpha. The other utility checks don’t report any additional pathologies. For my application `rho` and `phi` are not actually the parameters of interest (I’m interested in the `beta` values after conditioning on the potential for spatial dependence in the data). Which leads me back to trying to understand what “sufficient” means in this context:

> [@mattwilliamson13](#):
>
> n\_eff numbers, while low for rho and sigma, are sufficient

Am I to assume that estimates for `rho` and `phi` are sufficient if the model appears to be mixing and free of other pathologies given that my interest is in the `beta` estimates conditional on the existence of spatial structure? The mean’s and sd’s for the `beta` parameters seem reasonable given previous iterations of the model without spatial structure. Or am I missing some key element that determines “sufficiency”?

Thanks to both you and @mitzimorris for your help on this.

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [July 4, 2018, 8:23pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/116 "2018-07-04T20:23:58Z")

</div>

No other pathologies is great, so we can focus on effective sample size.

The worry in a model like this is that the latent parameters \rho and \sigma freeze and you end up sampling from the conditional distribution \pi(\theta \mid y, \rho, \sigma) instead of the joint distribution \pi(\theta, \rho, \sigma | y). But given all of the diagnostic information that doesn’t seem to be the case.

If you case about the \beta then what you have to figure out is how accurately do you need to know the mean of each \beta? Remember that MCMC only _estimates_ that mean, and the error in the estimate depends on the effective sample size. In `RStan` you see this as the `mcse`, the second column when you call `print(fit)`. If the `mcse` is small enough for your inferences then great, you’re done. If it’s not precise enough you’ll have to run longer or try to implement your model in a way that samples the latent parameters more efficiently (not that a superior implementation necessarily exists).

In particular, the threshold for “precise enough” depends on _your application_. We can’t tell you what is sufficient, you have to figure that out within the context of your domain expertise. That said, as a rule of thumb O(100) effective samples is typically sufficient for most applications.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 12, 2018, 11:08pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/117 "2018-07-12T23:08:27Z")

</div>

> [@mattwilliamson13](#):
>
> I’m running 4 chains with `warmup = 8000` and `iter=10000` which I presume means each chain is running that number of warmup and iterations.

That should run 8000 warmup iterations and 2000 sampling iterations.

[Previous page](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481.md?page=5)