diff --git a/scripts/simple-example.R b/scripts/simple-example.R
index b7c6a62..f6224ad 100644
--- a/scripts/simple-example.R
+++ b/scripts/simple-example.R
@@ -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
diff --git a/stan/functions/rt.stan b/stan/functions/rt.stan
index d456f3d..03f1b72 100644
--- a/stan/functions/rt.stan
+++ b/stan/functions/rt.stan
@@ -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;
@@ -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);
@@ -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,
diff --git a/stan/inc2prev_antibodies.stan b/stan/inc2prev_antibodies.stan
index f1214de..9a22bf6 100644
--- a/stan/inc2prev_antibodies.stan
+++ b/stan/inc2prev_antibodies.stan
@@ -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
@@ -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 {
@@ -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
@@ -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
@@ -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
@@ -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) {
@@ -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 {
@@ -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]);