# 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:** 20\
**Page:** 5

<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:** [August 22, 2017, 5:34pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/81 "2017-08-22T17:34:16Z")

</div>

Hi Katherine,

Mitzi’s model assumes that the graph has only one connected component. So it doesn’t work with “islands”. In the event that you do have islands, the scaling will be wrong, but more critically, the constraints are wrong. You need to make it sum to zero on _each_ connected component. INLA does all of this automatically, but for the moment you’ve got to do it manually in Stan.

Those alternate adjacency methods seem like a good idea.

---

<div class="post-metadata">

**Author:** ![MigueBeneito](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/miguebeneito/32/18418_2.png) [@MigueBeneito](https://discourse.mc-stan.org/u/MigueBeneito)\
**Post date:** [August 23, 2017, 11:37am UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/82 "2017-08-23T11:37:46Z")

</div>

Hi Bob

Thanks for your answer

> [@Bob\_Carpenter](#):
>
> You probably want to do much more warmup than 400 iterations. If you start with a small step size, set a high target acceptance rate, and warmup longer, you’ll get much more stable sampling times and probalby higher n\_eff overall per unit time.

This is interesting. I usually take a burn-in period of about 10% of the whole simulation. In my opinion this is a safer implementation of Geyer’s advice (Geyer (1992). Practical Markov Chain Monte-Carlo. Statistical Science, 7(4):473-511) who propose to set a burn-in periodo of 1-2% of the total simulation. In my experience with WinBUGS this is a reasonable advice, surely more efficient than taking 50% of the simulations as burn-in. But, is a 10% burn-in period an adequate choice for Stan? Does HMCMC requires/advices for longer burn-in periods? Does Stan tune its HMCMC algorithms during the burn-in period and setting a short warm-up period could harm Stan’s performance?

Apologizes if this is a bit off-topic but I have prefered to ask it here in response to Bob’s response :).

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 23, 2017, 12:08pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/83 "2017-08-23T12:08:57Z")

</div>

> [@MigueBeneito](#):
>
> Does Stan tune its HMCMC algorithms during the burn-in period and setting a short warm-up period could harm Stan’s performance?

Yes.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [August 23, 2017, 12:31pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/84 "2017-08-23T12:31:15Z")

</div>

With HMC/NUTS the warm-up period includes tuning the mass matrix and stepsize in addition to discovering the main mass of the distribution so yeah, for most interesting models the full warm-up is a good idea. The up side is that you only need a small number of iterations to get the required effective sample size. Even in complex models I often need only 1-2 iterations per effective sample (yeah, sometimes it’s not _that_ good, but those are cases where M-H would fail completely). so if you’re getting a preliminary estimate of a mean you only need 30-100 iterations post-warm-up, and maybe 300 for credible intervals. For final estimates you often need more to get the MCMC error down.

---

<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:** [August 23, 2017, 11:13pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/85 "2017-08-23T23:13:21Z")

</div>

Don’t follow Geyer on this, as he’s making a lot of assumptions that make sense theoretically, but not practically. If you have a chain that you know is going to converge (e.g., it satisfies geometric ergodicity) and you are going to run it for a bajillion iterations, then his advice is OK, but still dangerous. For real world problems, you should be running multiple chains.

Our reasoning was that you’ll lose at most a factor of two in speed by having warmup be half the iterations.

As @sakrejda hints, generally, people run for way more effective samples than they need, too.

---

<div class="post-metadata">

**Author:** ![kwheeler-martin](https://avatars.discourse-cdn.com/v4/letter/k/f4b2a3/32.png) [@kwheeler-martin](https://discourse.mc-stan.org/u/kwheeler-martin)\
**Post date:** [August 24, 2017, 8:55pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/86 "2017-08-24T20:55:40Z")

</div>

Hi Daniel and Mitzi,

Thank you for the clarification. I’m currently looking at the “5 nearest neighbors” approach with NYC Planning’s 2010 Census Tract (clipped to shoreline). I made one manual connection between Brooklyn and Staten Island by the Verrazano Bridge, yielding a fully connected map.

However, in the final step prior to running the model, where the adjacency matrix and scaling are set up, I’m now getting an error message as follows. I’m stumped as to what it means:

> This is INLA\_17.06.20 built 2017-06-20 03:42:30 UTC.  
> See [R-INLA Project - Contact us](http://www.r-inla.org/contact-us) for how to get help.  
> Error in validObject(r) :  
> invalid class “dsTMatrix” object: all row indices (slot ‘i’) must be between 0 and nrow-1 in a TsparseMatrix  
> In addition: Warning messages:  
> 1: package ‘INLA’ was built under R version 3.4.0  
> 2: In fun(libname, pkgname) : bytecode version mismatch; using eval

Here is the code I used:

> #Run ICAR model using STAN - developed by Mitzi Morris et al  
> library(rstan);  
> options(mc.cores = parallel::detectCores());  
> library(maptools);  
> library(spdep);  
> library(tripack);  
> library(rgdal)
> 
> nyckidped\_tract\<-readRDS(file=“D:/NYC/nyckidped\_tract.rds”)  
> #source(“nyckidped\_tract.rds”);  
> y = nyckidped\_tract$count;  
> x = nyckidped\_tract$expected;
> 
> nyc\_poptractIDs\<-nyckidped\_tract$ct2010full  
> head(nyc\_poptractIDs)
> 
> nyc\_all\_tracts.shp\<-readOGR(“D:/Census/nycTracts10/nyct2010.shp”, layer=“nyct2010”);  
> print(nyc\_all\_tracts.shp)  
> table(nyc\_all\_tracts.shp$CT2010)  
> table(nyc\_all\_tracts.shp$BoroCode, nyc\_all\_tracts.shp$BoroName)
> 
> nyc\_all\_tracts.shp$BoroFips\<-ifelse(nyc\_all\_tracts.shp$BoroCode==1, “061”, ifelse(nyc\_all\_tracts.shp$BoroCode==2, “005”, ifelse(nyc\_all\_tracts.shp$BoroCode==3, “047”, ifelse(nyc\_all\_tracts.shp$BoroCode==4, “081”, ifelse(nyc\_all\_tracts.shp$BoroCode==5, “085”, “”)))))  
> table(nyc\_all\_tracts.shp$BoroName, nyc\_all\_tracts.shp$BoroFips)
> 
> nyc\_all\_tracts.shp$CT2010full\<-paste0(“36”,nyc\_all\_tracts.shp$BoroFips,nyc\_all\_tracts.shp$CT2010)  
> table(nyc\_all\_tracts.shp$CT2010full)
> 
> pop\_tracts ← nyc\_all\_tracts.shp$CT2010full %in% nyc\_poptractIDs;  
> nyc\_poptracts.shp ← nyc\_all\_tracts.shp[pop\_tracts,];  
> nyc\_poptracts.shp ← nyc\_poptracts.shp[order(nyc\_poptracts.shp$BoroCT2010),];
> 
> #nearest neighbors  
> nb\_nyc\_5nn ← knn2nb(knearneigh(coords,k=5), row.names=IDs);  
> plot(nb\_nyc\_5nn, coords, pch = “.”, col = “red”);

> #Manually add connection between Brooklyn and Staten Island  
> nb\_nyc\_5nn1n \<-edit.nb(nb\_nyc\_5nn, coords)  
> #Identifying contiguity for deletion …  
> #No contiguity between chosen points  
> #Add contiguity? (y/n) y  
> #added contiguity between point 733 and 2018  
> #Options: quit[q] refresh[r] continue[c] r  
> #Options: quit[q] continue[c]q
> 
> plot(nb\_nyc\_5nn1n, coords, pch = “.”, col = “red”);
> 
> #source(“nbdata4stan.R”);  
> nbdata4stan = function(x) {  
> N = length(x);  
> n\_links = 0;  
> for (i in 1:N) {  
> if (x[[i]][1] != 0) {  
> n\_links = n\_links + length(x[[i]]);  
> }  
> }  
> N\_edges = n\_links / 2;  
> node1 = vector(mode=“numeric”, length=N\_edges);  
> node2 = vector(mode=“numeric”, length=N\_edges);  
> idx = 0;  
> for (i in 1:N) {  
> if (x[[i]][1] != 0) {  
> for (j in 1:length(x[[i]])) {  
> n2 = unlist(x[[i]][j]);  
> if (i \< n2) {  
> idx = idx + 1;  
> node1[idx] = i;  
> node2[idx] = n2;  
> }  
> }  
> }  
> }  
> return (list(“N”=N,“N\_edges”=N\_edges,“node1”=node1,“node2”=node2));  
> }

> nbs=nbdata4stan(nb\_nyc\_5nn1n);  
> N = nbs$N;  
> node1 = nbs$node1;  
> node2 = nbs$node2;  
> N\_edges = nbs$N\_edges;
> 
> DO\_CALC = TRUE;  
> #Calculate the scaling of the model. Requires the INLA package.
> 
> if(DO\_CALC) {  
> library(INLA)  
> #Build the adjacency matrix  
> adj.matrix = sparseMatrix(i=nbs$node1,j=nbs$node2,x=1,symmetric=TRUE)  
> #The ICAR precision matrix (note! This is singular)  
> Q= Diagonal(nbs$N, rowSums(adj.matrix)) - adj.matrix
> 
> #Add a small jitter to the diagonal for numerical stability (optional but recommended)  
> Q\_pert = Q + Diagonal(nbs$N) \* max(diag(Q)) \* sqrt(.Machine$double.eps)
> 
> # Compute the diagonal elements of the covariance matrix subject to the
> 
> # constraint that the entries of the ICAR sum to zero.
> 
> #See the function help for further details.  
> Q\_inv = inla.qinv(Q\_pert, constr=list(A = matrix(1,1,nbs$N),e=0))
> 
> #Compute the geometric mean of the variances, which are on the diagonal of Q.inv  
> scaling\_factor = exp(mean(log(diag(Q\_inv))))  
> } else{  
> scaling\_factor= 0.5;  
> }
> 
> bym2\_nyc\_stanfit = stan(“D:/NYC/bym2\_predictor\_only.stan”, data=list(N,N\_edges,node1,node2,y,x,scaling\_factor), iter=6000);  
> print(bym2\_nyc\_stanfit, digits=3, pars=c(“lp\_\_”, “beta0”, “beta1”, “rho”, “sigma”, “mu[1]”, “mu[2]”, “mu[3]”, “mu[2116]”, “mu[2117]”, “mu[2118]”), probs=c(0.025, 0.5, 0.975));

Katie

---

<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:** [August 24, 2017, 8:58pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/87 "2017-08-24T20:58:23Z")

</div>

Can you type `traceback()` and send the output ?

Edit: also `str(Q)`

---

<div class="post-metadata">

**Author:** ![kwheeler-martin](https://avatars.discourse-cdn.com/v4/letter/k/f4b2a3/32.png) [@kwheeler-martin](https://discourse.mc-stan.org/u/kwheeler-martin)\
**Post date:** [August 25, 2017, 1:07pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/88 "2017-08-25T13:07:11Z")

</div>

> traceback()
> 
> 3: stop(msg, ": ", errors, domain = NA)  
> 2: validObject(r)  
> 1: sparseMatrix(i = nbs$node1, j = nbs$node2, x = 1, symmetric = TRUE) at #4
> 
> str(Q)  
> Error in str(Q) : object ‘Q’ not found

---

<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:** [August 25, 2017, 1:17pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/89 "2017-08-25T13:17:22Z")

</div>

> [@kwheeler-martin](#):
>
> library(INLA)
> 
> #Build the adjacency matrix
> 
> adj.matrix = sparseMatrix(i=nbs$node1,j=nbs$node2,x=1,symmetric=TRUE)
> 
> #The ICAR precision matrix (note! This is singular)
> 
> Q= Diagonal(nbs$N, rowSums(adj.matrix)) - adj.matrix
> 
> #Add a small jitter to the diagonal for numerical stability (optional but recommended)
> 
> Q\_pert = Q + Diagonal(nbs$N) \* max(diag(Q)) \* sqrt(.Machine$double.eps)

Ok the Q wasn’t made which means that this bit of the code wasn’t run

---

<div class="post-metadata">

**Author:** ![kwheeler-martin](https://avatars.discourse-cdn.com/v4/letter/k/f4b2a3/32.png) [@kwheeler-martin](https://discourse.mc-stan.org/u/kwheeler-martin)\
**Post date:** [August 25, 2017, 3:23pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/90 "2017-08-25T15:23:37Z")

</div>

> [@anon75146577](#):
>
> library(INLA)
> 
> #Build the adjacency matrix
> 
> adj.matrix = sparseMatrix(i=nbs$node1,j=nbs$node2,x=1,symmetric=TRUE)
> 
> #The ICAR precision matrix (note! This is singular)
> 
> Q= Diagonal(nbs$N, rowSums(adj.matrix)) - adj.matrix
> 
> #Add a small jitter to the diagonal for numerical stability (optional but recommended)
> 
> Q\_pert = Q + Diagonal(nbs$N) \* max(diag(Q)) \* sqrt(.Machine$double.eps)

It gives the following error, again. Perhaps the output from the nearest neighbor function differs from poly2nb? I’ll have to check more closely at a later time.

> Error in validObject(r) :  
> invalid class “dsTMatrix” object: all row indices (slot ‘i’) must be between 0 and nrow-1 in a TsparseMatrix

---

<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:** [August 25, 2017, 5:08pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/91 "2017-08-25T17:08:27Z")

</div>

If you email me enough to run the code (dp(dot)simpson(at)gmail(dot)com) I’ll look at it for you.

---

<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:** [August 28, 2017, 6:46pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/92 "2017-08-28T18:46:43Z")

</div>

> [@anon75146577](#):
>
> dp(dot)simpson(at)gmail(dot)com

You really need to do this kind of thing with gmail? Is their spam filter that bad?

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 28, 2017, 7:12pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/93 "2017-08-28T19:12:49Z")

</div>

I’ve been told more than ten years ago that email harvesters can decode (dot) and (at)

---

<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:** [August 28, 2017, 10:32pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/94 "2017-08-28T22:32:37Z")

</div>

It wouldn’t be exactly rocket science to do, would it

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [August 29, 2017, 12:06pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/95 "2017-08-29T12:06:54Z")

</div>

This might be getting OT and maybe needs its own thread… ¯\_༼ ି ~ ି ༽_/¯

---

<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:** [August 29, 2017, 6:54pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/96 "2017-08-29T18:54:37Z")

</div>

> [@kwheeler-martin](#):
>
> It gives the following error, again. Perhaps the output from the nearest neighbor function differs from poly2nb? I’ll have to check more closely at a later time.

Ok. I know where the problem is. And I _think_ I know how to fix it.

The problem is that `nbdata4stan` assumes it gets a list of all the neighbours for each node (excluding the node number itself). It then makes an assumption about symmetry - namely that if node i has neighbour j, then node j has neighbour i and codes everything accordingly. _This is not true for output for knn2nb!_ You can see this by running

> is.symmetric.nb( knn2nb(knearneigh(coords,k=5), row.names=IDs)) #Output: FALSE

The fix is to replace

> nb\_nyc\_5nn ← knn2nb(knearneigh(coords,k=5), row.names=IDs); #bad

with

> nb\_nyc\_5nn ← knn2nb(knearneigh(coords,k=5), row.names=IDs, sym=TRUE); #good

With this change, the rest of the function works.

**NB: This was tested on `nb_nyc_5nn` rather than `nb_nyc_5nn1` because I couldn’t make the latter work.**

You can always use the functions `is.symmetric.nb` (desired output: `TRUE`) and `make.sym.nb` (in the event that the nb is not symmetric) to make sure everythign is still ok with `nb_nyc_5nn1`.)

EDIT: Here is the “safe” version of `nbdata4stan`

> nbdata4stan = function(x) {
> 
> if (! is.symmetric.nb(x)){  
> stop(“The neihbourhood you provided is not symmetric. Use make.sym.nb to fix this!”)  
> }
> 
> N = length(x);  
> n\_links = 0;  
> for (i in 1:N) {  
> if (x[[i]][1] != 0) {  
> n\_links = n\_links + length(x[[i]]);  
> }  
> }  
> N\_edges = n\_links / 2;  
> node1 = vector(mode=“numeric”, length=N\_edges);  
> node2 = vector(mode=“numeric”, length=N\_edges);  
> idx = 0;  
> for (i in 1:N) {  
> if (x[[i]][1] != 0) {  
> idx\_tmp=idx  
> for (j in 1:length(x[[i]])) {  
> n2 = unlist(x[[i]][j]);
> 
> ```
> if (i < n2) { 
> idx = idx + 1; 
> node1[idx] = i; 
> node2[idx] = n2; 
> } 
> } 
> if ( idx_tmp==idx) { print(i)}
> } 
> 
> ```
> 
> }  
> return (list(“N”=N,“N\_edges”=N\_edges,“node1”=node1,“node2”=node2));  
> }

---

<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:** [August 29, 2017, 8:31pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/97 "2017-08-29T20:31:03Z")

</div>

> Ok. I know where the problem is. And I think I know how to fix it.

many thanks for identifying and fixing the problem!  
I’ll update the code for this function in the case study accordingly!

---

<div class="post-metadata">

**Author:** ![kwheeler-martin](https://avatars.discourse-cdn.com/v4/letter/k/f4b2a3/32.png) [@kwheeler-martin](https://discourse.mc-stan.org/u/kwheeler-martin)\
**Post date:** [September 7, 2017, 4:16pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/98 "2017-09-07T16:16:51Z")

</div>

It indeed ran successfully, and the speed improved as well (~20 min). The results look much better:

> print(bym2\_nyc\_stanfit, digits=3, pars=c(“lp\_\_”, “beta0”, “beta1”, “rho”, “sigma”, “mu[1]”, “mu[2]”, “mu[3]”, “mu[4]”, “mu[5]”, “mu[1001]”, “mu[1002]”,“mu[1003]”,“mu[1004]”,“mu[1005]”, “mu[2114]”, “mu[2115]”, “mu[2116]”, “mu[2117]”, “mu[2118]”), probs=c(0.025, 0.5, 0.975));  
> Inference for Stan model: bym2\_predictor\_only.  
> 4 chains, each with iter=6000; warmup=3000; thin=1;  
> post-warmup draws per chain=3000, total post-warmup draws=12000.

```
> mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
> lp__ 15069.600 0.985 57.223 14960.410 15069.149 15182.101 3372 1.000
> beta0 1.073 0.001 0.036 1.002 1.074 1.145 3066 1.000
> beta1 0.078 0.000 0.004 0.070 0.078 0.087 2716 1.000
> rho 0.578 0.002 0.055 0.468 0.579 0.682 537 1.005
> sigma 0.877 0.002 0.038 0.807 0.875 0.956 605 1.004
> mu[1] 7.641 0.026 2.807 3.382 7.239 14.229 12000 1.000
> mu[2] 12.552 0.035 3.840 6.442 12.050 21.250 12000 1.000
> mu[3] 15.450 0.040 4.383 8.206 14.980 25.361 12000 1.000
> mu[4] 16.739 0.045 4.975 8.895 16.157 28.099 12000 1.000
> mu[5] 17.434 0.042 4.547 9.881 16.965 27.448 12000 1.000
> mu[1001] 5.174 0.020 2.230 1.971 4.815 10.450 12000 1.000
> mu[1002] 11.701 0.035 3.884 5.523 11.236 20.666 12000 1.000
> mu[1003] 12.255 0.035 3.857 6.135 11.833 21.126 12000 1.000
> mu[1004] 3.047 0.014 1.519 0.940 2.750 6.756 12000 1.000
> mu[1005] 5.419 0.020 2.243 2.092 5.064 10.702 12000 1.000
> mu[2114] 6.920 0.024 2.673 2.877 6.507 13.159 12000 1.000
> mu[2115] 16.810 0.042 4.597 9.260 16.331 27.131 12000 1.000
> mu[2116] 4.596 0.018 1.955 1.707 4.281 9.221 12000 1.000
> mu[2117] 7.212 0.024 2.677 3.132 6.819 13.523 12000 1.000
> mu[2118] 7.993 0.028 3.016 3.408 7.545 15.094 12000 1.000

```

Correspondingly below, the input count (y) and expected (x) values (i.e., citywide rate applied to census tract population) for the same 15 observations. Data usage note: original crash files were provided by NYC DOT. We assume responsibility for the accuracy of the aggregate data used for analysis, as shown.

```
print(nyckidped_tract[c(1:5,1001:1005,2114:2118),])
      ct2010full geoname years agegroup type count expected
2 36005000200 Census Tract 2, Bronx County, New York 2011to2014 0to19yrs Pedestrian 7 9.712549
3 36005000400 Census Tract 4, Bronx County, New York 2011to2014 0to19yrs Pedestrian 12 12.694855
4 36005001600 Census Tract 16, Bronx County, New York 2011to2014 0to19yrs Pedestrian 15 12.505271
5 36005001900 Census Tract 19, Bronx County, New York 2011to2014 0to19yrs Pedestrian 16 3.930228
6 36005002000 Census Tract 20, Bronx County, New York 2011to2014 0to19yrs Pedestrian 17 21.619900
1017 36047097000 Census Tract 970, Kings County, New York 2011to2014 0to19yrs Pedestrian 5 3.901061
1018 36047097400 Census Tract 974, Kings County, New York 2011to2014 0to19yrs Pedestrian 11 4.783357
1019 36047098200 Census Tract 982, Kings County, New York 2011to2014 0to19yrs Pedestrian 12 11.149014
1020 36047098400 Census Tract 984, Kings County, New York 2011to2014 0to19yrs Pedestrian 2 3.944811
1021 36047098600 Census Tract 986, Kings County, New York 2011to2014 0to19yrs Pedestrian 5 5.534403
2163 36085030301 Census Tract 303.01, Richmond County, New York 2011to2014 0to19yrs Pedestrian 7 9.406297
2164 36085030302 Census Tract 303.02, Richmond County, New York 2011to2014 0to19yrs Pedestrian 17 14.094862
2165 36085031901 Census Tract 319.01, Richmond County, New York 2011to2014 0to19yrs Pedestrian 4 8.866711
2166 36085031902 Census Tract 319.02, Richmond County, New York 2011to2014 0to19yrs Pedestrian 7 13.022982
2167 36085032300 Census Tract 323, Richmond County, New York 2011to2014 0to19yrs Pedestrian 9 2.501054

```

It looks like we are well on our way! As a sensitivity analysis I’ll try running with different neighborhood adjacency methods.

Cheers!

---

<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 10, 2017, 9:29am UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/99 "2017-09-10T09:29:28Z")

</div>

> [@kwheeler-martin](#):
>
> It indeed ran successfully, and the speed improved as well (~20 min). The results look much better:

Great news! Thanks for reporting back with successful run results. It makes us very happy.

---

<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 18, 2017, 1:15pm UTC](https://discourse.mc-stan.org/t/case-study-on-spatial-models-for-areal-data-poisson-car-iar/481/100 "2017-09-18T13:15:17Z")

</div>

> [@anon75146577](#):
>
> Mitzi’s model assumes that the graph has only one connected component. So it doesn’t work with “islands”. In the event that you do have islands, the scaling will be wrong, but more critically, the constraints are wrong. You need to make it sum to zero on each connected component. INLA does all of this automatically, but for the moment you’ve got to do it manually in Stan.

I’ve just found this paper on islands: A note on intrinsic Conditional Autoregressive models for disconnected graphs, Anna Freni-Sterrantino et al, [[1705.04854] A note on intrinsic Conditional Autoregressive models for disconnected graphs](https://arxiv.org/abs/1705.04854)  
can the Stan model be coded-up accordingly?

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

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