Thank you!
Here are my three models:
The model proposed at the beginning of this post (applied to age-period-gender-area data),There are a total of 14,900 parameters here:
functions {
real partial_sum_lpdf(array[] real logm,
int start, int end,
array[] int death) {
return poisson_log_lupmf(death[start:end]| logm);
}
/**
* Return the log probability of a proper conditional autoregressive (CAR) prior
* with a sparse representation for the adjacency matrix
*
* @param phi Vector containing the parameters with a CAR prior
* @param tau Precision parameter for the CAR prior (real)
* @param alpha Dependence (usually spatial) parameter for the CAR prior (real)
* @param W_sparse Sparse representation of adjacency matrix (int array)
* @param n Length of phi (int)
* @param W_n Number of adjacent pairs (int)
* @param D_sparse Number of neighbors for each location (vector)
* @param lambda Eigenvalues of D^{-1/2}*W*D^{-1/2} (vector)
*
* @return Log probability density of CAR prior up to additive constant
*/
real sparse_car_lpdf(vector phi, real tau, real alpha,
array[,] int W_sparse, vector D_sparse, vector lambda, int n, int W_n) {
row_vector[n] phit_D; // phi' * D
row_vector[n] phit_W; // phi' * W
vector[n] ldet_terms;
phit_D = (phi .* D_sparse)';
phit_W = rep_row_vector(0, n);
for (i in 1:W_n) {
phit_W[W_sparse[i, 1]] = phit_W[W_sparse[i, 1]] + phi[W_sparse[i, 2]];
phit_W[W_sparse[i, 2]] = phit_W[W_sparse[i, 2]] + phi[W_sparse[i, 1]];
}
for (i in 1:n) ldet_terms[i] = log1m(alpha * lambda[i]);
return 0.5 * (n * log(tau)
+ sum(ldet_terms)
- tau * (phit_D * phi - alpha * (phit_W * phi)));
}
}
data {
int<lower=0> T; // number of years = 7
int<lower=0> A; // number of age categories = 21
int<lower=0> P; // number of area categories = 46
int<lower=0> Tf; // number of forecast years = 2
array[A*T*P*2] int<lower=0> death; // deaths
array[A*T*P*2] real<lower=0> expos; // exposures
int n; // numver of iszero = 1
matrix<lower = 0, upper = 1>[(P-n),(P-n)] W;
array[n] int<lower=1> iszero;
int W_n; // number of adjacent region pairs
}
transformed data {
array[A*T*P*2] real loge = log(expos); // log exposures
int AT = A*T;
int ATG = AT*2;
int ATGP = ATG*P;
array[W_n, 2] int W_sparse; // adjacency pairs
vector[P-n] D_sparse; // diagonal of D (number of neigbors for each site)
vector[P-n] lambda; // eigenvalues of invsqrtD * W * invsqrtD
{ // generate sparse representation for W
int counter;
counter = 1;
// loop over upper triangular part of W to identify neighbor pairs
for (i in 1:(P-n-1)) {
for (j in (i + 1):(P-n)) {
if (W[i, j] == 1) {
W_sparse[counter, 1] = i;
W_sparse[counter, 2] = j;
counter = counter + 1;
}
}
}
}
for (i in 1:(P-n)) D_sparse[i] = sum(W[i]);
{
vector[(P-n)] invsqrtD;
for (i in 1:(P-n)) {
invsqrtD[i] = 1 / sqrt(D_sparse[i]);
}
lambda = eigenvalues_sym(quad_form(W, diag_matrix(invsqrtD)));
}
}
parameters {
array[2,P] vector[A] agp; // alpha_x
simplex[A] b1; // beta_x, strictly positive and sums to 1
simplex[A] b2;
array[2,P] simplex[A] bgp1;
array[2,P] simplex[A] bgp2;
vector[T] ks; // vector of k_t differences
array[2,P] vector[T] kgps;
vector[T+A-2] gs; // vector of k_t differences
array[2,P] vector[T+A-2] ggps;
vector[T] alpha_k;
vector[T] delta_k;
vector[T+A-2] alpha_g;
vector[T+A-2] delta_g;
array[2,P] vector[T] alpha_kgp;
array[2,P] vector[T] delta_kgp;
array[2,P] vector[T+A-2] alpha_ggp;
array[2,P] vector[T+A-2] delta_ggp;
real<lower = 0> sigma_k; // standard deviation of the random walk
array[2] vector<lower = 0>[P] sigma_kgp;
real<lower = 0> sigma_g; // standard deviation of the random walk
array[2] vector<lower = 0>[P] sigma_ggp;
real<lower=0> sigmak_a;
real<lower=0> sigmak_d;
real<lower=0> sigmag_a;
real<lower=0> sigmag_d;
array[2] vector<lower = 0>[P] sigmakgp_a;
array[2] vector<lower = 0>[P] sigmakgp_d;
array[2] vector<lower = 0>[P] sigmaggp_a;
array[2] vector<lower = 0>[P] sigmaggp_d;
vector[P-n] thetas; // spatial effects without zero
real<lower = 0> tau;
real<lower = 0, upper = 1> alpha;
}
transformed parameters {
vector[T] k;
array[2,P] vector[T] kgp;
real km;
array[2,P] real kgpm;
vector[T+A-2] g;
array[2,P] vector[T+A-2] ggp;
real gm;
array[2,P] real ggpm;
array[A*T*P*2] real logm;
vector[P] theta;// spatial effects with zero
km = mean(ks);
k = ks - km;
gm = mean(gs);
g = gs - gm;
for(s in 1:2) for(p in 1:P){
kgpm[s,p] = mean(kgps[s,p]);
kgp[s,p] = kgps[s,p] - kgpm[s,p];
ggpm[s,p] = mean(ggps[s,p]);
ggp[s,p] = ggps[s,p] - ggpm[s,p];
}
theta[1:46] = thetas;
theta[47] = 0;
for (p in 1:P) for (s in 1:2) for (j in 1:T) for (i in 1:A) {
if(i == 1){
logm[ATG*(p-1) + AT*(s-1) + A*(j-1) + i] = loge[ATG*(p-1) + AT*(s-1) + A*(j-1) + i] +
b1[i] * k[j] + b2[i] * g[j-i+A-1] +
agp[s,p][i] + bgp1[s,p][i] * kgp[s,p][j] + bgp2[s,p][i] * ggp[s,p][j-i+A-1] +
theta[p];
}
else{
logm[ATG*(p-1) + AT*(s-1) + A*(j-1) + i] = loge[ATG*(p-1) + AT*(s-1) + A*(j-1) + i] +
b1[i] * k[j] + b2[i] * g[j-i+A] +
agp[s,p][i] + bgp1[s,p][i] * kgp[s,p][j] + bgp2[s,p][i] * ggp[s,p][j-i+A] +
theta[p];
}
}
}
model {
int grainsize = ATGP%/%4;
target += reduce_sum(partial_sum_lupdf, logm, grainsize, death);
// k[1:2], g[1:2] ~ improper flat priors; or whatever
target += normal_lupdf(ks|alpha_k, sigma_k);
target += normal_lupdf(alpha_k[2:T]|alpha_k[1:(T-1)] + delta_k[2:T],sigmak_a);
target += normal_lupdf(delta_k[2:T]|delta_k[1:(T-1)], sigmak_d);
target += normal_lupdf(gs|alpha_g,sigma_g);
target += normal_lupdf(alpha_g[2:(T+A-2)]|alpha_g[1:(T+A-3)] + delta_g[2:(T+A-2)], sigmag_a);
target += normal_lupdf(delta_g[2:(T+A-2)]|delta_g[1:(T+A-3)], sigmag_d);
target += dirichlet_lupdf(b1|rep_vector(1, A));
target += dirichlet_lupdf(b2|rep_vector(1, A));
target += student_t_lupdf(sigma_k|4, 0, 1);
target += student_t_lupdf(sigma_g|4, 0, 1);
target += student_t_lupdf(sigmak_a|4, 0, 1);
target += student_t_lupdf(sigmak_d|4, 0, 1);
target += student_t_lupdf(sigmag_a|4, 0, 1);
target += student_t_lupdf(sigmag_d|4, 0, 1);
for(s in 1:2) for(p in 1:P){
target += normal_lupdf(kgps[s,p]|alpha_kgp[s,p], sigma_kgp[s][p]);
target += normal_lupdf(alpha_kgp[s,p][2:T]|alpha_kgp[s,p][1:(T-1)] + delta_kgp[s,p][2:T],sigmakgp_a[s][p]);
target += normal_lupdf(delta_kgp[s,p][2:T]|delta_kgp[s,p][1:(T-1)], sigmakgp_d[s][p]);
target += normal_lupdf(ggps[s,p]|alpha_ggp[s,p],sigma_ggp[s][p]);
target += normal_lupdf(alpha_ggp[s,p][2:(T+A-2)]|alpha_ggp[s,p][1:(T+A-3)] + delta_ggp[s,p][2:(T+A-2)], sigmaggp_a[s][p]);
target += normal_lupdf(delta_ggp[s,p][2:(T+A-2)]|delta_ggp[s,p][1:(T+A-3)], sigmaggp_d[s][p]);
target += normal_lupdf(agp[s,p]|0,10);
target += dirichlet_lupdf(bgp1[s,p]|rep_vector(1, A));
target += dirichlet_lupdf(bgp2[s,p]|rep_vector(1, A));
}
for(s in 1:2){
target += student_t_lupdf(sigma_kgp[s]|4, 0, 1);
target += student_t_lupdf(sigma_ggp[s]|4, 0, 1);
target += student_t_lupdf(sigmakgp_a[s]|4, 0, 1);
target += student_t_lupdf(sigmakgp_d[s]|4, 0, 1);
target += student_t_lupdf(sigmaggp_a[s]|4, 0, 1);
target += student_t_lupdf(sigmaggp_d[s]|4, 0, 1);
}
target += sparse_car_lpdf(thetas| tau, alpha, W_sparse, D_sparse, lambda, P-n, W_n);
target += gamma_lupdf(tau|2,2);
}
generated quantities {
vector[Tf] k_p;
vector[Tf] alpha_k_p;
vector[Tf] delta_k_p;
array[2,P] vector[Tf] kgp_p;
array[2,P] vector[Tf] alpha_kgp_p;
array[2,P] vector[Tf] delta_kgp_p;
vector[Tf] g_p;
vector[Tf] alpha_g_p;
vector[Tf] delta_g_p;
array[2,P] vector[Tf] ggp_p;
array[2,P] vector[Tf] alpha_ggp_p;
array[2,P] vector[Tf] delta_ggp_p;
vector[T+Tf+A-2] g_f;
array[2,P] vector[T+Tf+A-2] ggp_f;
vector[A*Tf*2*P] mfore; // predicted death rates
vector[ATGP] log_lik;
array[ATGP] real m;
int pos1 = 1;
delta_k_p[1] = delta_k[T] + sigmak_d * normal_rng(0,1);
alpha_k_p[1] = alpha_k[T] + delta_k_p[1] + sigmak_a * normal_rng(0,1);
for (t in 2:Tf){
delta_k_p[t] = delta_k_p[t-1] + sigmak_d * normal_rng(0,1);
alpha_k_p[t] = alpha_k_p[t-1] + delta_k_p[t] + sigmak_a * normal_rng(0,1);
}
for (t in 1:Tf){k_p[t] = alpha_k_p[t] - km + sigma_k * normal_rng(0,1);}
delta_g_p[1] = delta_g[T+A-2] + sigmag_d * normal_rng(0,1);
alpha_g_p[1] = alpha_g[T+A-2] + delta_g_p[1] + sigmag_a * normal_rng(0,1);
for (t in 2:Tf){
delta_g_p[t] = delta_g_p[t-1] + sigmag_d * normal_rng(0,1);
alpha_g_p[t] = alpha_g_p[t-1] + delta_g_p[t] + sigmag_a * normal_rng(0,1);
}
for (t in 1:Tf){g_p[t] = alpha_g_p[t] - gm + sigma_g * normal_rng(0,1);}
g_f = append_row(g, g_p);
for(s in 1:2) for(p in 1:P){
delta_kgp_p[s,p][1] = delta_kgp[s,p][T] + sigmakgp_d[s][p] * normal_rng(0,1);
alpha_kgp_p[s,p][1] = alpha_kgp[s,p][T] + delta_kgp_p[s,p][1] + sigmakgp_a[s][p] * normal_rng(0,1);
for (t in 2:Tf){
delta_kgp_p[s,p][t] = delta_kgp_p[s,p][t-1] + sigmakgp_d[s,p] * normal_rng(0,1);
alpha_kgp_p[s,p][t] = alpha_kgp_p[s,p][t-1] + delta_kgp_p[s,p][t] + sigmakgp_a[s][p] * normal_rng(0,1);
}
for (t in 1:Tf){kgp_p[s,p][t] = alpha_kgp_p[s,p][t] - kgpm[s,p] + sigma_kgp[s][p] * normal_rng(0,1);}
delta_ggp_p[s,p][1] = delta_ggp[s,p][T+A-2] + sigmaggp_d[s][p] * normal_rng(0,1);
alpha_ggp_p[s,p][1] = alpha_ggp[s,p][T+A-2] + delta_ggp_p[s,p][1] + sigmaggp_a[s][p] * normal_rng(0,1);
for (t in 2:Tf){
delta_ggp_p[s,p][t] = delta_ggp_p[s,p][t-1] + sigmaggp_d[s][p] * normal_rng(0,1);
alpha_ggp_p[s,p][t] = alpha_ggp_p[s,p][t-1] + delta_ggp_p[s,p][t] + sigmaggp_a[s][p] * normal_rng(0,1);
}
for (t in 1:Tf){ggp_p[s,p][t] = alpha_ggp_p[s,p][t] - ggpm[s,p] + sigma_ggp[s][p] * normal_rng(0,1);}
ggp_f[s,p] = append_row(ggp[s,p], ggp_p[s,p]);
}
for (p in 1:P) for (s in 1:2) for (j in 1:Tf) for (i in 1:A) {
if(i == 1){
mfore[pos1] = exp(b1[i] * k_p[j] + b2[i] * g_f[T+j-i+A-1] +
agp[s,p][i] + bgp1[s,p][i] * kgp_p[s,p][j] + bgp2[s,p][i] * ggp_f[s,p][T+j-i+A-1] +
theta[p]);
}
else{
mfore[pos1] = exp(b1[i] * k_p[j] + b2[i] * g_f[T+j-i+A] +
agp[s,p][i] + bgp1[s,p][i] * kgp_p[s,p][j] + bgp2[s,p][i] * ggp_f[s,p][T+j-i+A] +
theta[p]);
}
pos1 += 1;
}
for(i in 1:ATGP){
m[i] = exp(logm[i] - loge[i]);
log_lik[i] = poisson_log_lpmf(death[i]|logm[i]);
}
}
The function sparse_car_lpdf in the model is referenced from this page:Exact sparse CAR models in Stan
model A(applied to age-period data),There are a total of 825 parameters here:
functions {
real partial_sum_lpdf(array[] real logm,
int start, int end,
array[] int death) {
return poisson_log_lupmf(death[start:end] | logm);
}
}
data {
int<lower=0> T; // number of years = 36
int<lower=0> A; // number of age categories = 101
int<lower=0> Tf; // number of forecast years = 10
array[A*T] int<lower=0> death; // deaths
array[A*T] real<lower=0> expos; // exposures
}
transformed data {
array[A*T] real loge = log(expos); // log exposures
}
parameters {
vector[A] a; // alpha_x
simplex[A] b1; // beta_x, strictly positive and sums to 1
simplex[A] b2;
vector[T] ks; // vector of k_t(without first element)
vector[T+A-1] gs;
vector[T] alpha_k;
vector[T] delta_k;
vector[T+A-1] alpha_g;
vector[T+A-1] delta_g;
real<lower = 0> sigma_k; // standard deviation of the random walk
real<lower = 0> sigma_g;
real<lower=0> sigmak_a;
real<lower=0> sigmak_d;
real<lower=0> sigmag_a;
real<lower=0> sigmag_d;
}
transformed parameters {
vector[T] k;
vector[T+A-1] g;
array[A*T] real logm;
real km;
real gm;
km = mean(ks);
k = ks - km;
gm = mean(gs);
g = gs - gm;
for (j in 1:T) for (i in 1:A) {
logm[A*(j-1)+i] = loge[A*(j-1)+i] + a[i] + b1[i] * k[j] + b2[i] * g[j-i+A];
}
}
model {
int grainsize = A*T%/%3;
target += reduce_sum(partial_sum_lpdf, logm, grainsize, death);
target += normal_lpdf(ks|alpha_k,sigma_k);
target += normal_lpdf(alpha_k[1]|delta_k[1],sigmak_a);
target += normal_lpdf(alpha_k[2:T]|alpha_k[1:(T-1)]+delta_k[2:T],sigmak_a);
target += normal_lpdf(delta_k[1]|0,sigmak_d);
target += normal_lpdf(delta_k[2:T]|delta_k[1:(T-1)],sigmak_d);
target += normal_lpdf(gs|alpha_g,sigma_g);
target += normal_lpdf(alpha_g[1]|delta_g[1],sigmag_a);
target += normal_lpdf(alpha_g[2:(T+A-1)]|alpha_g[1:(T+A-2)]+delta_g[2:(T+A-1)],sigmag_a);
target += normal_lpdf(delta_g[1]|0,sigmag_d);
target += normal_lpdf(delta_g[2:(T+A-1)]|delta_g[1:(T+A-2)],sigmag_d);
target += normal_lpdf(a|0,10);
target += dirichlet_lpdf(b1|rep_vector(1, A));
target += dirichlet_lpdf(b2|rep_vector(1, A));
target += student_t_lpdf(sigma_k|4, 0, 1);
target += student_t_lpdf(sigma_g|4, 0, 1);
target += student_t_lpdf(sigmak_a|4, 0, 1);
target += student_t_lpdf(sigmak_d|4, 0, 1);
target += student_t_lpdf(sigmag_a|4, 0, 1);
target += student_t_lpdf(sigmag_d|4, 0, 1);
}
generated quantities {
vector[Tf] k_p;
vector[Tf] alpha_k_p;
vector[Tf] delta_k_p;
vector[Tf] g_p;
vector[T+Tf+A-1] g_f;
vector[Tf] alpha_g_p;
vector[Tf] delta_g_p;
vector[A*Tf] mfore; // predicted death rates
vector[A*T] log_lik;
array[A*T] real m;
int pos1 = 1;
delta_k_p[1] = delta_k[T] + sigmak_d * normal_rng(0,1);
alpha_k_p[1] = alpha_k[T] + delta_k_p[1] + sigmak_a * normal_rng(0,1);
for (t in 2:Tf){
delta_k_p[t] = delta_k_p[t-1] + sigmak_d * normal_rng(0,1);
alpha_k_p[t] = alpha_k_p[t-1] + delta_k_p[t] + sigmak_a * normal_rng(0,1);
}
for (t in 1:Tf){k_p[t] = alpha_k_p[t] - km + sigma_k * normal_rng(0,1);}
delta_g_p[1] = delta_g[T+A-1] + sigmag_d * normal_rng(0,1);
alpha_g_p[1] = alpha_g[T+A-1] + delta_g_p[1] + sigmag_a * normal_rng(0,1);
for (t in 2:Tf){
delta_g_p[t] = delta_g_p[t-1] + sigmag_d * normal_rng(0,1);
alpha_g_p[t] = alpha_g_p[t-1] + delta_g_p[t] + sigmag_a * normal_rng(0,1);
}
for (t in 1:Tf){g_p[t] = alpha_g_p[t] - gm + sigma_g * normal_rng(0,1);}
g_f = append_row(g, g_p);
for (t in 1:Tf) for (i in 1:A) {
mfore[pos1] = exp(a[i] + b1[i] * k_p[t] + b2[i] * g_f[T+t-i+A]);
pos1 += 1;
}
for(i in 1:A*T){
m[i] = exp(logm[i] - loge[i]);
log_lik[i] = poisson_log_lpmf(death[i]|logm[i]);
}
}
model B(applied to age-period data),There are a total of 553 parameters here:
functions {
real partial_sum_lpdf(array[] real logm,
int start, int end,
array[] int death) {
return poisson_log_lupmf(death[start:end] | logm);
}
}
data {
int<lower=0> T; // number of years = 36
int<lower=0> A; // number of age categories = 101
int<lower=0> Tf; // number of forecast years = 10
array[A*T] int<lower=0> death; // deaths
array[A*T] real<lower=0> expos; // exposures
}
transformed data {
array[A*T] real loge = log(expos); // log exposures
}
parameters {
vector[A] a; // alpha_x
simplex[A] b1; // beta_x, strictly positive and sums to 1
simplex[A] b2;
vector[T] ks; // vector of k_t(without first element)
vector[T+A-1] gs;
real c;
real<lower = -1, upper = 1> psi;
vector[T] alpha_k;
vector[T] delta_k;
real<lower = 0> sigma_k; // standard deviation of the random walk
real<lower = 0> sigma_g;
real<lower=0> sigmak_a;
real<lower=0> sigmak_d;
}
transformed parameters {
vector[T] k;
vector[T+A-1] g;
array[A*T] real logm;
real km;
real gm;
km = mean(ks);
k = ks - km;
gm = mean(gs);
g = gs - gm;
for (j in 1:T) for (i in 1:A) {
logm[A*(j-1)+i] = loge[A*(j-1)+i] + a[i] + b1[i] * k[j] + b2[i] * g[j-i+A];
}
}
model {
int grainsize = A*T%/%3;
target += reduce_sum(partial_sum_lpdf, logm, grainsize, death);
target += normal_lpdf(ks|alpha_k,sigma_k);
target += normal_lpdf(alpha_k[1]|delta_k[1],sigmak_a);
target += normal_lpdf(alpha_k[2:T]|alpha_k[1:(T-1)]+delta_k[2:T],sigmak_a);
target += normal_lpdf(delta_k[1]|0,sigmak_d);
target += normal_lpdf(delta_k[2:T]|delta_k[1:(T-1)],sigmak_d);
target += normal_lpdf(gs[3:(T+A-1)]|c + gs[2:(T+A-2)] + psi * (gs[2:(T+A-2)] - gs[1:(T+A-3)] - c),sigma_g);
target += normal_lpdf(a|0,10);
target += normal_lpdf(c|0,10);
target += normal_lpdf(psi|0,1);
target += dirichlet_lpdf(b1|rep_vector(1, A));
target += dirichlet_lpdf(b2|rep_vector(1, A));
target += student_t_lpdf(sigma_k|4, 0, 1);
target += student_t_lpdf(sigma_g|4, 0, 1);
target += student_t_lpdf(sigmak_a|4, 0, 1);
target += student_t_lpdf(sigmak_d|4, 0, 1);
}
generated quantities {
vector[Tf] k_p;
vector[Tf] alpha_k_p;
vector[Tf] delta_k_p;
vector[Tf] g_p;
vector[T+Tf+A-1] g_f;
vector[A*Tf] mfore; // predicted death rates
vector[A*T] log_lik;
array[A*T] real m;
int pos1 = 1;
delta_k_p[1] = delta_k[T] + sigmak_d * normal_rng(0,1);
alpha_k_p[1] = alpha_k[T] + delta_k_p[1] + sigmak_a * normal_rng(0,1);
for (t in 2:Tf){
delta_k_p[t] = delta_k_p[t-1] + sigmak_d * normal_rng(0,1);
alpha_k_p[t] = alpha_k_p[t-1] + delta_k_p[t] + sigmak_a * normal_rng(0,1);
}
for (t in 1:Tf){k_p[t] = alpha_k_p[t] - km + sigma_k * normal_rng(0,1);}
g_p[1] = c + g[T+A-1] + psi * (g[T+A-1] - g[T+A-2] - c) + sigma_g * normal_rng(0,1);
g_p[2] = c + g_p[1] + psi * (g_p[1] - g[T+A-1] - c) + sigma_g * normal_rng(0,1);
for (t in 3:Tf){g_p[t] = c + g_p[t-1] + psi * (g_p[t-1] - g_p[t-2] - c) + sigma_g * normal_rng(0,1);}
g_f = append_row(g, g_p);
for (t in 1:Tf) for (i in 1:A) {
mfore[pos1] = exp(a[i] + b1[i] * k_p[t] + b2[i] * g_f[T+t-i+A]);
pos1 += 1;
}
for(i in 1:A*T){
m[i] = exp(logm[i] - loge[i]);
log_lik[i] = poisson_log_lpmf(death[i]|logm[i]);
}
}