checkpoint work on bym2 et al

Author: mitzimorris

knitr/car-iar-poisson/update_2021_02/BYM2 islands.ipynb

@@ -385,7 +385,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -407,9 +407,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "num nodes: 56, num edges: 126\n",
+      "num components: 4\n",
+      "scaling factors: [0.4504, 1, 1, 1]\n",
+      "nodes per component: [53, 1, 1, 1]\n",
+      "node indices: [[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]]\n"
+     ]
+    }
+   ],
    "source": [
     "with open('scotland_islands.data.json') as fd:\n",
     "    islands_data = json.load(fd)\n",

knitr/car-iar-poisson/update_2021_02/bym2.stan

@@ -67,10 +67,11 @@ parameters {
 }
 transformed parameters {
   // spatial effects (combine heterogeneous and spatially smoothed)
-  vector[I] gamma = (sqrt(1 - rho) * theta + sqrt(rho / tau) * phi) * sigma;
+  vector[I] gamma = sigma * (sqrt(1 - rho) * theta
+			     + sqrt(rho / tau) * phi);
 }
 model {
-  y ~ poisson_log(log_E + alpha + x * beta + gamma * sigma);  // co-variates
+  y ~ poisson_log(log_E + alpha + x * beta + gamma);  // likelihood
 
   alpha ~ normal(0, 1);
   beta ~ normal(0, 1);
@@ -82,20 +83,6 @@ model {
   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];
-  //   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));
+  vector[I] eta = log_E + alpha + x * beta + gamma;
+  vector[I] y_est = exp(eta);
 }

knitr/car-iar-poisson/update_2021_02/bym2_islands.stan

@@ -52,15 +52,16 @@ functions {
 	     " elements in phi and columns of node_idxs, and ",
 	     num_edges,
 	     " columns of edge_idxs.");
-
     real total = 0;
     for (n in 1:num_comps) {
-      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]);
-      else
-	total += normal_lpdf(phi[node_idxs[n, 1]] | 0, 1);
+      int nodes[node_cts[n]] = node_idxs[n, 1:node_cts[n]];
+      if (node_cts[n] > 1) {
+	int edges[edge_cts[n]] = edge_idxs[n, 1:edge_cts[n]];
+	total += -0.5 * dot_self(phi[adjacency[1, edges]] - phi[adjacency[2, edges]])
+	  + normal_lpdf(sum(phi[nodes]) | 0, 0.001 * node_cts[n]);
+      } else {
+	total += normal_lpdf(phi[nodes] | 0, 1);
+      }
     }
     return total;
   }
@@ -77,36 +78,37 @@ data {
   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
 
-  vector[K] tau; // scaling factor
+  vector[K] tau; // per-component scaling factor
 
   int<lower=0> y[I];              // count outcomes
   vector<lower=0>[I] E;           // exposure
   vector[I] x;                 // predictor
 }
 transformed data {
   vector[I] log_E = log(E);
+  vector[I] taus;
+  for (k in 1:K) {
+    int node_idxs[K_node_cts[k]] = K_node_idxs[k, 1:K_node_cts[k]];
+    for (j in 1:size(node_idxs))
+      taus[node_idxs[j]] = tau[k];
+  }
 }
 parameters {
-  real alpha;            // intercept
-  real beta;       // covariates
+  real alpha;      // intercept
+  real beta;       // single covariate (Scotland dataset)
 
   // 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> sigma;  // scale of spatial effects
   vector[I] theta;  // standardized heterogeneous spatial effects
   vector[I] phi;  // standardized spatially smoothed spatial effects
 }
 transformed parameters {
-  vector[I] gamma;
-  for (k in 1: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 / tau[k]) * phi[K_node_idxs[k, 1:K_node_cts[k]]])
-      * sigma;
-}
+  vector[I] gamma = sigma * (sqrt(1 - rho) * theta
+			     + sqrt(rho ./ taus) .* phi);
+}  
 model {
-  y ~ poisson_log(log_E + alpha + x * beta + gamma * sigma);  // co-variates
+  y ~ poisson_log(log_E + alpha + x * beta + gamma);  // likelihood
 
   alpha ~ normal(0, 1);
   beta ~ normal(0, 1);
@@ -118,20 +120,6 @@ model {
   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));
+  vector[I] eta = log_E + alpha + x * beta + gamma;
+  vector[I] y_est = exp(eta);
 }

knitr/car-iar-poisson/update_2021_02/scot_3_comp.png

@@ -1,18 +1,17 @@
 {
-  "I": 56,
-  "J": 132,
+    "I": 56,
+    "J": 132,
 
   "edges": [
     [ 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 ],
     [ 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 ]
   ],
 
-  "tau": 0.4853,
+    "tau": 0.4853,
 
-  "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  ],
-
-  "E": [ 1.4, 8.7, 3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5, 6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4, 7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7, 19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8  ],
-
-  "x": [1.6, 1.6, 1.0, 2.4, 1.0, 2.4, 1.0, 0.7, 0.7, 1.6, 0.7, 1.6, 1.0, 2.4, 0.7, 1.6, 1.0, 0.7, 0.7, 1.0, 0.7, 1.6, 1.0, 0.7, 0.1, 0.1, 0.7, 0.7, 1.0, 1.0, 0.7, 2.4, 1.0, 0.7, 0.7, 0, 1.0, 0.1, 1.6, 0, 0.1, 1.6, 1.6, 0, 0.1, 0.7, 0.1, 0.1, 0, 0.1, 0.1, 0, 0.1, 0.1, 1.6, 1.0]
+    "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  ],
+    
+    "E": [ 1.4, 8.7, 3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5, 6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4, 7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7, 19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8  ],
 
+    "x": [16, 16, 10, 24, 10, 24, 10, 7, 7, 16, 7, 16, 10, 24, 7, 16, 10, 7, 7, 10, 7, 16, 10, 7, 1, 1, 7, 7, 10, 10, 7, 24, 10, 7, 7, 0, 10, 1, 16, 0, 1, 16, 16, 0, 1, 7, 1, 1, 0, 1, 1, 0, 1, 1, 16, 10]
 }

knitr/car-iar-poisson/update_2021_02/scot_islands.png

@@ -1,4 +1,4 @@
-{
+xo{
   "I": 56,
   "J": 132,
 

knitr/car-iar-poisson/update_2021_02/scotland_connected.data.json

@@ -22,5 +22,5 @@
 
     "E": [ 1.4, 8.7, 3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5, 6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4, 7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7, 19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8  ],
 
-  "x": [1.6, 1.6, 1.0, 2.4, 1.0, 2.4, 1.0, 0.7, 0.7, 1.6, 0.7, 1.6, 1.0, 2.4, 0.7, 1.6, 1.0, 0.7, 0.7, 1.0, 0.7, 1.6, 1.0, 0.7, 0.1, 0.1, 0.7, 0.7, 1.0, 1.0, 0.7, 2.4, 1.0, 0.7, 0.7, 0, 1.0, 0.1, 1.6, 0, 0.1, 1.6, 1.6, 0, 0.1, 0.7, 0.1, 0.1, 0, 0.1, 0.1, 0, 0.1, 0.1, 1.6, 1.0]
+    "x": [16, 16, 10, 24, 10, 24, 10, 7, 7, 16, 7, 16, 10, 24, 7, 16, 10, 7, 7, 10, 7, 16, 10, 7, 1, 1, 7, 7, 10, 10, 7, 24, 10, 7, 7, 0, 10, 1, 16, 0, 1, 16, 16, 0, 1, 7, 1, 1, 0, 1, 1, 0, 1, 1, 16, 10]
 }