Skip to content

Outline of approach to fitting counts with scaling #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
seabbs wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Outline of approach to fitting counts with scaling #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6b7566c
add outline of count model
seabbs Jan 17, 2022
46a0ae4
updated indexing
seabbs Jan 17, 2022
e1596ad
better indexing
seabbs Jan 17, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion scripts/simple-example.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ fit <- incidence(
),
prob_detect = prob_detect, parallel_chains = 2, iter_warmup = 200,
chains = 2, model = mod, adapt_delta = 0.85, max_treedepth = 15,
data_args = list(gp_tune_model = tune, horizon = 14, differencing = 0),
data_args = list(gp_tune_model = tune, horizon = 14, differencing = 1),
keep_fit = TRUE
)
fit
Expand Down
10 changes: 5 additions & 5 deletions stan/functions/rt.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
// Code from:
// https://github.com/epiforecasts/EpiNow2/tree/master/inst/stan/functions
// discretised truncated gamma pmf
vector discretised_gamma_pmf(int[] y, real mu, real sigma, int max_val) {
vector discretised_gamma_pmf(int[] y, real mu, real sigma, int max_val, int off) {
int n = num_elements(y);
vector[n] pmf;
real trunc_pmf;
Expand All @@ -18,9 +18,9 @@ vector discretised_gamma_pmf(int[] y, real mu, real sigma, int max_val) {
beta = beta < small ? small : beta;
beta = beta > large ? large : beta;
// calculate pmf
trunc_pmf = gamma_cdf(max_val + 1 | alpha, beta) - gamma_cdf(1 | alpha, beta);
trunc_pmf = gamma_cdf(max_val + 1 | alpha, beta) - gamma_cdf(off | alpha, beta);
for (i in 1:n){
pmf[i] = (gamma_cdf(y[i] + 1 | alpha, beta) - gamma_cdf(y[i] | alpha, beta)) /
pmf[i] = (gamma_cdf(y[i] | alpha, beta) - gamma_cdf(y[i] - 1 | alpha, beta)) /
trunc_pmf;
}
return(pmf);
Expand Down Expand Up @@ -48,9 +48,9 @@ vector calculate_Rt(vector infections, int seeding_time,
vector[ot] infectiousness = rep_vector(1e-5, ot);
// calculate PMF of the generation time
for (i in 1:(max_gt)) {
gt_indexes[i] = max_gt - i + 1;
gt_indexes[i] = max_gt - i + 2;
}
gt_pmf = discretised_gamma_pmf(gt_indexes, gt_mean, gt_sd, max_gt);
gt_pmf = discretised_gamma_pmf(gt_indexes, gt_mean, gt_sd, max_gt, 1);
// calculate Rt using Cori et al. approach
for (s in 1:ot) {
infectiousness[s] += update_infectiousness(infections, gt_pmf, seeding_time,
Expand Down
56 changes: 49 additions & 7 deletions stan/inc2prev_antibodies.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@ data {
int prev_etime[obs]; // end times of positivity prevalence observations
int ab_stime[ab_obs]; // starting times of antibody prevalence observations
int ab_etime[ab_obs]; // end times of antibody prevalence observations
int c_obs; // Number of count observations
int c_grps; // Number of groups of counts
int c_obs_stime[c_grps]; // Group starting membership for each count
int c_obs_etime[c_grps]; // Group ending membership for each count
int c_time[c_obs]; // Occurance times for counts
vector[c_obs] counts;
vector[t] vacc; // vaccinations
int pbt; // maximum detection time
vector[pbt] prob_detect_mean; // at each time since infection, probability of detection
Expand Down Expand Up @@ -47,6 +53,7 @@ data {
vector[lvacc_ab_delay] vacc_ab_delay;
int prev_likelihood; // Should the likelihood for prevalence data be included
int ab_likelihood; // Should the likelihood for antibody data be included
int count_likelihood; // Should the likelihood for count data be included
}

transformed data {
Expand All @@ -73,7 +80,8 @@ parameters {
}

transformed parameters {
vector[t] gp; // value of gp at time t + initialisation
vector[t] inf_gp; // value of gp at time t + initialisation
vector[t] cscale[c_grps]; //time-varying scaling factor of infections to counts
vector[t] infections; // incident infections at time t
vector[t] infs_with_potential_abs; // Infections with the potential to have ab
vector<lower = 0, upper = 1>[t] dcases; // detectable cases at time t
Expand All @@ -83,18 +91,18 @@ transformed parameters {
vector[obs] combined_sigma;
vector[ab_obs] combined_ab_sigma;
// update gaussian process
gp[(1 + diff_order):t] = update_gp(PHI, M, L, alpha, rho, eta, 0);
inf_gp[(1 + diff_order):t] = update_gp(PHI, M, L, alpha[1], rho[1], eta[1:M],
0);
// setup differencing of the GP
if (diff_order) {
gp[1:diff_order] = init_growth;
inf_gp[1:diff_order] = init_growth;
for (i in 1:diff_order) {
gp = cumulative_sum(gp);
inf_gp = cumulative_sum(inf_gp);
}
}
// relative probability of infection
// inc_init is the mean incidence
infections = inv_logit(init_inc + gp);

infections = inv_logit(init_inc + inf_gp);
// calculate detectable cases
dcases = convolve(infections, prob_detect);
// calculate observed detectable cases
Expand All @@ -109,12 +117,26 @@ transformed parameters {
//combined standard error
combined_sigma = sqrt(square(sigma) + prev_sd2);
combined_ab_sigma = sqrt(square(ab_sigma) + ab_sd2);
// calculate infections reported as counts
if (c_grps) {
for (i in 1:c_grps) {
cscale[i] = update_gp(cPHI[c_grps], M, L, alpha[1+c_grps],
rho[1+c_grps],
eta[((c_grps-1)*M + 1):(c_grps*M)], 0);
cscale[i] = inv_logit(cscale_mean[i] + cscale[i]);
cdist = discretised_gamma_pmf(cdist_indexes, cdist_mean[i], cdist_sd[i],
cdist_max, 0);
ecounts[c_grps] = convolve(cscale[c_grps] .* infections, cdist[c_grps]);
}
}
}

model {
// gaussian process priors
rho ~ inv_gamma(lengthscale_alpha, lengthscale_beta);
alpha ~ std_normal() T[0,];
for (i in 1:gps) {
alpha[i] ~ std_normal() T[0,];
}
eta ~ std_normal();

// Initial infections
Expand All @@ -133,6 +155,16 @@ model {
logit(gamma) ~ normal(pgamma_mean, pgamma_sd);
logit(delta) ~ normal(pdelta, pdelta);


// Priors for count data model
if (c_grps) {
cscale_mean ~ normal(0, 5);
for (i in c_grps) {
cdist_mean[i] ~ normal(5, 10) T[0,];
cdist_sd[i] ~ normal(5, 10) T[0,];
}
}

sigma ~ normal(0.005, 0.0025) T[0,];
ab_sigma ~ normal(0.025, 0.025) T[0,];
if (prev_likelihood) {
Expand All @@ -141,6 +173,11 @@ model {
if (ab_likelihood) {
ab ~ normal(odab, combined_ab_sigma);
}
if (count_likelihood & c_grps) {
for (i in 1:c_grps) {
counts[c_obs_stime[i]:c_obs_etime[i]] ~ neg_binomial_2(ecounts[i], phi[i]);
}
}
}

generated quantities {
Expand All @@ -151,6 +188,11 @@ generated quantities {
// sample estimated prevalence
est_prev = normal_rng(odcases, combined_sigma);
est_ab = normal_rng(odab, combined_ab_sigma);
if (c_grps) {
for (i in 1:c_grps) {
est_counts[c_obs_stime[i]:c_obs_etime[i]] ~ neg_binomial_2(ecounts[i], phi[i]);
}
}
// sample generation time
real gtm_sample = normal_rng(gtm[1], gtm[2]);
real gtsd_sample = normal_rng(gtsd[1], gtsd[2]);
Expand Down