Skip to content

Commit cc81a5a

Browse files
mitzimorrismitzimorris
committed
model works, questions remain
1 parent 7c75532 commit cc81a5a

13 files changed

+420
-158
lines changed

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

Lines changed: 317 additions & 80 deletions
Large diffs are not rendered by default.

‎knitr/car-iar-poisson/update_2021_02/BYM2_islands_notebook.pdf‎

679 KB
Binary file not shown.

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

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ parameters {
6767
}
6868
transformed parameters {
6969
// spatial effects (combine heterogeneous and spatially smoothed)
70-
vector[I] gamma = (sqrt(1 - rho) * theta + sqrt(rho) * sqrt(1 / tau) * phi) * sigma;
70+
vector[I] gamma = (sqrt(1 - rho) * theta + sqrt(rho / tau) * phi) * sigma;
7171
}
7272
model {
7373
y ~ poisson_log(log_E + alpha + x * beta + gamma * sigma); // co-variates
@@ -77,25 +77,25 @@ model {
7777

7878
// spatial hyperpriors and priors
7979
sigma ~ normal(0, 1);
80-
rho ~ beta(0.5, 0.5);
80+
rho ~ normal(0, 1);
8181
theta ~ normal(0, 1);
8282
phi ~ standard_icar(edges);
8383
}
8484
generated quantities {
8585
// posterior predictive checks
8686
vector[I] eta = log_E + alpha + x * beta + gamma * sigma;
8787
vector[I] y_prime = exp(eta);
88-
int y_rep[I,10];
89-
for (j in 1:10) {
90-
if (max(eta) > 20) {
91-
// avoid overflow in poisson_log_rng
92-
print("max eta too big: ", max(eta));
93-
for (i in 1:I)
94-
y_rep[i,j] = -1;
95-
} else {
96-
for (i in 1:I)
97-
y_rep[i,j] = poisson_log_rng(eta[i]);
98-
}
99-
}
88+
// int y_rep[I,10];
89+
// for (j in 1:10) {
90+
// if (max(eta) > 20) {
91+
// // avoid overflow in poisson_log_rng
92+
// print("max eta too big: ", max(eta));
93+
// for (i in 1:I)
94+
// y_rep[i,j] = -1;
95+
// } else {
96+
// for (i in 1:I)
97+
// y_rep[i,j] = poisson_log_rng(eta[i]);
98+
// }
99+
// }
100100
real logit_rho = log(rho / (1.0 - rho));
101101
}

‎knitr/car-iar-poisson/update_2021_02/bym2_helpers.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ index_components <- function(nb_obj) {
9191
edges = nb_to_edge_array(sub_nb)
9292
scaling_factors[k] = scaling_factor(edges)
9393
} else {
94-
scaling_factors[k] = 1.0
94+
scaling_factors[k] = 1.0/num_comps
9595
}
9696
}
9797

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

Lines changed: 36 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,9 @@ functions {
55
* adjacency is determined by the adjacency array and the spatial
66
* structure is a disconnected graph which has at least one
77
* connected component. Each connected component has a
8-
* soft sum-to-zero constraint.
9-
8+
* soft sum-to-zero constraint. Singleton nodes have
9+
* distribution normal(0, 1/sqrt(K))
10+
*
1011
* The spatial structure is described by a 2-D adjacency array
1112
* over the all edges in the areal map and a arrays of the
1213
* indices of per-component nodes and edges which are used as
@@ -34,29 +35,32 @@ functions {
3435
int[ ] edge_cts,
3536
int[ , ] node_idxs,
3637
int[ , ] edge_idxs) {
38+
int num_nodes = size(phi);
39+
int num_edges = dims(adjacency)[2];
40+
int num_comps = size(edge_cts);
3741
if (size(adjacency) != 2)
3842
reject("require 2 rows for adjacency array;",
3943
" found rows = ", size(adjacency));
40-
41-
if (!(size(phi) == dims(node_idxs)[2]
44+
if (!(num_nodes == dims(node_idxs)[2]
4245
&& size(node_cts) == size(edge_cts)
4346
&& size(node_cts) == size(node_idxs)
4447
&& size(edge_cts) == size(edge_idxs)))
45-
reject("bad graph indexes, expecting ",
46-
size(node_cts),
47-
" rows for node and edge index matrices;",
48-
" inputs have ",
49-
size(node_cts),
50-
" and ",
51-
size(edge_cts),
52-
" respectively.");
48+
reject("arguments have size mismatch, expecting ",
49+
num_comps,
50+
" rows for node_cts edge_cts, node_idxs, and edge_idxs,",
51+
num_nodes,
52+
" elements in phi and columns of node_idxs, and ",
53+
num_edges,
54+
" columns of edge_idxs.");
5355

5456
real total = 0;
55-
for (n in 1:size(node_cts)) {
57+
for (n in 1:num_comps) {
5658
if (node_cts[n] > 1)
5759
total += -0.5 * dot_self(phi[adjacency[1, edge_idxs[n, 1:edge_cts[n]]]] -
5860
phi[adjacency[2, edge_idxs[n, 1:edge_cts[n]]]])
5961
+ normal_lpdf(sum(phi[node_idxs[n, 1:node_cts[n]]]) | 0, 0.001 * node_cts[n]);
62+
else
63+
total += normal_lpdf(phi[node_idxs[n, 1]] | 0, 1);
6064
}
6165
return total;
6266
}
@@ -94,19 +98,12 @@ parameters {
9498
}
9599
transformed parameters {
96100
vector[I] gamma;
97-
// each component has its own spatial smoothing
98-
// singleton nodes have distribution normal(0, 1/sqrt(K))
99-
for (k in 1:K) {
100-
if (K_node_cts[k] == 1) {
101-
gamma[K_node_idxs[k,1]] =
102-
theta[K_node_idxs[k,1]] + normal_lpdf(phi[K_node_idxs[k,1]] | 0, inv_sqrt(K));
103-
} else {
104-
gamma[K_node_idxs[k, 1:K_node_cts[k]]] =
105-
(sqrt(1 - rho) * theta[K_node_idxs[k, 1:K_node_cts[k]]]
106-
+ (sqrt(rho) * sqrt(1 / tau[k])
107-
* phi[K_node_idxs[k, 1:K_node_cts[k]]]) * sigma);
108-
}
109-
}
101+
for (k in 1:K)
102+
gamma[K_node_idxs[k, 1:K_node_cts[k]]] =
103+
(sqrt(1 - rho) * theta[K_node_idxs[k, 1:K_node_cts[k]]]
104+
+
105+
sqrt(rho / tau[k]) * phi[K_node_idxs[k, 1:K_node_cts[k]]])
106+
* sigma;
110107
}
111108
model {
112109
y ~ poisson_log(log_E + alpha + x * beta + gamma * sigma); // co-variates
@@ -116,25 +113,25 @@ model {
116113

117114
// spatial hyperpriors and priors
118115
sigma ~ normal(0, 1);
119-
rho ~ beta(0.5, 0.5);
116+
rho ~ normal(0, 1);
120117
theta ~ normal(0, 1);
121118
phi ~ standard_icar_disconnected(edges, K_node_cts, K_edge_cts, K_node_idxs, K_edge_idxs);
122119
}
123120
generated quantities {
124121
// posterior predictive checks
125122
vector[I] eta = log_E + alpha + x * beta + gamma * sigma;
126123
vector[I] y_prime = exp(eta);
127-
int y_rep[I,10];
128-
for (j in 1:10) {
129-
if (max(eta) > 20) {
130-
// avoid overflow in poisson_log_rng
131-
print("max eta too big: ", max(eta));
132-
for (i in 1:I)
133-
y_rep[i,j] = -1;
134-
} else {
135-
for (i in 1:I)
136-
y_rep[i,j] = poisson_log_rng(eta[i]);
137-
}
138-
}
124+
// int y_rep[I,10];
125+
// for (j in 1:10) {
126+
// if (max(eta) > 20) {
127+
// // avoid overflow in poisson_log_rng
128+
// print("max eta too big: ", max(eta));
129+
// for (i in 1:I)
130+
// y_rep[i,j] = -1;
131+
// } else {
132+
// for (i in 1:I)
133+
// y_rep[i,j] = poisson_log_rng(eta[i]);
134+
// }
135+
// }
139136
real logit_rho = log(rho / (1.0 - rho));
140137
}

‎knitr/car-iar-poisson/update_2021_02/fig_1a_nbs.R‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
#spdep nb object list of integer vectors
2+
fig_1a_nb = list(
3+
as.vector(c(2,3), mode="integer"), # 1
4+
as.vector(c(1,3), mode="integer"), # 2
5+
as.vector(c(1,2,5,6), mode="integer"), # 3
6+
as.vector(c(5,6), mode="integer"), # 4
7+
as.vector(c(4,3), mode="integer"), # 5
8+
as.vector(c(4,3), mode="integer")) # 6
9+
10+
attr(fig_1a_nb, "class") <- "nb"
11+
12+
attr(fig_1a_nb, "region.id") <-
13+
c("1", "2", "3", "4", "5", "6")
14+

‎knitr/car-iar-poisson/update_2021_02/fig_1b_nbs.R‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
#spdep nb object list of integer vectors
2+
fig_1b_nb = list(
3+
as.vector(c(2,3), mode="integer"), # 1
4+
as.vector(c(1,3), mode="integer"), # 2
5+
as.vector(c(1,2), mode="integer"), # 3
6+
as.vector(c(5), mode="integer"), # 4
7+
as.vector(c(4), mode="integer"), # 5
8+
as.vector(c(0), mode="integer")) # 6
9+
10+
attr(fig_1b_nb, "class") <- "nb"
11+
12+
attr(fig_1b_nb, "region.id") <-
13+
c("1", "2", "3", "4", "5", "6")
14+

‎knitr/car-iar-poisson/update_2021_02/munge_scotland.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
library(SpatialEpi)
22
data("scotland")
33

4-
source("scotland_nbs.data.R")
4+
source("scotland_connected_nbs.data.R")
55
source("scotland_3_comp_nbs.data.R")
66
source("scotland_islands_nbs.data.R")
77

knitr/car-iar-poisson/update_2021_02/scot_3_comp.data.json renamed to knitr/car-iar-poisson/update_2021_02/scotland_3_comp.data.json

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,11 +13,11 @@
1313
[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],
1414
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],
1515

16-
"tau": [0.4504, 0.25, 1],
16+
"tau": [0.4504, 0.25, 0.3333],
1717

1818
"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 ],
1919

2020
"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 ],
2121

22-
"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 ]
22+
"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]
2323
}

knitr/car-iar-poisson/update_2021_02/scot_connected.data.json renamed to knitr/car-iar-poisson/update_2021_02/scotland_connected.data.json

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,5 +13,6 @@
1313

1414
"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 ],
1515

16-
"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 ]
16+
"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]
17+
1718
}

0 commit comments

Comments
 (0)