updated notebook, improved model code

mitzimorris · mitzimorris · commit 7c7553213f10 · 2021-02-19T20:38:30.000-05:00
diff --git a/knitr/car-iar-poisson/update_2021_02/BYM2 islands.ipynb b/knitr/car-iar-poisson/update_2021_02/BYM2 islands.ipynb
@@ -108,7 +108,22 @@
     "- I: `int<lower = 0> I;  // number of nodes`\n",
     "- J: `int<lower = 0> J;  // number of edges`\n",
     "- edges: `int<lower = 1, upper = I> edges[2, J];  // node[1, j] adjacent to node[2, j]`\n",
-    "- tau: `real tau; // scaling factor`"
+    "- tau: `real tau; // scaling factor`\n",
+    "\n",
+    "The helper function `nb_to_edge_array` takes the `nb` object and returns the 2 $\\times$ J edge array; the helper function `scaling_factor` uses the edge array to compute the geometric mean of the corresponding adjacency matrix.\n",
+    "For the fully connected Scotland graph, the spatial data inputs are:\n",
+    "\n",
+    "```\n",
+    "  \"I\": 56,\n",
+    "  \"J\": 132,\n",
+    "\n",
+    "  \"edges\": [\n",
+    "    [ 1, 1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 7, 7, 7, 7, 9, 9, 9, 9, 9, 10, 10, 11, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 16, 17, 17, 18, 18, 18, 18, 18, 20, 21, 21, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 26, 26, 26, 27, 27, 27, 28, 28, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 32, 33, 33, 34, 34, 34, 34, 34, 34, 34, 35, 35, 36, 36, 36, 37, 37, 38, 38, 38, 38, 38, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 44, 44, 45, 46, 46, 47, 47, 47, 48, 49, 49, 49, 51, 52, 55 ],\n",
+    "    [ 5, 9, 11, 19, 7, 10, 6, 12, 18, 20, 28, 11, 12, 13, 19, 8, 10, 13, 16, 17, 11, 17, 19, 23, 29, 16, 22, 12, 17, 19, 31, 32, 35, 25, 29, 50, 17, 21, 22, 29, 19, 29, 20, 28, 33, 55, 56, 55, 29, 50, 29, 34, 36, 37, 39, 27, 30, 31, 44, 47, 48, 55, 56, 26, 29, 29, 42, 43, 31, 32, 55, 33, 45, 34, 43, 50, 38, 42, 44, 45, 56, 32, 35, 46, 47, 35, 45, 56, 39, 40, 42, 43, 51, 52, 54, 37, 46, 37, 39, 41, 41, 46, 42, 44, 49, 51, 54, 40, 41, 41, 49, 52, 46, 49, 53, 43, 51, 48, 49, 56, 47, 53, 48, 49, 53, 49, 52, 53, 54, 54, 54, 56 ]\n",
+    "  ],\n",
+    "\n",
+    "  \"tau\": 0.4853,\n",
+    "```"
    ]
   },
   {
@@ -131,22 +146,62 @@
     "bym2_fit.summary()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bym2_fit.diagnose()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Fit disconnected graphs on Scotland Lip cancer dataset with BYM2 model implemented in Stan, following Freni-Sterrantino\n"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "All components of the BYM2 model remain the same.  The same 2 $\\times$ J edges array encodes the spatial structure of the graph, and as before, the spatial prior `phi` is a vector, one element per node (areal unit).  In order to compute the elements of `phi` on a per-component basis, we use Stan's multi-indexing operators.  For each component, we provide a vector of indices into `phi` and into the edges_array.  Because Stan doesn't have ragged arrays (yet), we construct two square matrices, where the number of rows in each are the number of components (including singletons) in the graph, and the number of columns are the number of nodes and edges, respectively.  In addition we compute the per-component scaling factor.\n",
+    "\n",
+    "The helper function `index_components` takes an `nb` object as input and returns the list of data structures needed by the `BYM2_islands.stan` model.  We use the R package `jsonlite` to save this in JSON format.  For the Scotland map with 3 components, in file `scotland_3_comps_nbs.data.R`, function `index_components` and `write_json` produce the following input data, used in addition to the number of nodes, edges, and edge array:\n",
+    "\n",
+    "```\n",
+    "    \"K\":3,\n",
+    "    \"K_node_cts\":[53,2,1],\n",
+    "    \"K_edge_cts\":[126,1,0],\n",
+    "    \"K_node_idxs\":[[1,2,3,4,5,7,9,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,0,0,0],\n",
+    "\t\t   [6,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],\n",
+    "\t\t   [11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],\n",
+    "    \"K_edge_idxs\":[[1,2,3,4,5,6,7,8,9,10,11,12,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,0],\n",
+    "\t\t   [13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],\n",
+    "\t\t   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],\n",
+    "\n",
+    "  \"tau\": [0.4504, 0.25, 1],\n",
+    "```"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "from cmdstanpy import cmdstan_path, CmdStanModel\n",
-    "bym2_model = CmdStanModel(stan_file='bym2_islands.stan')\n",
-    "print(bym2_model.code())"
+    "bym2_islands_model = CmdStanModel(stan_file='bym2_islands.stan')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(bym2_islands_model.code())"
    ]
   },
   {
@@ -158,7 +213,7 @@
     "import json\n",
     "with open('scot_3_comp.data.json') as fd:\n",
     "    scot_data = json.load(fd)\n",
-    "print(scot_data)"
+    "# print(scot_data)"
    ]
   },
   {
@@ -167,16 +222,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bym2_fit = bym2_model.sample(data=scot_data)\n",
-    "bym2_fit.summary()"
+    "bym2_islands_fit = bym2_islands_model.sample(data=scot_data, max_treedepth=11)\n",
+    "bym2_islands_fit.summary()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "bym2_islands_fit.diagnose()"
+   ]
   },
   {
    "cell_type": "code",
diff --git a/knitr/car-iar-poisson/update_2021_02/bym2.stan b/knitr/car-iar-poisson/update_2021_02/bym2.stan
@@ -60,13 +60,12 @@ parameters {
   real beta;       // covariate
 
   // spatial effects
-  real logit_rho;  // proportion of spatial effect that's spatially smoothed
+  real<lower=0, upper=1> rho; // proportion of spatial effect that's spatially smoothed
   real<lower = 0> sigma;  // scale of spatial effects
   vector[I] theta;  // standardized heterogeneous spatial effects
   vector[I] phi;  // standardized spatially smoothed spatial effects
 }
 transformed parameters {
-  real<lower=0, upper=1> rho = inv_logit(logit_rho);
   // spatial effects (combine heterogeneous and spatially smoothed)
   vector[I] gamma = (sqrt(1 - rho) * theta + sqrt(rho) * sqrt(1 / tau) * phi) * sigma;
 }
@@ -78,11 +77,12 @@ model {
 
   // spatial hyperpriors and priors
   sigma ~ normal(0, 1);
-  logit_rho ~ normal(0, 1);
+  rho ~ beta(0.5, 0.5);
   theta ~ normal(0, 1);
   phi ~ standard_icar(edges);
 }
 generated quantities {
+  // posterior predictive checks
   vector[I] eta = log_E + alpha + x * beta + gamma * sigma;
   vector[I] y_prime = exp(eta);
   int y_rep[I,10];
@@ -97,4 +97,5 @@ generated quantities {
         y_rep[i,j] = poisson_log_rng(eta[i]);
     }
   }
+  real logit_rho = log(rho / (1.0 - rho));
 }
diff --git a/knitr/car-iar-poisson/update_2021_02/bym2_helpers.R b/knitr/car-iar-poisson/update_2021_02/bym2_helpers.R
@@ -96,10 +96,10 @@ index_components <- function(nb_obj) {
     }   
 
     return(list("K"=num_comps,
-                "K_num_nodes"=as.vector(table(comp_idxs)),
-                "K_num_edges"=edges_per_comp,
+                "K_node_cts"=as.vector(table(comp_idxs)),
+                "K_edge_cts"=edges_per_comp,
                 "K_node_idxs"=comp_node_idxs,
                 "K_edge_idxs"=comp_edge_idxs,
-                "K_scaling_factors"=scaling_factors))
+                "tau"=scaling_factors))
 }
 
diff --git a/knitr/car-iar-poisson/update_2021_02/bym2_islands.stan b/knitr/car-iar-poisson/update_2021_02/bym2_islands.stan
@@ -4,62 +4,59 @@ functions {
    * coefficients under the ICAR model with unit variance, where
    * adjacency is determined by the adjacency array and the spatial
    * structure is a disconnected graph which has at least one
-   * connected component.  The spatial structure is described by
-   * a 2-D adjacency array over the all edges in the areal map
-   * and a series of arrays of per-component indexes.
-   * The node arrays provide per-component masks into phi, the
-   * per-node ICAR component, and corresponding masks into the
-   * edge array.  Because the Stan language lacks ragged arrays,
-   * these are all zero-padded square matrices, and additional vectors
-   * record the number of nodes and edges in each component.
-   *
-   *
-   * Each connected component has a soft sum-to-zero constraint.
-   * Singleton components don't contribute to the ICAR model.
+   * connected component.  Each connected component has a
+   * soft sum-to-zero constraint.
+
+   * The spatial structure is described by a 2-D adjacency array
+   * over the all edges in the areal map and a arrays of the
+   * indices of per-component nodes and edges which are used as
+   * masks into phi and the adjacency matrix.   Because the Stan
+   * language lacks ragged arrays, these are all square matrices,
+   * padded out with zeros; additional vectors record the number
+   * of nodes and edges in each component.
    *
    * @param phi vector of varying effects
    * @param adjacency parallel arrays of indexes of adjacent elements of phi
-   * @param num_nodes array of sizes of per-component nodes
-   * @param num_edges array of sizes of per-component edges
+   * @param node_cts array of sizes of per-component nodes
+   * @param edge_cts array of sizes of per-component edges
    * @param node_idxs array of arrays of per_component node indexes.
    * @param edge_idxs array of arrays of per_component edge indexes.
    *
    * @return ICAR log probability density
    *
    * @reject if the the adjacency matrix does not have two rows
    * @reject if size mismatch between indexing arrays
-   * @reject if size mismatch between phi and dimension 2 of comp_members
+   * @reject if size mismatch between phi and node indexes columns.
    */
   real standard_icar_disconnected_lpdf(vector phi,
 				       int[ , ] adjacency,
-				       int[ ] num_nodes,
-				       int[ ] num_edges,
+				       int[ ] node_cts,
+				       int[ ] edge_cts,
 				       int[ , ] node_idxs,
 				       int[ , ] edge_idxs) {
     if (size(adjacency) != 2)
       reject("require 2 rows for adjacency array;",
              " found rows = ", size(adjacency));
       
-    if (!(size(num_nodes) == size(num_edges)
-	  && size(num_nodes) == size(node_idxs)
-	  && size(num_edges) == size(edge_idxs)))
+    if (!(size(phi) == dims(node_idxs)[2]
+	  && size(node_cts) == size(edge_cts)
+	  && size(node_cts) == size(node_idxs)
+	  && size(edge_cts) == size(edge_idxs)))
       reject("bad graph indexes, expecting ",
-	     size(num_nodes),
+	     size(node_cts),
 	     " rows for node and edge index matrices;",
              " inputs have ",
-	     size(num_nodes),
+	     size(node_cts),
 	     " and ",
-	     size(num_edges),
+	     size(edge_cts),
 	     " respectively.");
 
     real total = 0;
-    for (n in 1:size(num_nodes)) {
-      if (num_nodes[n] > 1)
-	total += -0.5 * dot_self(phi[adjacency[1, edge_idxs[n, 1:num_edges[n]]]] -
-				 phi[adjacency[2, edge_idxs[n, 1:num_edges[n]]]])
-	  + normal_lpdf(sum(phi[node_idxs[n, 1:num_nodes[n]]]) | 0, 0.001 * num_nodes[n]);
-      else
-	total += normal_lpdf(phi[node_idxs[n, 1]] | 0, inv_sqrt(1.0 * size(num_nodes)));
+    for (n in 1:size(node_cts)) {
+      if (node_cts[n] > 1)
+	total += -0.5 * dot_self(phi[adjacency[1, edge_idxs[n, 1:edge_cts[n]]]] -
+				 phi[adjacency[2, edge_idxs[n, 1:edge_cts[n]]]])
+	  + normal_lpdf(sum(phi[node_idxs[n, 1:node_cts[n]]]) | 0, 0.001 * node_cts[n]);
     }
     return total;
   }
@@ -71,8 +68,8 @@ data {
   int<lower = 1, upper = I> edges[2, J];  // node[1, j] adjacent to node[2, j]
 
   int<lower=0, upper=I> K;  // number of components in spatial graph
-  int<lower=0, upper=I> K_num_nodes[K];   // per-component nodes
-  int<lower=0, upper=J> K_num_edges[K];   // per-component edges
+  int<lower=0, upper=I> K_node_cts[K];   // per-component nodes
+  int<lower=0, upper=J> K_edge_cts[K];   // per-component edges
   int<lower=0, upper=I> K_node_idxs[K, I];  // rows contain per-component node indexes
   int<lower=0, upper=J> K_edge_idxs[K, J];  // rows contain per-component edge indexes
 
@@ -90,23 +87,24 @@ parameters {
   real beta;       // covariates
 
   // spatial effects
+  real<lower=0, upper=1> rho; // proportion unstructured vs. spatially structured variance
   real<lower = 0> sigma;  // scale of spatial effects
-  real<lower = 0, upper = 1> rho;  // proportion of spatial effect that's spatially smoothed
   vector[I] theta;  // standardized heterogeneous spatial effects
   vector[I] phi;  // standardized spatially smoothed spatial effects
 }
 transformed parameters {
   vector[I] gamma;
   // each component has its own spatial smoothing
+  // singleton nodes have distribution normal(0, 1/sqrt(K))
   for (k in 1:K) {
-    if (K_num_nodes[k] == 1) {
+    if (K_node_cts[k] == 1) {
       gamma[K_node_idxs[k,1]] =
 	theta[K_node_idxs[k,1]] + normal_lpdf(phi[K_node_idxs[k,1]] | 0, inv_sqrt(K));
     } else {
-      gamma[K_node_idxs[k, 1:K_num_nodes[k]]] = 
-	    (sqrt(1 - rho) * theta[K_node_idxs[k, 1:K_num_nodes[k]]]
+      gamma[K_node_idxs[k, 1:K_node_cts[k]]] = 
+	    (sqrt(1 - rho) * theta[K_node_idxs[k, 1:K_node_cts[k]]]
 	     + (sqrt(rho) * sqrt(1 / tau[k])
-		* phi[K_node_idxs[k, 1:K_num_nodes[k]]]) * sigma);
+		* phi[K_node_idxs[k, 1:K_node_cts[k]]]) * sigma);
     }
   }
 }
@@ -120,5 +118,23 @@ model {
   sigma ~ normal(0, 1);
   rho ~ beta(0.5, 0.5);
   theta ~ normal(0, 1);
-  phi ~ standard_icar_disconnected(edges, K_num_nodes, K_num_edges, K_node_idxs, K_edge_idxs);
+  phi ~ standard_icar_disconnected(edges, K_node_cts, K_edge_cts, K_node_idxs, K_edge_idxs);
+}
+generated quantities {
+  // posterior predictive checks
+  vector[I] eta = log_E + alpha + x * beta + gamma * sigma;
+  vector[I] y_prime = exp(eta);
+  int y_rep[I,10];
+  for (j in 1:10) {
+    if (max(eta) > 20) {
+      // avoid overflow in poisson_log_rng
+      print("max eta too big: ", max(eta));  
+      for (i in 1:I)
+	y_rep[i,j] = -1;
+    } else {
+      for (i in 1:I)
+        y_rep[i,j] = poisson_log_rng(eta[i]);
+    }
+  }
+  real logit_rho = log(rho / (1.0 - rho));
 }
diff --git a/knitr/car-iar-poisson/update_2021_02/scot_3_comp.data.json b/knitr/car-iar-poisson/update_2021_02/scot_3_comp.data.json
@@ -4,8 +4,8 @@
     "edges":[[1,1,1,2,2,3,4,4,4,5,5,5,6,7,7,7,7,9,9,9,9,10,10,13,13,14,14,14,15,15,15,16,16,16,16,17,17,18,18,18,18,18,20,21,21,23,23,23,23,23,24,24,24,24,24,24,24,24,25,25,26,26,26,27,27,27,28,28,29,29,29,30,30,30,30,30,31,31,31,31,32,33,33,34,34,34,34,34,34,34,35,35,36,36,36,37,37,38,38,38,38,38,39,39,40,40,40,41,41,41,42,42,44,44,45,46,46,47,47,47,48,49,49,49,51,52,55],
 	     [5,9,19,7,10,12,18,20,28,12,13,19,8,10,13,16,17,17,19,23,29,16,22,17,19,31,32,35,25,29,50,17,21,22,29,19,29,20,28,33,55,56,55,29,50,29,34,36,37,39,27,30,31,44,47,48,55,56,26,29,29,42,43,31,32,55,33,45,34,43,50,38,42,44,45,56,32,35,46,47,35,45,56,39,40,42,43,51,52,54,37,46,37,39,41,41,46,42,44,49,51,54,40,41,41,49,52,46,49,53,43,51,48,49,56,47,53,48,49,53,49,52,53,54,54,54,56]],
     "K":3,
-    "K_num_nodes":[53,2,1],
-    "K_num_edges":[126,1,0],
+    "K_node_cts":[53,2,1],
+    "K_edge_cts":[126,1,0],
     "K_node_idxs":[[1,2,3,4,5,7,9,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,0,0,0],
 		   [6,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 		   [11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],
diff --git a/knitr/car-iar-poisson/update_2021_02/scot_islands.data.json b/knitr/car-iar-poisson/update_2021_02/scot_islands.data.json
@@ -1,19 +1,21 @@
 {
-  "I": 56,
-  "J": 126,
+    "I": 56,
+    "J": 126,
 
-  "edges" :[[1,1,1,2,2,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,10,10,11,11,11,12,12,12,13,13,13,13,14,14,15,15,15,15,15,17,18,18,20,20,20,20,20,21,21,21,21,21,21,21,21,22,22,23,23,23,24,24,24,25,25,26,26,26,27,27,27,27,27,28,28,28,28,29,30,30,31,31,31,31,31,31,31,32,32,33,33,33,34,34,35,35,35,35,35,36,36,37,37,37,38,38,38,39,39,41,41,42,43,43,44,44,44,45,46,46,46,48,49,52],[5,7,16,6,8,9,15,17,25,9,10,16,8,10,13,14,14,16,20,26,13,19,14,16,28,29,32,22,26,47,14,18,19,26,16,26,17,25,30,52,53,52,26,47,26,31,33,34,36,24,27,28,41,44,45,52,53,23,26,26,39,40,28,29,52,30,42,31,40,47,35,39,41,42,53,29,32,43,44,32,42,53,36,37,39,40,48,49,51,34,43,34,36,38,38,43,39,41,46,48,51,37,38,38,46,49,43,46,50,40,48,45,46,53,44,50,45,46,50,46,49,50,51,51,51,53]], 
+    "edges" :[[1,1,1,2,2,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,10,10,11,11,11,12,12,12,13,13,13,13,14,14,15,15,15,15,15,17,18,18,20,20,20,20,20,21,21,21,21,21,21,21,21,22,22,23,23,23,24,24,24,25,25,26,26,26,27,27,27,27,27,28,28,28,28,29,30,30,31,31,31,31,31,31,31,32,32,33,33,33,34,34,35,35,35,35,35,36,36,37,37,37,38,38,38,39,39,41,41,42,43,43,44,44,44,45,46,46,46,48,49,52],[5,7,16,6,8,9,15,17,25,9,10,16,8,10,13,14,14,16,20,26,13,19,14,16,28,29,32,22,26,47,14,18,19,26,16,26,17,25,30,52,53,52,26,47,26,31,33,34,36,24,27,28,41,44,45,52,53,23,26,26,39,40,28,29,52,30,42,31,40,47,35,39,41,42,53,29,32,43,44,32,42,53,36,37,39,40,48,49,51,34,43,34,36,38,38,43,39,41,46,48,51,37,38,38,46,49,43,46,50,40,48,45,46,53,44,50,45,46,50,46,49,50,51,51,51,53]], 
 
     "K": 4,   
-    "K_size" : [53, 1, 1, 1 ],
-    "K_members" : [
-	[1,2,3,4,5,7,9,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,0,0,0],
-	[6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
-	[8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
-	[11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
-    ],
-	
-    "tau": [0.4504, 0.25, 1],
+    "K_node_cts":[53,1,1,1],
+    "K_edge_cts":[126,0,0,0],
+    "K_node_idxs":[[1,2,3,4,5,7,9,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,0,0,0],
+		   [6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+		   [8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+		   [11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],
+    "K_edge_idxs":[[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126],
+		   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+		   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+		   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],
+    "tau":[0.4504,1,1,1],
 
     "y": [ 9, 39, 11, 9, 15, 8, 26, 7, 6, 20, 13, 5, 3, 8, 17, 9, 2, 7, 9, 7, 16, 31, 11, 7, 19, 15, 7, 10, 16, 11, 5, 3, 7, 8, 11, 9, 11, 8, 6, 4, 10, 8, 2, 6, 19, 3, 2, 3, 28, 6, 1, 1, 1, 1, 0, 0  ],
     

Original file line number	Diff line number	Diff line change
`@@ -60,13 +60,12 @@ parameters {`
`60`	`60`	`real beta; // covariate`
`61`	`61`
`62`	`62`	`// spatial effects`
`63`		`- real logit_rho; // proportion of spatial effect that's spatially smoothed`
	`63`	`+ real<lower=0, upper=1> rho; // proportion of spatial effect that's spatially smoothed`
`64`	`64`	`real<lower = 0> sigma; // scale of spatial effects`
`65`	`65`	`vector[I] theta; // standardized heterogeneous spatial effects`
`66`	`66`	`vector[I] phi; // standardized spatially smoothed spatial effects`
`67`	`67`	`}`
`68`	`68`	`transformed parameters {`
`69`		`- real<lower=0, upper=1> rho = inv_logit(logit_rho);`
`70`	`69`	`// spatial effects (combine heterogeneous and spatially smoothed)`
`71`	`70`	`vector[I] gamma = (sqrt(1 - rho) * theta + sqrt(rho) * sqrt(1 / tau) * phi) * sigma;`
`72`	`71`	`}`
`@@ -78,11 +77,12 @@ model {`
`78`	`77`
`79`	`78`	`// spatial hyperpriors and priors`
`80`	`79`	`sigma ~ normal(0, 1);`
`81`		`- logit_rho ~ normal(0, 1);`
	`80`	`+ rho ~ beta(0.5, 0.5);`
`82`	`81`	`theta ~ normal(0, 1);`
`83`	`82`	`phi ~ standard_icar(edges);`
`84`	`83`	`}`
`85`	`84`	`generated quantities {`
	`85`	`+ // posterior predictive checks`
`86`	`86`	`vector[I] eta = log_E + alpha + x * beta + gamma * sigma;`
`87`	`87`	`vector[I] y_prime = exp(eta);`
`88`	`88`	`int y_rep[I,10];`
`@@ -97,4 +97,5 @@ generated quantities {`
`97`	`97`	`y_rep[i,j] = poisson_log_rng(eta[i]);`
`98`	`98`	`}`
`99`	`99`	`}`
	`100`	`+ real logit_rho = log(rho / (1.0 - rho));`
`100`	`101`	`}`
Original file line number	Diff line number	Diff line change
`@@ -96,10 +96,10 @@ index_components <- function(nb_obj) {`
`96`	`96`	`}`
`97`	`97`
`98`	`98`	`return(list("K"=num_comps,`
`99`		`- "K_num_nodes"=as.vector(table(comp_idxs)),`
`100`		`- "K_num_edges"=edges_per_comp,`
	`99`	`+ "K_node_cts"=as.vector(table(comp_idxs)),`
	`100`	`+ "K_edge_cts"=edges_per_comp,`
`101`	`101`	`"K_node_idxs"=comp_node_idxs,`
`102`	`102`	`"K_edge_idxs"=comp_edge_idxs,`
`103`		`- "K_scaling_factors"=scaling_factors))`
	`103`	`+ "tau"=scaling_factors))`
`104`	`104`	`}`
`105`	`105`