-
Notifications
You must be signed in to change notification settings - Fork 1
/
icar-functions.stan
172 lines (164 loc) · 6.47 KB
/
icar-functions.stan
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
/**
* Log probability of the intrinsic conditional autoregressive (ICAR) prior,
* excluding additive constants.
*
* @param phi Vector of parameters for spatial smoothing (on unit scale)
* @param spatial_scale Scale parameter for the ICAR model
* @param node1
* @param node2
* @param k number of groups
* @param group_size number of observational units in each group
* @param group_idx index of observations in order of their group membership
* @param has_theta If the model contains an independent partial pooling term, phi for singletons can be zeroed out; otherwise, they require a standard normal prior. Both BYM and BYM2 have theta.
*
* @return Log probability density of ICAR prior up to additive constant
**/
real icar_normal_lpdf(vector phi, real spatial_scale,
int[] node1, int[] node2,
int k, int[] group_size, int[] group_idx,
int has_theta) {
real lp;
int pos=1;
lp = -0.5 * dot_self(phi[node1] - phi[node2]);
if (has_theta) {
for (j in 1:k) {
/* sum to zero constraint for each connected group; singletons zero out */
lp += normal_lpdf(sum(phi[segment(group_idx, pos, group_size[j])]) | 0, 0.001 * group_size[j]);
pos += group_size[j];
}
} else {
/* does not have theta */
for (j in 1:k) {
if (group_size[j] > 1) {
/* same as above for non-singletons: sum to zero constraint */
lp += normal_lpdf(sum(phi[segment(group_idx, pos, group_size[j])]) | 0, 0.001 * group_size[j]);
} else {
/* its a singleton: independent Gaussian prior on phi */
lp += normal_lpdf(phi[ segment(group_idx, pos, group_size[j]) ] | 0, spatial_scale);
}
pos += group_size[j];
}
}
return lp;
}
/**
* Combine local and global partial-pooling components into the convolved BYM term.
*
* @param phi local (spatially autocorrelated) component (not phi_tilde!)
* @param theta global component (not theta_tilde!)
* @param n number of spatial units
* @param k number of connected groups
* @param group_size number of observational units in each group
* @param group_idx index of observations in order of their group membership
*
* @return BYM convolution vector
*/
vector convolve_bym(vector phi, vector theta,
int n, int k,
int[] group_size, int[] group_idx
) {
vector[n] convolution;
int pos=1;
for (j in 1:k) {
if (group_size[j] == 1) {
convolution[ segment(group_idx, pos, group_size[j]) ] = theta[ segment(group_idx, pos, group_size[j]) ];
} else {
convolution[ segment(group_idx, pos, group_size[j]) ] =
phi[ segment(group_idx, pos, group_size[j]) ] + theta[ segment(group_idx, pos, group_size[j]) ];
}
pos += group_size[j];
}
return convolution;
}
/**
* Combine local and global partial-pooling components into the convolved BYM2 term.
*
* @param phi_tilde local (spatially autocorrelated) component
* @param theta_tilde global component
* @param spatial_scale scale parameter for the convolution term
* @param n number of spatial units
* @param k number of connected groups
* @param group_size number of observational units in each group
* @param group_idx index of observations in order of their group membership
* @param inv_sqrt_scale_factor The scaling factor for the ICAR variance (see scale_c R function, using R-INLA);
* transformed from 1/scale^2 --> scale. Or, a vector of ones.
* @param rho proportion of convolution that is spatially autocorrelated
*
* @return BYM2 convolution vector
*/
vector convolve_bym2(vector phi_tilde, vector theta_tilde,
real spatial_scale,
int n, int k,
int[] group_size, int[] group_idx,
real rho, vector inv_sqrt_scale_factor
) {
vector[n] convolution;
int pos=1;
for (j in 1:k) {
if (group_size[j] == 1) {
convolution[ segment(group_idx, pos, group_size[j]) ] = spatial_scale * theta_tilde[ segment(group_idx, pos, group_size[j]) ];
} else {
convolution[ segment(group_idx, pos, group_size[j]) ] = spatial_scale * (
sqrt(rho) * inv_sqrt_scale_factor[j] * phi_tilde[ segment(group_idx, pos, group_size[j]) ] +
sqrt(1 - rho) * theta_tilde[ segment(group_idx, pos, group_size[j]) ]
);
}
pos += group_size[j];
}
return convolution;
}
/**
* Create phi from phi_tilde, inv_sqrt_scale_factor, and spatial_scale.
*
* @param phi_tilde local component (spatially autocorrelated)
* @param phi_scale scale parameter for ICAR Normal model
* @param inv_sqrt_scale_factor The scaling factor for the ICAR variance (see scale_c R function, using R-INLA);
* transformed from 1/scale^2 --> scale. Or, a vector of ones.
* @param n number of spatial units
* @param k number of connected groups
* @param group_size number of observational units in each group
* @param group_idx index of observations in order of their group membership
*
* @return phi vector of spatially autocorrelated coefficients
*/
vector make_phi(vector phi_tilde, real phi_scale,
vector inv_sqrt_scale_factor,
int n, int k,
int[] group_size, int[] group_idx
) {
vector[n] phi;
int pos=1;
for (j in 1:k) {
phi[ segment(group_idx, pos, group_size[j]) ] = phi_scale * inv_sqrt_scale_factor[j] * phi_tilde[ segment(group_idx, pos, group_size[j]) ];
pos += group_size[j];
}
return phi;
}
/**
* Create phi from phi_tilde, inv_sqrt_scale_factor, and spatial_scale:
* with separate scale parameters per component.
*
* @param phi_tilde local component (spatially autocorrelated)
* @param phi_scale scale parameters for ICAR Normal model, per connected component
* @param inv_sqrt_scale_factor The scaling factor for the ICAR variance (see scale_c R function, using R-INLA);
* transformed from 1/scale^2 --> scale. Or, a vector of ones.
* @param n number of spatial units
* @param k number of connected groups
* @param group_size number of observational units in each group
* @param group_idx index of observations in order of their group membership
*
* @return phi vector of spatially autocorrelated coefficients
*/
vector make_phi2(vector phi_tilde, vector phi_scale,
vector inv_sqrt_scale_factor,
int n, int k,
int[] group_size, int[] group_idx
) {
vector[n] phi;
int pos=1;
for (j in 1:k) {
phi[ segment(group_idx, pos, group_size[j]) ] = phi_scale[j] * inv_sqrt_scale_factor[j] * phi_tilde[ segment(group_idx, pos, group_size[j]) ];
pos += group_size[j];
}
return phi;
}