diff --git a/NAMESPACE b/NAMESPACE
index ba16632..bb3653b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -29,7 +29,6 @@ export(calculate_parameter_estimates)
 export(create_sample_linelist)
 export(default_model)
 export(detection)
-export(estimate_severity)
 export(fit)
 export(has_detection)
 export(has_outcome)
diff --git a/NEWS.md b/NEWS.md
index 9c13fa6..6a1ae37 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -125,6 +125,7 @@ fit_result
 - Added `summary` generic for `SeverityEstimateFit` objects. [#30](https://github.com/ACCIDDA/SeverityEstimate/issues/30).
 - Added an integration-style notebook under `inst/notebooks/` for synthetic parameter-recovery checks, along with local render commands for heavier non-vignette R Markdown analyses. [#116](https://github.com/ACCIDDA/SeverityEstimate/issues/116).
 - Added `calculate_hazard` to get the underlying force of infection as a `data.frame` for diagnostics/additional analysis. [#9](https://github.com/ACCIDDA/SeverityEstimate/issues/9).
+- Removed the `estimate_severity` API. [#88](https://github.com/ACCIDDA/SeverityEstimate/issues/88).
 
 # SeverityEstimate 0.0.1
 
diff --git a/R/estimate_severity.R b/R/estimate_severity.R
deleted file mode 100644
index 20d7df6..0000000
--- a/R/estimate_severity.R
+++ /dev/null
@@ -1,196 +0,0 @@
-#' @title
-#' Estimate Severity From A Line List
-#'
-#' @description
-#' Fit a severity estimate model accounting for under reporting of asymptomatic
-#' and mildly symptomatic cases.
-#'
-#' @param linelist A `data.frame` (or `data.frame` extending object like a
-#' `tibble`) of the line list data.
-#' @param population A `data.frame` (or `data.frame` extending object like a
-#' `tibble`) of the population data.
-#' @param time_period A character of columns describing the time period such as
-#' 'week' or 'day'
-#' @param strata A character of columns describing the attributes to stratify
-#' the data on.
-# TODO: `surveillance` and `outcome` are special, need to add those details
-#' @param surveillance A character of columns describing the surveillance
-#' methods. The values in this column must be coercible to 'Active' and
-#' 'Passive'.
-#' @param outcome A character of columns describing the outcome of the line list
-#' entry. The values in this column must be coercible to 'Asymptomatic',
-#' 'Death', and 'Symptomatic'.
-#' @param population_value A unit length character corresponding to the column
-#' in the `population` data.frame (or data.frame like object) for the population
-#' value.
-#' @param surveillance_reference Either `NULL`, a `data.frame`, or a vector of
-#' values to use as a reference for the surveillance dimension. If a vector then
-#' it is expected that `surveillance` is unit length.
-#' @param outcome_reference Either `NULL`, a `data.frame`, or a vector of
-#' values to use as a reference for the outcome dimension. If a vector then
-#' it is expected that `outcome` is unit length.
-#' @param time_period_reference Either `NULL`, a `data.frame`, or a vector of
-#' values to use as a reference for the time period dimension. If a vector then
-#' it is expected that `time_period` is unit length.
-#' @param strata_reference Either `NULL`, a `data.frame`, or a vector of
-#' values to use as a reference for the strata dimension. If a vector then
-#' it is expected that `strata` is unit length.
-#' @param hazard_std A single length numeric greater than zero for the standard
-#' deviation to use in the community hazard prior distribution.
-#' @param degrees_of_freedom A single length integer greater than zero for the
-#' degrees of freedom to use in the model.
-#' @param active_prior The parameters for the prior beta distribution for the
-#' active detection probability. Can be specified as 'alpha'/'beta',
-#' 'mean'/'var', 'mean'/'sd', or 'mean'/'concentration'.
-#' @param passive_asymptomatic_prior The parameters for the prior beta
-#' distribution for the passive asymptomatic detection probability. Can be
-#' specified as 'alpha'/'beta', 'mean'/'var', 'mean'/'sd', or
-#' 'mean'/'concentration'.
-#' @param passive_symptomatic_prior The parameters for the prior beta
-#' distribution for the passive symptomatic detection probability. Can be
-#' specified as 'alpha'/'beta', 'mean'/'var', 'mean'/'sd', or
-#' 'mean'/'concentration'.
-#' @param ... Further optional args that are eventually given to
-#' [rstan::sampling()] related to fitting.
-#'
-#' @returns
-#' A \linkS4class{SeverityEstimateFit} S4 object.
-#'
-#' @importFrom checkmate assert_integerish
-#' @importFrom checkmate assert_number
-#' @importFrom checkmate assert_string
-#' @importFrom methods new
-#' @importFrom utils tail
-#' @export
-estimate_severity <- function(
-  linelist,
-  population,
-  surveillance,
-  outcome,
-  time_period = character(),
-  strata = character(),
-  population_value = utils::tail(names(population), n = 1L),
-  surveillance_reference = NULL,
-  outcome_reference = NULL,
-  time_period_reference = NULL,
-  strata_reference = NULL,
-  hazard_std = 3.0,
-  degrees_of_freedom = 1L,
-  active_prior = c("alpha" = 1.0, "beta" = 25.0),
-  passive_asymptomatic_prior = c("alpha" = 1.0, "beta" = 2.0),
-  passive_symptomatic_prior = c("alpha" = 25.0, "beta" = 1.0),
-  ...
-) {
-  # Input validation
-  checkmate::assert_string(surveillance)
-  checkmate::assert_string(outcome)
-  linelist <- is_data_frame(
-    linelist,
-    has_string_columns = c(surveillance, outcome)
-  )
-  population <- is_data_frame(population)
-  # Check on hazard variance
-  checkmate::assert_number(
-    hazard_std,
-    lower = .Machine$double.eps,
-    finite = TRUE
-  )
-  # Check on degrees of freedom
-  checkmate::assert_integerish(degrees_of_freedom, lower = 0L)
-  degrees_of_freedom <- as.integer(degrees_of_freedom)
-  # Checks on prior parameters
-  active_prior <- beta_parameterization(active_prior)
-  passive_asymptomatic_prior <- beta_parameterization(
-    passive_asymptomatic_prior
-  )
-  passive_symptomatic_prior <- beta_parameterization(passive_symptomatic_prior)
-  # Construct the incidence array
-  arrays <- incidence_population_arrays(
-    linelist,
-    population,
-    time_period,
-    strata,
-    surveillance,
-    outcome,
-    population_value,
-    time_period_reference,
-    strata_reference,
-    surveillance_reference,
-    outcome_reference
-  )
-  surveillance_df <- format_surveillance_data_frame(arrays$surveillance)
-  outcome_df <- format_outcome_data_frame(arrays$outcome)
-  # Compile together the data given and format it
-  surveillance_ind <- match(
-    c("Active", "Passive", "Unknown"),
-    surveillance_df[, surveillance]
-  )
-  outcome_ind <- match(
-    c("Asymptomatic", "Symptomatic", "Death"),
-    outcome_df[, outcome]
-  )
-  incidence_without_outcome <- rowSums(arrays$incidence, dims = 3L)
-  active_ind <- which(
-    arrays$linelist_ind[, "surveillance", drop = TRUE] == surveillance_ind[1L]
-  )
-  passive_ind <- which(
-    arrays$linelist_ind[, "surveillance", drop = TRUE] == surveillance_ind[2L]
-  )
-  data <- list(
-    strata_groups = nrow(arrays$strata),
-    time_groups = nrow(arrays$time_period),
-    I_passive = incidence_without_outcome[,, surveillance_ind[2L]],
-    I_active = incidence_without_outcome[,, surveillance_ind[1L]],
-    population = arrays$population,
-    observed_active = length(active_ind),
-    observed_passive = length(passive_ind),
-    strata_active = arrays$linelist_ind[active_ind, "strata", drop = TRUE],
-    symptoms_active = (arrays$linelist_ind[
-      active_ind,
-      "outcome",
-      drop = TRUE
-    ] %in%
-      outcome_ind[2L:3L]),
-    dead_active = (arrays$linelist_ind[active_ind, "outcome", drop = TRUE] %in%
-      outcome_ind[3L]),
-    strata_passive = arrays$linelist_ind[passive_ind, "strata", drop = TRUE],
-    symptoms_passive = (arrays$linelist_ind[
-      passive_ind,
-      "outcome",
-      drop = TRUE
-    ] %in%
-      outcome_ind[2L:3L]),
-    dead_passive = (arrays$linelist_ind[
-      passive_ind,
-      "outcome",
-      drop = TRUE
-    ] %in%
-      outcome_ind[3L]),
-    hazard_std = hazard_std,
-    degrees_of_freedom = degrees_of_freedom,
-    active_detection_alpha = active_prior["alpha"],
-    active_detection_beta = active_prior["beta"],
-    passive_asymptomatic_alpha = passive_asymptomatic_prior["alpha"],
-    passive_asymptomatic_beta = passive_asymptomatic_prior["beta"],
-    passive_symptomatic_alpha = passive_symptomatic_prior["alpha"],
-    passive_symptomatic_beta = passive_symptomatic_prior["beta"]
-  )
-  # Pass along everything to the model
-  model_fit <- stan_model(
-    "severity_estimate_by_surveillance",
-    data = data,
-    ...
-  )
-  # Create and return severity fit
-  severity_estimate_fit <- new(
-    "SeverityEstimateFit",
-    model_fit = model_fit,
-    population = arrays$population,
-    incidence = arrays$incidence,
-    time_period = arrays$time_period,
-    strata = arrays$strata,
-    surveillance = surveillance_df,
-    outcome = outcome_df
-  )
-  severity_estimate_fit
-}
diff --git a/R/stanmodels.R b/R/stanmodels.R
index 358c82e..cbdf916 100644
--- a/R/stanmodels.R
+++ b/R/stanmodels.R
@@ -1,11 +1,10 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("estimate_severity", "severity_estimate_by_surveillance")
+stanmodels <- c("estimate_severity")
 
 # load each stan module
 Rcpp::loadModule("stan_fit4estimate_severity_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4severity_estimate_by_surveillance_mod", what = TRUE)
 
 # instantiate each stanmodel object
 stanmodels <- sapply(stanmodels, function(model_name) {
diff --git a/inst/stan/severity_estimate_by_surveillance.stan b/inst/stan/severity_estimate_by_surveillance.stan
deleted file mode 100644
index 877080f..0000000
--- a/inst/stan/severity_estimate_by_surveillance.stan
+++ /dev/null
@@ -1,170 +0,0 @@
-data {
-  // *Dimensions of data*
-  // The number of strata
-  int <lower=0> strata_groups;
-  // The number of times
-  int <lower=0> time_groups;
-  // The number of cases observed through active surveillance
-  int <lower=0> observed_active;
-  // The number of cases observed through passive surveillance
-  int <lower=0> observed_passive;
-  // *Matrices of data*
-  // The number of incidence detected through active surveillance
-  array[time_groups, strata_groups] int<lower=0> I_active;
-  // The number of incidence detected through passive surveillance
-  array[time_groups, strata_groups] int<lower=0> I_passive;
-  // The total population
-  array[strata_groups] int<lower=0> population;
-  // *Vectors of data*
-  // The strata number of the actively observed cases
-  array[observed_active] int<lower=0> strata_active;
-  // Indicator if the actively observed case presented with symptoms
-  array[observed_active] int<lower=0> symptoms_active;
-  // Indicator if the actively observed case died
-  array[observed_active] int<lower=0> dead_active;
-  // The strata number of the passively observed cases
-  array[observed_passive] int<lower=0> strata_passive;
-  // Indicator if the passively observed case presented with symptoms
-  array[observed_passive] int<lower=0> symptoms_passive;
-  // Indicator if the passively observed case died
-  array[observed_passive] int<lower=0> dead_passive;
-  // *Model parameters and priors*
-  // The stdev of the community hazard
-  real <lower=0> hazard_std;
-  // Spline degrees of freedom for mortality and symptom terms
-  int <lower=1> degrees_of_freedom;
-  // Active detection probability prior
-  real <lower=0> active_detection_alpha;
-  real <lower=0> active_detection_beta;
-  // Passive detection probability prior
-  real <lower=0> passive_asymptomatic_alpha;
-  real <lower=0> passive_asymptomatic_beta;
-  real <lower=0> passive_symptomatic_alpha;
-  real <lower=0> passive_symptomatic_beta;
-}
-parameters {
-   // Symptom development/mortality spline coefficients
-  array[1 + degrees_of_freedom] real alpha;
-  array[1 + degrees_of_freedom] real mort_coef;
-  // The hazard of infection in each time step
-  array[time_groups, strata_groups] real logit_hzd;
-  // Active detection probability
-  real <lower=0, upper=1> active_detection;
-  // Passive detection probabilities
-  real <lower=0, upper=1> passive_asymptomatic_detection;
-  real <lower=0, upper=1> passive_symptomatic_detection;
-}
-transformed parameters {
-  // *Transformed parameters*
-  // Strata specific symptom/mortality rate
-  array[strata_groups] real<lower=0, upper=1> xi;
-  array[strata_groups] real<lower=0, upper=1> mortality;
-  // The susceptibles/casesat at each time by strata
-  array[time_groups, strata_groups] real<lower=0> S;
-  array[time_groups, strata_groups] real<lower=0> C;
-  // Intermediates
-  array[strata_groups] real<lower=0> theta;
-  array[strata_groups] real<lower=0> passive_denom;
-  real xi_tmp;
-  real mort_tmp;
-  // Calculate IFR/SIR for each strata group
-  for (i in 1:strata_groups) {
-    // Calculate the symptomatic/mortality rate from the spline
-    xi_tmp = alpha[1];
-    mort_tmp = mort_coef[1];
-    for (j in 1:degrees_of_freedom) {
-      xi_tmp = xi_tmp + (alpha[j + 1] * pow(i - 1.0, j));
-      mort_tmp = mort_tmp + (mort_coef[j + 1] * pow(i - 1.0, j));
-    }
-    xi[i] = inv_logit(xi_tmp);
-    mortality[i] = inv_logit(mort_tmp);
-    // For first time step assume the population at risk is the full population
-    S[1, i] = population[i];
-    C[1, i] = population[i] * inv_logit(logit_hzd[1, i]);
-  }
-  // For each subsequent time step we assume the number of passive cases is
-  // based on the passive hazard and the number susceptible
-  for (i in 2:time_groups) {
-    for (j in 1:strata_groups) {
-      S[i, j] = S[i - 1, j] - C[i - 1, j];
-      C[i, j] = S[i, j] * inv_logit(logit_hzd[i, j]);
-    }
-  }
-  // Helpful commonly reused variable
-  for (i in 1:strata_groups) {
-    theta[i] = (passive_asymptomatic_detection * (1.0 - xi[i]))
-      + (passive_symptomatic_detection * xi[i]);
-    passive_denom[i] = 1.0 - ((1.0 - mortality[i]) * (1.0 - theta[i]));
-  }
-}
-model {
-  // *Model priors*
-  // Spline coefficients
-  for (i in 1:(degrees_of_freedom + 1)) {
-    alpha[i] ~ normal(0.0, 10000.0);
-    mort_coef[i] ~ normal(0.0, 10000.0);
-  }
-  // Priors for detection probabilities
-  active_detection ~ beta(active_detection_alpha, active_detection_beta);
-  passive_asymptomatic_detection ~ beta(
-    passive_asymptomatic_alpha, passive_asymptomatic_beta
-  );
-  passive_symptomatic_detection ~ beta(
-    passive_symptomatic_alpha, passive_symptomatic_beta
-  );
-  // Prior for community hazard
-  for (i in 1:time_groups) {
-    for (j in 1:strata_groups) {
-      logit_hzd[i, j] ~ normal(
-        logit((1.0 * I_passive[i, j])/population[j]
-          + (100.0 * machine_precision())),
-        hazard_std
-      );
-      I_active[i, j] ~ poisson(active_detection * C[i, j]);
-      I_passive[i, j] ~ poisson((1.0 - active_detection) * theta[j] * C[i, j]);
-    }
-  }
-  // Symtpomatic proability in active cases
-  for (i in 1:observed_active) {
-    symptoms_active[i] ~ bernoulli(xi[strata_active[i]]);
-    dead_active[i] ~ bernoulli(mortality[strata_active[i]]);
-  }
-  // Reporting delay in active cases
-  for (i in 1:observed_passive) {
-    symptoms_passive[i] ~ bernoulli(
-      (1.0 - ((1.0 - mortality[strata_passive[i]]))
-        * (1.0 - (passive_symptomatic_detection * xi[strata_passive[i]])))
-          / passive_denom[strata_passive[i]]
-		);
-		dead_passive[i] ~ bernoulli(
-		  mortality[strata_passive[i]] / passive_denom[strata_passive[i]]
-		);
-  }
-}
-generated quantities {
-  // *Additional helpful quantities*
-  // Distribution of unseen additional active/passive cases
-  array[time_groups, strata_groups] int<lower=0> C_active_additional;
-  array[time_groups, strata_groups] int<lower=0> C_passive_additional;
-  real gq_tmp;
-  // Loop over each time step and generate cases
-  for (i in 1:time_groups) {
-    for (j in 1:strata_groups) {
-      // First generate the active cases, used in passive cases
-      gq_tmp = C[i, j] - I_active[i, j];
-      if (gq_tmp > 0) {
-        C_active_additional[i, j] = poisson_rng(active_detection * gq_tmp);
-      } else {
-        C_active_additional[i, j] = 0;
-      }
-      // Next generate the passive cases
-      gq_tmp = C[i, j] - I_passive[i, j] - I_active[i, j]
-        - C_active_additional[i, j];
-      if (gq_tmp > 0) {
-        C_passive_additional[i, j] = poisson_rng(gq_tmp);
-      } else {
-        C_passive_additional[i, j] = 0;
-      }
-    }
-  }
-}
diff --git a/man/estimate_severity.Rd b/man/estimate_severity.Rd
deleted file mode 100644
index f603aa9..0000000
--- a/man/estimate_severity.Rd
+++ /dev/null
@@ -1,97 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/estimate_severity.R
-\name{estimate_severity}
-\alias{estimate_severity}
-\title{Estimate Severity From A Line List}
-\usage{
-estimate_severity(
-  linelist,
-  population,
-  surveillance,
-  outcome,
-  time_period = character(),
-  strata = character(),
-  population_value = utils::tail(names(population), n = 1L),
-  surveillance_reference = NULL,
-  outcome_reference = NULL,
-  time_period_reference = NULL,
-  strata_reference = NULL,
-  hazard_std = 3,
-  degrees_of_freedom = 1L,
-  active_prior = c(alpha = 1, beta = 25),
-  passive_asymptomatic_prior = c(alpha = 1, beta = 2),
-  passive_symptomatic_prior = c(alpha = 25, beta = 1),
-  ...
-)
-}
-\arguments{
-\item{linelist}{A \code{data.frame} (or \code{data.frame} extending object like a
-\code{tibble}) of the line list data.}
-
-\item{population}{A \code{data.frame} (or \code{data.frame} extending object like a
-\code{tibble}) of the population data.}
-
-\item{surveillance}{A character of columns describing the surveillance
-methods. The values in this column must be coercible to 'Active' and
-'Passive'.}
-
-\item{outcome}{A character of columns describing the outcome of the line list
-entry. The values in this column must be coercible to 'Asymptomatic',
-'Death', and 'Symptomatic'.}
-
-\item{time_period}{A character of columns describing the time period such as
-'week' or 'day'}
-
-\item{strata}{A character of columns describing the attributes to stratify
-the data on.}
-
-\item{population_value}{A unit length character corresponding to the column
-in the \code{population} data.frame (or data.frame like object) for the population
-value.}
-
-\item{surveillance_reference}{Either \code{NULL}, a \code{data.frame}, or a vector of
-values to use as a reference for the surveillance dimension. If a vector then
-it is expected that \code{surveillance} is unit length.}
-
-\item{outcome_reference}{Either \code{NULL}, a \code{data.frame}, or a vector of
-values to use as a reference for the outcome dimension. If a vector then
-it is expected that \code{outcome} is unit length.}
-
-\item{time_period_reference}{Either \code{NULL}, a \code{data.frame}, or a vector of
-values to use as a reference for the time period dimension. If a vector then
-it is expected that \code{time_period} is unit length.}
-
-\item{strata_reference}{Either \code{NULL}, a \code{data.frame}, or a vector of
-values to use as a reference for the strata dimension. If a vector then
-it is expected that \code{strata} is unit length.}
-
-\item{hazard_std}{A single length numeric greater than zero for the standard
-deviation to use in the community hazard prior distribution.}
-
-\item{degrees_of_freedom}{A single length integer greater than zero for the
-degrees of freedom to use in the model.}
-
-\item{active_prior}{The parameters for the prior beta distribution for the
-active detection probability. Can be specified as 'alpha'/'beta',
-'mean'/'var', 'mean'/'sd', or 'mean'/'concentration'.}
-
-\item{passive_asymptomatic_prior}{The parameters for the prior beta
-distribution for the passive asymptomatic detection probability. Can be
-specified as 'alpha'/'beta', 'mean'/'var', 'mean'/'sd', or
-'mean'/'concentration'.}
-
-\item{passive_symptomatic_prior}{The parameters for the prior beta
-distribution for the passive symptomatic detection probability. Can be
-specified as 'alpha'/'beta', 'mean'/'var', 'mean'/'sd', or
-'mean'/'concentration'.}
-
-\item{...}{Further optional args that are eventually given to
-\code{\link[rstan:sampling]{rstan::sampling()}} related to fitting.}
-}
-\value{
-A \linkS4class{SeverityEstimateFit} S4 object.
-}
-\description{
-Fit a severity estimate model accounting for under reporting of asymptomatic
-and mildly symptomatic cases.
-}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 974e4ef..4d4fe3c 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -13,11 +13,9 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 
 
 RcppExport SEXP _rcpp_module_boot_stan_fit4estimate_severity_mod();
-RcppExport SEXP _rcpp_module_boot_stan_fit4severity_estimate_by_surveillance_mod();
 
 static const R_CallMethodDef CallEntries[] = {
     {"_rcpp_module_boot_stan_fit4estimate_severity_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4estimate_severity_mod, 0},
-    {"_rcpp_module_boot_stan_fit4severity_estimate_by_surveillance_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4severity_estimate_by_surveillance_mod, 0},
     {NULL, NULL, 0}
 };
 
diff --git a/src/stanExports_severity_estimate_by_surveillance.cc b/src/stanExports_severity_estimate_by_surveillance.cc
deleted file mode 100644
index bee239d..0000000
--- a/src/stanExports_severity_estimate_by_surveillance.cc
+++ /dev/null
@@ -1,32 +0,0 @@
-// Generated by rstantools.  Do not edit by hand.
-
-#include <Rcpp.h>
-using namespace Rcpp ;
-#include "stanExports_severity_estimate_by_surveillance.h"
-
-RCPP_MODULE(stan_fit4severity_estimate_by_surveillance_mod) {
-
-
-    class_<rstan::stan_fit<stan_model, boost::ecuyer1988> >("rstantools_model_severity_estimate_by_surveillance")
-
-    .constructor<SEXP,SEXP,SEXP>()
-
-
-    .method("call_sampler", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::call_sampler)
-    .method("param_names", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_names)
-    .method("param_names_oi", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_names_oi)
-    .method("param_fnames_oi", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_fnames_oi)
-    .method("param_dims", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_dims)
-    .method("param_dims_oi", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_dims_oi)
-    .method("update_param_oi", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::update_param_oi)
-    .method("param_oi_tidx", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::param_oi_tidx)
-    .method("grad_log_prob", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::grad_log_prob)
-    .method("log_prob", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::log_prob)
-    .method("unconstrain_pars", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::unconstrain_pars)
-    .method("constrain_pars", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::constrain_pars)
-    .method("num_pars_unconstrained", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::num_pars_unconstrained)
-    .method("unconstrained_param_names", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::unconstrained_param_names)
-    .method("constrained_param_names", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::constrained_param_names)
-    .method("standalone_gqs", &rstan::stan_fit<stan_model, boost::ecuyer1988> ::standalone_gqs)
-    ;
-}
diff --git a/src/stanExports_severity_estimate_by_surveillance.h b/src/stanExports_severity_estimate_by_surveillance.h
deleted file mode 100644
index dfadc91..0000000
--- a/src/stanExports_severity_estimate_by_surveillance.h
+++ /dev/null
@@ -1,1572 +0,0 @@
-// Generated by rstantools.  Do not edit by hand.
-
-#ifndef MODELS_HPP
-#define MODELS_HPP
-#define STAN__SERVICES__COMMAND_HPP
-#ifndef USE_STANC3
-#define USE_STANC3
-#endif
-#include <rstan/rstaninc.hpp>
-// Code generated by stanc v2.32.2
-#include <stan/model/model_header.hpp>
-namespace model_severity_estimate_by_surveillance_namespace {
-using stan::model::model_base_crtp;
-using namespace stan::math;
-stan::math::profile_map profiles__;
-static constexpr std::array<const char*, 126> locations_array__ =
-  {" (found before start of program)",
-  " (in 'severity_estimate_by_surveillance', line 47, column 2 to column 43)",
-  " (in 'severity_estimate_by_surveillance', line 48, column 2 to column 47)",
-  " (in 'severity_estimate_by_surveillance', line 50, column 2 to column 51)",
-  " (in 'severity_estimate_by_surveillance', line 52, column 2 to column 43)",
-  " (in 'severity_estimate_by_surveillance', line 54, column 2 to column 57)",
-  " (in 'severity_estimate_by_surveillance', line 55, column 2 to column 56)",
-  " (in 'severity_estimate_by_surveillance', line 60, column 2 to column 49)",
-  " (in 'severity_estimate_by_surveillance', line 61, column 2 to column 56)",
-  " (in 'severity_estimate_by_surveillance', line 63, column 2 to column 52)",
-  " (in 'severity_estimate_by_surveillance', line 64, column 2 to column 52)",
-  " (in 'severity_estimate_by_surveillance', line 66, column 2 to column 43)",
-  " (in 'severity_estimate_by_surveillance', line 67, column 2 to column 51)",
-  " (in 'severity_estimate_by_surveillance', line 68, column 2 to column 14)",
-  " (in 'severity_estimate_by_surveillance', line 69, column 2 to column 16)",
-  " (in 'severity_estimate_by_surveillance', line 147, column 2 to column 69)",
-  " (in 'severity_estimate_by_surveillance', line 148, column 2 to column 70)",
-  " (in 'severity_estimate_by_surveillance', line 149, column 2 to column 14)",
-  " (in 'severity_estimate_by_surveillance', line 73, column 4 to column 22)",
-  " (in 'severity_estimate_by_surveillance', line 74, column 4 to column 28)",
-  " (in 'severity_estimate_by_surveillance', line 76, column 6 to column 57)",
-  " (in 'severity_estimate_by_surveillance', line 77, column 6 to column 65)",
-  " (in 'severity_estimate_by_surveillance', line 75, column 36 to line 78, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 75, column 4 to line 78, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 79, column 4 to column 30)",
-  " (in 'severity_estimate_by_surveillance', line 80, column 4 to column 39)",
-  " (in 'severity_estimate_by_surveillance', line 82, column 4 to column 28)",
-  " (in 'severity_estimate_by_surveillance', line 83, column 4 to column 57)",
-  " (in 'severity_estimate_by_surveillance', line 71, column 29 to line 84, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 71, column 2 to line 84, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 89, column 6 to column 42)",
-  " (in 'severity_estimate_by_surveillance', line 90, column 6 to column 53)",
-  " (in 'severity_estimate_by_surveillance', line 88, column 31 to line 91, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 88, column 4 to line 91, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 87, column 27 to line 92, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 87, column 2 to line 92, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 95, column 4 to line 96, column 48)",
-  " (in 'severity_estimate_by_surveillance', line 97, column 4 to column 71)",
-  " (in 'severity_estimate_by_surveillance', line 94, column 29 to line 98, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 94, column 2 to line 98, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 154, column 6 to column 40)",
-  " (in 'severity_estimate_by_surveillance', line 158, column 8 to column 38)",
-  " (in 'severity_estimate_by_surveillance', line 157, column 13 to line 159, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 156, column 8 to column 75)",
-  " (in 'severity_estimate_by_surveillance', line 155, column 22 to line 157, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 155, column 6 to line 159, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 161, column 6 to line 162, column 36)",
-  " (in 'severity_estimate_by_surveillance', line 166, column 8 to column 39)",
-  " (in 'severity_estimate_by_surveillance', line 165, column 13 to line 167, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 164, column 8 to column 57)",
-  " (in 'severity_estimate_by_surveillance', line 163, column 22 to line 165, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 163, column 6 to line 167, column 7)",
-  " (in 'severity_estimate_by_surveillance', line 152, column 31 to line 168, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 152, column 4 to line 168, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 151, column 27 to line 169, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 151, column 2 to line 169, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 104, column 4 to column 36)",
-  " (in 'severity_estimate_by_surveillance', line 105, column 4 to column 40)",
-  " (in 'severity_estimate_by_surveillance', line 103, column 40 to line 106, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 103, column 2 to line 106, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 108, column 2 to column 73)",
-  " (in 'severity_estimate_by_surveillance', line 109, column 2 to line 111, column 4)",
-  " (in 'severity_estimate_by_surveillance', line 112, column 2 to line 114, column 4)",
-  " (in 'severity_estimate_by_surveillance', line 118, column 6 to line 122, column 8)",
-  " (in 'severity_estimate_by_surveillance', line 123, column 6 to column 59)",
-  " (in 'severity_estimate_by_surveillance', line 124, column 6 to column 79)",
-  " (in 'severity_estimate_by_surveillance', line 117, column 31 to line 125, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 117, column 4 to line 125, column 5)",
-  " (in 'severity_estimate_by_surveillance', line 116, column 27 to line 126, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 116, column 2 to line 126, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 129, column 4 to column 57)",
-  " (in 'severity_estimate_by_surveillance', line 130, column 4 to column 60)",
-  " (in 'severity_estimate_by_surveillance', line 128, column 31 to line 131, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 128, column 2 to line 131, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 134, column 4 to line 138, column 4)",
-  " (in 'severity_estimate_by_surveillance', line 139, column 2 to line 141, column 4)",
-  " (in 'severity_estimate_by_surveillance', line 133, column 32 to line 142, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 133, column 2 to line 142, column 3)",
-  " (in 'severity_estimate_by_surveillance', line 4, column 2 to column 30)",
-  " (in 'severity_estimate_by_surveillance', line 6, column 2 to column 28)",
-  " (in 'severity_estimate_by_surveillance', line 8, column 2 to column 32)",
-  " (in 'severity_estimate_by_surveillance', line 10, column 2 to column 33)",
-  " (in 'severity_estimate_by_surveillance', line 13, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 13, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 13, column 2 to column 58)",
-  " (in 'severity_estimate_by_surveillance', line 15, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 15, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 15, column 2 to column 59)",
-  " (in 'severity_estimate_by_surveillance', line 17, column 8 to column 21)",
-  " (in 'severity_estimate_by_surveillance', line 17, column 2 to column 47)",
-  " (in 'severity_estimate_by_surveillance', line 20, column 8 to column 23)",
-  " (in 'severity_estimate_by_surveillance', line 20, column 2 to column 52)",
-  " (in 'severity_estimate_by_surveillance', line 22, column 8 to column 23)",
-  " (in 'severity_estimate_by_surveillance', line 22, column 2 to column 54)",
-  " (in 'severity_estimate_by_surveillance', line 24, column 8 to column 23)",
-  " (in 'severity_estimate_by_surveillance', line 24, column 2 to column 50)",
-  " (in 'severity_estimate_by_surveillance', line 26, column 8 to column 24)",
-  " (in 'severity_estimate_by_surveillance', line 26, column 2 to column 54)",
-  " (in 'severity_estimate_by_surveillance', line 28, column 8 to column 24)",
-  " (in 'severity_estimate_by_surveillance', line 28, column 2 to column 56)",
-  " (in 'severity_estimate_by_surveillance', line 30, column 8 to column 24)",
-  " (in 'severity_estimate_by_surveillance', line 30, column 2 to column 52)",
-  " (in 'severity_estimate_by_surveillance', line 33, column 2 to column 28)",
-  " (in 'severity_estimate_by_surveillance', line 35, column 2 to column 35)",
-  " (in 'severity_estimate_by_surveillance', line 37, column 2 to column 40)",
-  " (in 'severity_estimate_by_surveillance', line 38, column 2 to column 39)",
-  " (in 'severity_estimate_by_surveillance', line 40, column 2 to column 44)",
-  " (in 'severity_estimate_by_surveillance', line 41, column 2 to column 43)",
-  " (in 'severity_estimate_by_surveillance', line 42, column 2 to column 43)",
-  " (in 'severity_estimate_by_surveillance', line 43, column 2 to column 42)",
-  " (in 'severity_estimate_by_surveillance', line 47, column 8 to column 30)",
-  " (in 'severity_estimate_by_surveillance', line 48, column 8 to column 30)",
-  " (in 'severity_estimate_by_surveillance', line 50, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 50, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 60, column 8 to column 21)",
-  " (in 'severity_estimate_by_surveillance', line 61, column 8 to column 21)",
-  " (in 'severity_estimate_by_surveillance', line 63, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 63, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 64, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 64, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 66, column 8 to column 21)",
-  " (in 'severity_estimate_by_surveillance', line 67, column 8 to column 21)",
-  " (in 'severity_estimate_by_surveillance', line 147, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 147, column 21 to column 34)",
-  " (in 'severity_estimate_by_surveillance', line 148, column 8 to column 19)",
-  " (in 'severity_estimate_by_surveillance', line 148, column 21 to column 34)"};
-#include <stan_meta_header.hpp>
-class model_severity_estimate_by_surveillance final : public model_base_crtp<model_severity_estimate_by_surveillance> {
-private:
-  int strata_groups;
-  int time_groups;
-  int observed_active;
-  int observed_passive;
-  std::vector<std::vector<int>> I_active;
-  std::vector<std::vector<int>> I_passive;
-  std::vector<int> population;
-  std::vector<int> strata_active;
-  std::vector<int> symptoms_active;
-  std::vector<int> dead_active;
-  std::vector<int> strata_passive;
-  std::vector<int> symptoms_passive;
-  std::vector<int> dead_passive;
-  double hazard_std;
-  int degrees_of_freedom;
-  double active_detection_alpha;
-  double active_detection_beta;
-  double passive_asymptomatic_alpha;
-  double passive_asymptomatic_beta;
-  double passive_symptomatic_alpha;
-  double passive_symptomatic_beta;
-  int alpha_1dim__;
-  int mort_coef_1dim__;
-public:
-  ~model_severity_estimate_by_surveillance() {}
-  model_severity_estimate_by_surveillance(stan::io::var_context& context__,
-                                          unsigned int random_seed__ = 0,
-                                          std::ostream* pstream__ = nullptr)
-      : model_base_crtp(0) {
-    int current_statement__ = 0;
-    using local_scalar_t__ = double;
-    boost::ecuyer1988 base_rng__ =
-      stan::services::util::create_rng(random_seed__, 0);
-    // suppress unused var warning
-    (void) base_rng__;
-    static constexpr const char* function__ =
-      "model_severity_estimate_by_surveillance_namespace::model_severity_estimate_by_surveillance";
-    // suppress unused var warning
-    (void) function__;
-    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
-    // suppress unused var warning
-    (void) DUMMY_VAR__;
-    try {
-      int pos__ = std::numeric_limits<int>::min();
-      pos__ = 1;
-      current_statement__ = 78;
-      context__.validate_dims("data initialization", "strata_groups", "int",
-        std::vector<size_t>{});
-      strata_groups = std::numeric_limits<int>::min();
-      current_statement__ = 78;
-      strata_groups = context__.vals_i("strata_groups")[(1 - 1)];
-      current_statement__ = 78;
-      stan::math::check_greater_or_equal(function__, "strata_groups",
-        strata_groups, 0);
-      current_statement__ = 79;
-      context__.validate_dims("data initialization", "time_groups", "int",
-        std::vector<size_t>{});
-      time_groups = std::numeric_limits<int>::min();
-      current_statement__ = 79;
-      time_groups = context__.vals_i("time_groups")[(1 - 1)];
-      current_statement__ = 79;
-      stan::math::check_greater_or_equal(function__, "time_groups",
-        time_groups, 0);
-      current_statement__ = 80;
-      context__.validate_dims("data initialization", "observed_active",
-        "int", std::vector<size_t>{});
-      observed_active = std::numeric_limits<int>::min();
-      current_statement__ = 80;
-      observed_active = context__.vals_i("observed_active")[(1 - 1)];
-      current_statement__ = 80;
-      stan::math::check_greater_or_equal(function__, "observed_active",
-        observed_active, 0);
-      current_statement__ = 81;
-      context__.validate_dims("data initialization", "observed_passive",
-        "int", std::vector<size_t>{});
-      observed_passive = std::numeric_limits<int>::min();
-      current_statement__ = 81;
-      observed_passive = context__.vals_i("observed_passive")[(1 - 1)];
-      current_statement__ = 81;
-      stan::math::check_greater_or_equal(function__, "observed_passive",
-        observed_passive, 0);
-      current_statement__ = 82;
-      stan::math::validate_non_negative_index("I_active", "time_groups",
-        time_groups);
-      current_statement__ = 83;
-      stan::math::validate_non_negative_index("I_active", "strata_groups",
-        strata_groups);
-      current_statement__ = 84;
-      context__.validate_dims("data initialization", "I_active", "int",
-        std::vector<size_t>{static_cast<size_t>(time_groups),
-          static_cast<size_t>(strata_groups)});
-      I_active = std::vector<std::vector<int>>(time_groups,
-                   std::vector<int>(strata_groups,
-                     std::numeric_limits<int>::min()));
-      {
-        std::vector<int> I_active_flat__;
-        current_statement__ = 84;
-        I_active_flat__ = context__.vals_i("I_active");
-        current_statement__ = 84;
-        pos__ = 1;
-        current_statement__ = 84;
-        for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-          current_statement__ = 84;
-          for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-            current_statement__ = 84;
-            stan::model::assign(I_active, I_active_flat__[(pos__ - 1)],
-              "assigning variable I_active", stan::model::index_uni(sym2__),
-              stan::model::index_uni(sym1__));
-            current_statement__ = 84;
-            pos__ = (pos__ + 1);
-          }
-        }
-      }
-      current_statement__ = 84;
-      stan::math::check_greater_or_equal(function__, "I_active", I_active, 0);
-      current_statement__ = 85;
-      stan::math::validate_non_negative_index("I_passive", "time_groups",
-        time_groups);
-      current_statement__ = 86;
-      stan::math::validate_non_negative_index("I_passive", "strata_groups",
-        strata_groups);
-      current_statement__ = 87;
-      context__.validate_dims("data initialization", "I_passive", "int",
-        std::vector<size_t>{static_cast<size_t>(time_groups),
-          static_cast<size_t>(strata_groups)});
-      I_passive = std::vector<std::vector<int>>(time_groups,
-                    std::vector<int>(strata_groups,
-                      std::numeric_limits<int>::min()));
-      {
-        std::vector<int> I_passive_flat__;
-        current_statement__ = 87;
-        I_passive_flat__ = context__.vals_i("I_passive");
-        current_statement__ = 87;
-        pos__ = 1;
-        current_statement__ = 87;
-        for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-          current_statement__ = 87;
-          for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-            current_statement__ = 87;
-            stan::model::assign(I_passive, I_passive_flat__[(pos__ - 1)],
-              "assigning variable I_passive", stan::model::index_uni(sym2__),
-              stan::model::index_uni(sym1__));
-            current_statement__ = 87;
-            pos__ = (pos__ + 1);
-          }
-        }
-      }
-      current_statement__ = 87;
-      stan::math::check_greater_or_equal(function__, "I_passive", I_passive,
-        0);
-      current_statement__ = 88;
-      stan::math::validate_non_negative_index("population", "strata_groups",
-        strata_groups);
-      current_statement__ = 89;
-      context__.validate_dims("data initialization", "population", "int",
-        std::vector<size_t>{static_cast<size_t>(strata_groups)});
-      population = std::vector<int>(strata_groups,
-                     std::numeric_limits<int>::min());
-      current_statement__ = 89;
-      population = context__.vals_i("population");
-      current_statement__ = 89;
-      stan::math::check_greater_or_equal(function__, "population",
-        population, 0);
-      current_statement__ = 90;
-      stan::math::validate_non_negative_index("strata_active",
-        "observed_active", observed_active);
-      current_statement__ = 91;
-      context__.validate_dims("data initialization", "strata_active", "int",
-        std::vector<size_t>{static_cast<size_t>(observed_active)});
-      strata_active = std::vector<int>(observed_active,
-                        std::numeric_limits<int>::min());
-      current_statement__ = 91;
-      strata_active = context__.vals_i("strata_active");
-      current_statement__ = 91;
-      stan::math::check_greater_or_equal(function__, "strata_active",
-        strata_active, 0);
-      current_statement__ = 92;
-      stan::math::validate_non_negative_index("symptoms_active",
-        "observed_active", observed_active);
-      current_statement__ = 93;
-      context__.validate_dims("data initialization", "symptoms_active",
-        "int", std::vector<size_t>{static_cast<size_t>(observed_active)});
-      symptoms_active = std::vector<int>(observed_active,
-                          std::numeric_limits<int>::min());
-      current_statement__ = 93;
-      symptoms_active = context__.vals_i("symptoms_active");
-      current_statement__ = 93;
-      stan::math::check_greater_or_equal(function__, "symptoms_active",
-        symptoms_active, 0);
-      current_statement__ = 94;
-      stan::math::validate_non_negative_index("dead_active",
-        "observed_active", observed_active);
-      current_statement__ = 95;
-      context__.validate_dims("data initialization", "dead_active", "int",
-        std::vector<size_t>{static_cast<size_t>(observed_active)});
-      dead_active = std::vector<int>(observed_active,
-                      std::numeric_limits<int>::min());
-      current_statement__ = 95;
-      dead_active = context__.vals_i("dead_active");
-      current_statement__ = 95;
-      stan::math::check_greater_or_equal(function__, "dead_active",
-        dead_active, 0);
-      current_statement__ = 96;
-      stan::math::validate_non_negative_index("strata_passive",
-        "observed_passive", observed_passive);
-      current_statement__ = 97;
-      context__.validate_dims("data initialization", "strata_passive", "int",
-        std::vector<size_t>{static_cast<size_t>(observed_passive)});
-      strata_passive = std::vector<int>(observed_passive,
-                         std::numeric_limits<int>::min());
-      current_statement__ = 97;
-      strata_passive = context__.vals_i("strata_passive");
-      current_statement__ = 97;
-      stan::math::check_greater_or_equal(function__, "strata_passive",
-        strata_passive, 0);
-      current_statement__ = 98;
-      stan::math::validate_non_negative_index("symptoms_passive",
-        "observed_passive", observed_passive);
-      current_statement__ = 99;
-      context__.validate_dims("data initialization", "symptoms_passive",
-        "int", std::vector<size_t>{static_cast<size_t>(observed_passive)});
-      symptoms_passive = std::vector<int>(observed_passive,
-                           std::numeric_limits<int>::min());
-      current_statement__ = 99;
-      symptoms_passive = context__.vals_i("symptoms_passive");
-      current_statement__ = 99;
-      stan::math::check_greater_or_equal(function__, "symptoms_passive",
-        symptoms_passive, 0);
-      current_statement__ = 100;
-      stan::math::validate_non_negative_index("dead_passive",
-        "observed_passive", observed_passive);
-      current_statement__ = 101;
-      context__.validate_dims("data initialization", "dead_passive", "int",
-        std::vector<size_t>{static_cast<size_t>(observed_passive)});
-      dead_passive = std::vector<int>(observed_passive,
-                       std::numeric_limits<int>::min());
-      current_statement__ = 101;
-      dead_passive = context__.vals_i("dead_passive");
-      current_statement__ = 101;
-      stan::math::check_greater_or_equal(function__, "dead_passive",
-        dead_passive, 0);
-      current_statement__ = 102;
-      context__.validate_dims("data initialization", "hazard_std", "double",
-        std::vector<size_t>{});
-      hazard_std = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 102;
-      hazard_std = context__.vals_r("hazard_std")[(1 - 1)];
-      current_statement__ = 102;
-      stan::math::check_greater_or_equal(function__, "hazard_std",
-        hazard_std, 0);
-      current_statement__ = 103;
-      context__.validate_dims("data initialization", "degrees_of_freedom",
-        "int", std::vector<size_t>{});
-      degrees_of_freedom = std::numeric_limits<int>::min();
-      current_statement__ = 103;
-      degrees_of_freedom = context__.vals_i("degrees_of_freedom")[(1 - 1)];
-      current_statement__ = 103;
-      stan::math::check_greater_or_equal(function__, "degrees_of_freedom",
-        degrees_of_freedom, 1);
-      current_statement__ = 104;
-      context__.validate_dims("data initialization",
-        "active_detection_alpha", "double", std::vector<size_t>{});
-      active_detection_alpha = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 104;
-      active_detection_alpha = context__.vals_r("active_detection_alpha")[(1
-        - 1)];
-      current_statement__ = 104;
-      stan::math::check_greater_or_equal(function__,
-        "active_detection_alpha", active_detection_alpha, 0);
-      current_statement__ = 105;
-      context__.validate_dims("data initialization", "active_detection_beta",
-        "double", std::vector<size_t>{});
-      active_detection_beta = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 105;
-      active_detection_beta = context__.vals_r("active_detection_beta")[(1 -
-        1)];
-      current_statement__ = 105;
-      stan::math::check_greater_or_equal(function__, "active_detection_beta",
-        active_detection_beta, 0);
-      current_statement__ = 106;
-      context__.validate_dims("data initialization",
-        "passive_asymptomatic_alpha", "double", std::vector<size_t>{});
-      passive_asymptomatic_alpha = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 106;
-      passive_asymptomatic_alpha = context__.vals_r("passive_asymptomatic_alpha")[(1
-        - 1)];
-      current_statement__ = 106;
-      stan::math::check_greater_or_equal(function__,
-        "passive_asymptomatic_alpha", passive_asymptomatic_alpha, 0);
-      current_statement__ = 107;
-      context__.validate_dims("data initialization",
-        "passive_asymptomatic_beta", "double", std::vector<size_t>{});
-      passive_asymptomatic_beta = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 107;
-      passive_asymptomatic_beta = context__.vals_r("passive_asymptomatic_beta")[(1
-        - 1)];
-      current_statement__ = 107;
-      stan::math::check_greater_or_equal(function__,
-        "passive_asymptomatic_beta", passive_asymptomatic_beta, 0);
-      current_statement__ = 108;
-      context__.validate_dims("data initialization",
-        "passive_symptomatic_alpha", "double", std::vector<size_t>{});
-      passive_symptomatic_alpha = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 108;
-      passive_symptomatic_alpha = context__.vals_r("passive_symptomatic_alpha")[(1
-        - 1)];
-      current_statement__ = 108;
-      stan::math::check_greater_or_equal(function__,
-        "passive_symptomatic_alpha", passive_symptomatic_alpha, 0);
-      current_statement__ = 109;
-      context__.validate_dims("data initialization",
-        "passive_symptomatic_beta", "double", std::vector<size_t>{});
-      passive_symptomatic_beta = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 109;
-      passive_symptomatic_beta = context__.vals_r("passive_symptomatic_beta")[(1
-        - 1)];
-      current_statement__ = 109;
-      stan::math::check_greater_or_equal(function__,
-        "passive_symptomatic_beta", passive_symptomatic_beta, 0);
-      current_statement__ = 110;
-      alpha_1dim__ = std::numeric_limits<int>::min();
-      current_statement__ = 110;
-      alpha_1dim__ = (1 + degrees_of_freedom);
-      current_statement__ = 110;
-      stan::math::validate_non_negative_index("alpha",
-        "1 + degrees_of_freedom", alpha_1dim__);
-      current_statement__ = 111;
-      mort_coef_1dim__ = std::numeric_limits<int>::min();
-      current_statement__ = 111;
-      mort_coef_1dim__ = (1 + degrees_of_freedom);
-      current_statement__ = 111;
-      stan::math::validate_non_negative_index("mort_coef",
-        "1 + degrees_of_freedom", mort_coef_1dim__);
-      current_statement__ = 112;
-      stan::math::validate_non_negative_index("logit_hzd", "time_groups",
-        time_groups);
-      current_statement__ = 113;
-      stan::math::validate_non_negative_index("logit_hzd", "strata_groups",
-        strata_groups);
-      current_statement__ = 114;
-      stan::math::validate_non_negative_index("xi", "strata_groups",
-        strata_groups);
-      current_statement__ = 115;
-      stan::math::validate_non_negative_index("mortality", "strata_groups",
-        strata_groups);
-      current_statement__ = 116;
-      stan::math::validate_non_negative_index("S", "time_groups", time_groups);
-      current_statement__ = 117;
-      stan::math::validate_non_negative_index("S", "strata_groups",
-        strata_groups);
-      current_statement__ = 118;
-      stan::math::validate_non_negative_index("C", "time_groups", time_groups);
-      current_statement__ = 119;
-      stan::math::validate_non_negative_index("C", "strata_groups",
-        strata_groups);
-      current_statement__ = 120;
-      stan::math::validate_non_negative_index("theta", "strata_groups",
-        strata_groups);
-      current_statement__ = 121;
-      stan::math::validate_non_negative_index("passive_denom",
-        "strata_groups", strata_groups);
-      current_statement__ = 122;
-      stan::math::validate_non_negative_index("C_active_additional",
-        "time_groups", time_groups);
-      current_statement__ = 123;
-      stan::math::validate_non_negative_index("C_active_additional",
-        "strata_groups", strata_groups);
-      current_statement__ = 124;
-      stan::math::validate_non_negative_index("C_passive_additional",
-        "time_groups", time_groups);
-      current_statement__ = 125;
-      stan::math::validate_non_negative_index("C_passive_additional",
-        "strata_groups", strata_groups);
-    } catch (const std::exception& e) {
-      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
-    }
-    num_params_r__ = alpha_1dim__ + mort_coef_1dim__ + (time_groups *
-      strata_groups) + 1 + 1 + 1;
-  }
-  inline std::string model_name() const final {
-    return "model_severity_estimate_by_surveillance";
-  }
-  inline std::vector<std::string> model_compile_info() const noexcept {
-    return std::vector<std::string>{"stanc_version = stanc3 v2.32.2",
-             "stancflags = --allow-undefined"};
-  }
-  template <bool propto__, bool jacobian__, typename VecR, typename VecI,
-            stan::require_vector_like_t<VecR>* = nullptr,
-            stan::require_vector_like_vt<std::is_integral, VecI>* = nullptr>
-  inline stan::scalar_type_t<VecR>
-  log_prob_impl(VecR& params_r__, VecI& params_i__, std::ostream*
-                pstream__ = nullptr) const {
-    using T__ = stan::scalar_type_t<VecR>;
-    using local_scalar_t__ = T__;
-    T__ lp__(0.0);
-    stan::math::accumulator<T__> lp_accum__;
-    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
-    int current_statement__ = 0;
-    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
-    // suppress unused var warning
-    (void) DUMMY_VAR__;
-    static constexpr const char* function__ =
-      "model_severity_estimate_by_surveillance_namespace::log_prob";
-    // suppress unused var warning
-    (void) function__;
-    try {
-      std::vector<local_scalar_t__> alpha =
-        std::vector<local_scalar_t__>(alpha_1dim__, DUMMY_VAR__);
-      current_statement__ = 1;
-      alpha = in__.template read<std::vector<local_scalar_t__>>(alpha_1dim__);
-      std::vector<local_scalar_t__> mort_coef =
-        std::vector<local_scalar_t__>(mort_coef_1dim__, DUMMY_VAR__);
-      current_statement__ = 2;
-      mort_coef = in__.template read<
-                    std::vector<local_scalar_t__>>(mort_coef_1dim__);
-      std::vector<std::vector<local_scalar_t__>> logit_hzd =
-        std::vector<std::vector<local_scalar_t__>>(time_groups,
-          std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__));
-      current_statement__ = 3;
-      logit_hzd = in__.template read<
-                    std::vector<std::vector<local_scalar_t__>>>(time_groups,
-                    strata_groups);
-      local_scalar_t__ active_detection = DUMMY_VAR__;
-      current_statement__ = 4;
-      active_detection = in__.template read_constrain_lub<local_scalar_t__,
-                           jacobian__>(0, 1, lp__);
-      local_scalar_t__ passive_asymptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 5;
-      passive_asymptomatic_detection = in__.template read_constrain_lub<
-                                         local_scalar_t__, jacobian__>(0, 1,
-                                         lp__);
-      local_scalar_t__ passive_symptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 6;
-      passive_symptomatic_detection = in__.template read_constrain_lub<
-                                        local_scalar_t__, jacobian__>(0, 1,
-                                        lp__);
-      std::vector<local_scalar_t__> xi =
-        std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__);
-      std::vector<local_scalar_t__> mortality =
-        std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__);
-      std::vector<std::vector<local_scalar_t__>> S =
-        std::vector<std::vector<local_scalar_t__>>(time_groups,
-          std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__));
-      std::vector<std::vector<local_scalar_t__>> C =
-        std::vector<std::vector<local_scalar_t__>>(time_groups,
-          std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__));
-      std::vector<local_scalar_t__> theta =
-        std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__);
-      std::vector<local_scalar_t__> passive_denom =
-        std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__);
-      local_scalar_t__ xi_tmp = DUMMY_VAR__;
-      local_scalar_t__ mort_tmp = DUMMY_VAR__;
-      current_statement__ = 29;
-      for (int i = 1; i <= strata_groups; ++i) {
-        current_statement__ = 18;
-        xi_tmp = stan::model::rvalue(alpha, "alpha",
-                   stan::model::index_uni(1));
-        current_statement__ = 19;
-        mort_tmp = stan::model::rvalue(mort_coef, "mort_coef",
-                     stan::model::index_uni(1));
-        current_statement__ = 23;
-        for (int j = 1; j <= degrees_of_freedom; ++j) {
-          current_statement__ = 20;
-          xi_tmp = (xi_tmp +
-            (stan::model::rvalue(alpha, "alpha",
-               stan::model::index_uni((j + 1)))
-            * stan::math::pow((i - 1.0), j)));
-          current_statement__ = 21;
-          mort_tmp = (mort_tmp +
-            (stan::model::rvalue(mort_coef, "mort_coef",
-               stan::model::index_uni((j + 1)))
-            * stan::math::pow((i - 1.0), j)));
-        }
-        current_statement__ = 24;
-        stan::model::assign(xi, stan::math::inv_logit(xi_tmp),
-          "assigning variable xi", stan::model::index_uni(i));
-        current_statement__ = 25;
-        stan::model::assign(mortality, stan::math::inv_logit(mort_tmp),
-          "assigning variable mortality", stan::model::index_uni(i));
-        current_statement__ = 26;
-        stan::model::assign(S,
-          stan::model::rvalue(population, "population",
-            stan::model::index_uni(i)), "assigning variable S",
-          stan::model::index_uni(1), stan::model::index_uni(i));
-        current_statement__ = 27;
-        stan::model::assign(C,
-          (stan::model::rvalue(population, "population",
-             stan::model::index_uni(i)) *
-          stan::math::inv_logit(
-            stan::model::rvalue(logit_hzd, "logit_hzd",
-              stan::model::index_uni(1), stan::model::index_uni(i)))),
-          "assigning variable C", stan::model::index_uni(1),
-          stan::model::index_uni(i));
-      }
-      current_statement__ = 35;
-      for (int i = 2; i <= time_groups; ++i) {
-        current_statement__ = 33;
-        for (int j = 1; j <= strata_groups; ++j) {
-          current_statement__ = 30;
-          stan::model::assign(S,
-            (stan::model::rvalue(S, "S", stan::model::index_uni((i - 1)),
-               stan::model::index_uni(j)) -
-            stan::model::rvalue(C, "C", stan::model::index_uni((i - 1)),
-              stan::model::index_uni(j))), "assigning variable S",
-            stan::model::index_uni(i), stan::model::index_uni(j));
-          current_statement__ = 31;
-          stan::model::assign(C,
-            (stan::model::rvalue(S, "S", stan::model::index_uni(i),
-               stan::model::index_uni(j)) *
-            stan::math::inv_logit(
-              stan::model::rvalue(logit_hzd, "logit_hzd",
-                stan::model::index_uni(i), stan::model::index_uni(j)))),
-            "assigning variable C", stan::model::index_uni(i),
-            stan::model::index_uni(j));
-        }
-      }
-      current_statement__ = 39;
-      for (int i = 1; i <= strata_groups; ++i) {
-        current_statement__ = 36;
-        stan::model::assign(theta, ((passive_asymptomatic_detection * (1.0 -
-          stan::model::rvalue(xi, "xi", stan::model::index_uni(i)))) +
-          (passive_symptomatic_detection *
-          stan::model::rvalue(xi, "xi", stan::model::index_uni(i)))),
-          "assigning variable theta", stan::model::index_uni(i));
-        current_statement__ = 37;
-        stan::model::assign(passive_denom, (1.0 - ((1.0 -
-          stan::model::rvalue(mortality, "mortality",
-            stan::model::index_uni(i))) * (1.0 -
-          stan::model::rvalue(theta, "theta", stan::model::index_uni(i))))),
-          "assigning variable passive_denom", stan::model::index_uni(i));
-      }
-      current_statement__ = 7;
-      stan::math::check_greater_or_equal(function__, "xi", xi, 0);
-      current_statement__ = 7;
-      stan::math::check_less_or_equal(function__, "xi", xi, 1);
-      current_statement__ = 8;
-      stan::math::check_greater_or_equal(function__, "mortality", mortality,
-        0);
-      current_statement__ = 8;
-      stan::math::check_less_or_equal(function__, "mortality", mortality, 1);
-      current_statement__ = 9;
-      stan::math::check_greater_or_equal(function__, "S", S, 0);
-      current_statement__ = 10;
-      stan::math::check_greater_or_equal(function__, "C", C, 0);
-      current_statement__ = 11;
-      stan::math::check_greater_or_equal(function__, "theta", theta, 0);
-      current_statement__ = 12;
-      stan::math::check_greater_or_equal(function__, "passive_denom",
-        passive_denom, 0);
-      {
-        current_statement__ = 59;
-        for (int i = 1; i <= (degrees_of_freedom + 1); ++i) {
-          current_statement__ = 56;
-          lp_accum__.add(stan::math::normal_lpdf<propto__>(
-                           stan::model::rvalue(alpha, "alpha",
-                             stan::model::index_uni(i)), 0.0, 10000.0));
-          current_statement__ = 57;
-          lp_accum__.add(stan::math::normal_lpdf<propto__>(
-                           stan::model::rvalue(mort_coef, "mort_coef",
-                             stan::model::index_uni(i)), 0.0, 10000.0));
-        }
-        current_statement__ = 60;
-        lp_accum__.add(stan::math::beta_lpdf<propto__>(active_detection,
-                         active_detection_alpha, active_detection_beta));
-        current_statement__ = 61;
-        lp_accum__.add(stan::math::beta_lpdf<propto__>(
-                         passive_asymptomatic_detection,
-                         passive_asymptomatic_alpha,
-                         passive_asymptomatic_beta));
-        current_statement__ = 62;
-        lp_accum__.add(stan::math::beta_lpdf<propto__>(
-                         passive_symptomatic_detection,
-                         passive_symptomatic_alpha, passive_symptomatic_beta));
-        current_statement__ = 69;
-        for (int i = 1; i <= time_groups; ++i) {
-          current_statement__ = 67;
-          for (int j = 1; j <= strata_groups; ++j) {
-            current_statement__ = 63;
-            lp_accum__.add(stan::math::normal_lpdf<propto__>(
-                             stan::model::rvalue(logit_hzd, "logit_hzd",
-                               stan::model::index_uni(i),
-                               stan::model::index_uni(j)),
-                             stan::math::logit((((1.0 *
-                               stan::model::rvalue(I_passive, "I_passive",
-                                 stan::model::index_uni(i),
-                                 stan::model::index_uni(j))) /
-                               stan::model::rvalue(population, "population",
-                                 stan::model::index_uni(j))) + (100.0 *
-                               stan::math::machine_precision()))), hazard_std));
-            current_statement__ = 64;
-            lp_accum__.add(stan::math::poisson_lpmf<propto__>(
-                             stan::model::rvalue(I_active, "I_active",
-                               stan::model::index_uni(i),
-                               stan::model::index_uni(j)), (active_detection
-                             *
-                             stan::model::rvalue(C, "C",
-                               stan::model::index_uni(i),
-                               stan::model::index_uni(j)))));
-            current_statement__ = 65;
-            lp_accum__.add(stan::math::poisson_lpmf<propto__>(
-                             stan::model::rvalue(I_passive, "I_passive",
-                               stan::model::index_uni(i),
-                               stan::model::index_uni(j)), (((1.0 -
-                             active_detection) *
-                             stan::model::rvalue(theta, "theta",
-                               stan::model::index_uni(j))) *
-                             stan::model::rvalue(C, "C",
-                               stan::model::index_uni(i),
-                               stan::model::index_uni(j)))));
-          }
-        }
-        current_statement__ = 73;
-        for (int i = 1; i <= observed_active; ++i) {
-          current_statement__ = 70;
-          lp_accum__.add(stan::math::bernoulli_lpmf<propto__>(
-                           stan::model::rvalue(symptoms_active,
-                             "symptoms_active", stan::model::index_uni(i)),
-                           stan::model::rvalue(xi, "xi",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_active,
-                                 "strata_active", stan::model::index_uni(i))))));
-          current_statement__ = 71;
-          lp_accum__.add(stan::math::bernoulli_lpmf<propto__>(
-                           stan::model::rvalue(dead_active, "dead_active",
-                             stan::model::index_uni(i)),
-                           stan::model::rvalue(mortality, "mortality",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_active,
-                                 "strata_active", stan::model::index_uni(i))))));
-        }
-        current_statement__ = 77;
-        for (int i = 1; i <= observed_passive; ++i) {
-          current_statement__ = 74;
-          lp_accum__.add(stan::math::bernoulli_lpmf<propto__>(
-                           stan::model::rvalue(symptoms_passive,
-                             "symptoms_passive", stan::model::index_uni(i)),
-                           ((1.0 - ((1.0 -
-                           stan::model::rvalue(mortality, "mortality",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_passive,
-                                 "strata_passive", stan::model::index_uni(i)))))
-                           * (1.0 - (passive_symptomatic_detection *
-                           stan::model::rvalue(xi, "xi",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_passive,
-                                 "strata_passive", stan::model::index_uni(i))))))))
-                           /
-                           stan::model::rvalue(passive_denom,
-                             "passive_denom",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_passive,
-                                 "strata_passive", stan::model::index_uni(i)))))));
-          current_statement__ = 75;
-          lp_accum__.add(stan::math::bernoulli_lpmf<propto__>(
-                           stan::model::rvalue(dead_passive, "dead_passive",
-                             stan::model::index_uni(i)),
-                           (stan::model::rvalue(mortality, "mortality",
-                              stan::model::index_uni(
-                                stan::model::rvalue(strata_passive,
-                                  "strata_passive", stan::model::index_uni(i))))
-                           /
-                           stan::model::rvalue(passive_denom,
-                             "passive_denom",
-                             stan::model::index_uni(
-                               stan::model::rvalue(strata_passive,
-                                 "strata_passive", stan::model::index_uni(i)))))));
-        }
-      }
-    } catch (const std::exception& e) {
-      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
-    }
-    lp_accum__.add(lp__);
-    return lp_accum__.sum();
-  }
-  template <typename RNG, typename VecR, typename VecI, typename VecVar,
-            stan::require_vector_like_vt<std::is_floating_point,
-            VecR>* = nullptr, stan::require_vector_like_vt<std::is_integral,
-            VecI>* = nullptr, stan::require_vector_vt<std::is_floating_point,
-            VecVar>* = nullptr>
-  inline void
-  write_array_impl(RNG& base_rng__, VecR& params_r__, VecI& params_i__,
-                   VecVar& vars__, const bool
-                   emit_transformed_parameters__ = true, const bool
-                   emit_generated_quantities__ = true, std::ostream*
-                   pstream__ = nullptr) const {
-    using local_scalar_t__ = double;
-    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
-    stan::io::serializer<local_scalar_t__> out__(vars__);
-    static constexpr bool propto__ = true;
-    // suppress unused var warning
-    (void) propto__;
-    double lp__ = 0.0;
-    // suppress unused var warning
-    (void) lp__;
-    int current_statement__ = 0;
-    stan::math::accumulator<double> lp_accum__;
-    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
-    // suppress unused var warning
-    (void) DUMMY_VAR__;
-    constexpr bool jacobian__ = false;
-    static constexpr const char* function__ =
-      "model_severity_estimate_by_surveillance_namespace::write_array";
-    // suppress unused var warning
-    (void) function__;
-    try {
-      std::vector<double> alpha =
-        std::vector<double>(alpha_1dim__,
-          std::numeric_limits<double>::quiet_NaN());
-      current_statement__ = 1;
-      alpha = in__.template read<std::vector<local_scalar_t__>>(alpha_1dim__);
-      std::vector<double> mort_coef =
-        std::vector<double>(mort_coef_1dim__,
-          std::numeric_limits<double>::quiet_NaN());
-      current_statement__ = 2;
-      mort_coef = in__.template read<
-                    std::vector<local_scalar_t__>>(mort_coef_1dim__);
-      std::vector<std::vector<double>> logit_hzd =
-        std::vector<std::vector<double>>(time_groups,
-          std::vector<double>(strata_groups,
-            std::numeric_limits<double>::quiet_NaN()));
-      current_statement__ = 3;
-      logit_hzd = in__.template read<
-                    std::vector<std::vector<local_scalar_t__>>>(time_groups,
-                    strata_groups);
-      double active_detection = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 4;
-      active_detection = in__.template read_constrain_lub<local_scalar_t__,
-                           jacobian__>(0, 1, lp__);
-      double passive_asymptomatic_detection =
-        std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 5;
-      passive_asymptomatic_detection = in__.template read_constrain_lub<
-                                         local_scalar_t__, jacobian__>(0, 1,
-                                         lp__);
-      double passive_symptomatic_detection =
-        std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 6;
-      passive_symptomatic_detection = in__.template read_constrain_lub<
-                                        local_scalar_t__, jacobian__>(0, 1,
-                                        lp__);
-      std::vector<double> xi =
-        std::vector<double>(strata_groups,
-          std::numeric_limits<double>::quiet_NaN());
-      std::vector<double> mortality =
-        std::vector<double>(strata_groups,
-          std::numeric_limits<double>::quiet_NaN());
-      std::vector<std::vector<double>> S =
-        std::vector<std::vector<double>>(time_groups,
-          std::vector<double>(strata_groups,
-            std::numeric_limits<double>::quiet_NaN()));
-      std::vector<std::vector<double>> C =
-        std::vector<std::vector<double>>(time_groups,
-          std::vector<double>(strata_groups,
-            std::numeric_limits<double>::quiet_NaN()));
-      std::vector<double> theta =
-        std::vector<double>(strata_groups,
-          std::numeric_limits<double>::quiet_NaN());
-      std::vector<double> passive_denom =
-        std::vector<double>(strata_groups,
-          std::numeric_limits<double>::quiet_NaN());
-      double xi_tmp = std::numeric_limits<double>::quiet_NaN();
-      double mort_tmp = std::numeric_limits<double>::quiet_NaN();
-      out__.write(alpha);
-      out__.write(mort_coef);
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          out__.write(logit_hzd[(sym2__ - 1)][(sym1__ - 1)]);
-        }
-      }
-      out__.write(active_detection);
-      out__.write(passive_asymptomatic_detection);
-      out__.write(passive_symptomatic_detection);
-      if (stan::math::logical_negation(
-            (stan::math::primitive_value(emit_transformed_parameters__) ||
-            stan::math::primitive_value(emit_generated_quantities__)))) {
-        return ;
-      }
-      current_statement__ = 29;
-      for (int i = 1; i <= strata_groups; ++i) {
-        current_statement__ = 18;
-        xi_tmp = stan::model::rvalue(alpha, "alpha",
-                   stan::model::index_uni(1));
-        current_statement__ = 19;
-        mort_tmp = stan::model::rvalue(mort_coef, "mort_coef",
-                     stan::model::index_uni(1));
-        current_statement__ = 23;
-        for (int j = 1; j <= degrees_of_freedom; ++j) {
-          current_statement__ = 20;
-          xi_tmp = (xi_tmp +
-            (stan::model::rvalue(alpha, "alpha",
-               stan::model::index_uni((j + 1)))
-            * stan::math::pow((i - 1.0), j)));
-          current_statement__ = 21;
-          mort_tmp = (mort_tmp +
-            (stan::model::rvalue(mort_coef, "mort_coef",
-               stan::model::index_uni((j + 1)))
-            * stan::math::pow((i - 1.0), j)));
-        }
-        current_statement__ = 24;
-        stan::model::assign(xi, stan::math::inv_logit(xi_tmp),
-          "assigning variable xi", stan::model::index_uni(i));
-        current_statement__ = 25;
-        stan::model::assign(mortality, stan::math::inv_logit(mort_tmp),
-          "assigning variable mortality", stan::model::index_uni(i));
-        current_statement__ = 26;
-        stan::model::assign(S,
-          stan::model::rvalue(population, "population",
-            stan::model::index_uni(i)), "assigning variable S",
-          stan::model::index_uni(1), stan::model::index_uni(i));
-        current_statement__ = 27;
-        stan::model::assign(C,
-          (stan::model::rvalue(population, "population",
-             stan::model::index_uni(i)) *
-          stan::math::inv_logit(
-            stan::model::rvalue(logit_hzd, "logit_hzd",
-              stan::model::index_uni(1), stan::model::index_uni(i)))),
-          "assigning variable C", stan::model::index_uni(1),
-          stan::model::index_uni(i));
-      }
-      current_statement__ = 35;
-      for (int i = 2; i <= time_groups; ++i) {
-        current_statement__ = 33;
-        for (int j = 1; j <= strata_groups; ++j) {
-          current_statement__ = 30;
-          stan::model::assign(S,
-            (stan::model::rvalue(S, "S", stan::model::index_uni((i - 1)),
-               stan::model::index_uni(j)) -
-            stan::model::rvalue(C, "C", stan::model::index_uni((i - 1)),
-              stan::model::index_uni(j))), "assigning variable S",
-            stan::model::index_uni(i), stan::model::index_uni(j));
-          current_statement__ = 31;
-          stan::model::assign(C,
-            (stan::model::rvalue(S, "S", stan::model::index_uni(i),
-               stan::model::index_uni(j)) *
-            stan::math::inv_logit(
-              stan::model::rvalue(logit_hzd, "logit_hzd",
-                stan::model::index_uni(i), stan::model::index_uni(j)))),
-            "assigning variable C", stan::model::index_uni(i),
-            stan::model::index_uni(j));
-        }
-      }
-      current_statement__ = 39;
-      for (int i = 1; i <= strata_groups; ++i) {
-        current_statement__ = 36;
-        stan::model::assign(theta, ((passive_asymptomatic_detection * (1.0 -
-          stan::model::rvalue(xi, "xi", stan::model::index_uni(i)))) +
-          (passive_symptomatic_detection *
-          stan::model::rvalue(xi, "xi", stan::model::index_uni(i)))),
-          "assigning variable theta", stan::model::index_uni(i));
-        current_statement__ = 37;
-        stan::model::assign(passive_denom, (1.0 - ((1.0 -
-          stan::model::rvalue(mortality, "mortality",
-            stan::model::index_uni(i))) * (1.0 -
-          stan::model::rvalue(theta, "theta", stan::model::index_uni(i))))),
-          "assigning variable passive_denom", stan::model::index_uni(i));
-      }
-      current_statement__ = 7;
-      stan::math::check_greater_or_equal(function__, "xi", xi, 0);
-      current_statement__ = 7;
-      stan::math::check_less_or_equal(function__, "xi", xi, 1);
-      current_statement__ = 8;
-      stan::math::check_greater_or_equal(function__, "mortality", mortality,
-        0);
-      current_statement__ = 8;
-      stan::math::check_less_or_equal(function__, "mortality", mortality, 1);
-      current_statement__ = 9;
-      stan::math::check_greater_or_equal(function__, "S", S, 0);
-      current_statement__ = 10;
-      stan::math::check_greater_or_equal(function__, "C", C, 0);
-      current_statement__ = 11;
-      stan::math::check_greater_or_equal(function__, "theta", theta, 0);
-      current_statement__ = 12;
-      stan::math::check_greater_or_equal(function__, "passive_denom",
-        passive_denom, 0);
-      if (emit_transformed_parameters__) {
-        out__.write(xi);
-        out__.write(mortality);
-        for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-          for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-            out__.write(S[(sym2__ - 1)][(sym1__ - 1)]);
-          }
-        }
-        for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-          for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-            out__.write(C[(sym2__ - 1)][(sym1__ - 1)]);
-          }
-        }
-        out__.write(theta);
-        out__.write(passive_denom);
-        out__.write(xi_tmp);
-        out__.write(mort_tmp);
-      }
-      if (stan::math::logical_negation(emit_generated_quantities__)) {
-        return ;
-      }
-      std::vector<std::vector<int>> C_active_additional =
-        std::vector<std::vector<int>>(time_groups,
-          std::vector<int>(strata_groups, std::numeric_limits<int>::min()));
-      std::vector<std::vector<int>> C_passive_additional =
-        std::vector<std::vector<int>>(time_groups,
-          std::vector<int>(strata_groups, std::numeric_limits<int>::min()));
-      double gq_tmp = std::numeric_limits<double>::quiet_NaN();
-      current_statement__ = 55;
-      for (int i = 1; i <= time_groups; ++i) {
-        current_statement__ = 53;
-        for (int j = 1; j <= strata_groups; ++j) {
-          current_statement__ = 40;
-          gq_tmp = (stan::model::rvalue(C, "C", stan::model::index_uni(i),
-                      stan::model::index_uni(j))
-            -
-            stan::model::rvalue(I_active, "I_active",
-              stan::model::index_uni(i), stan::model::index_uni(j)));
-          current_statement__ = 45;
-          if (stan::math::logical_gt(gq_tmp, 0)) {
-            current_statement__ = 43;
-            stan::model::assign(C_active_additional,
-              stan::math::poisson_rng((active_detection * gq_tmp), base_rng__),
-              "assigning variable C_active_additional",
-              stan::model::index_uni(i), stan::model::index_uni(j));
-          } else {
-            current_statement__ = 41;
-            stan::model::assign(C_active_additional, 0,
-              "assigning variable C_active_additional",
-              stan::model::index_uni(i), stan::model::index_uni(j));
-          }
-          current_statement__ = 46;
-          gq_tmp = (((stan::model::rvalue(C, "C", stan::model::index_uni(i),
-                        stan::model::index_uni(j))
-            -
-            stan::model::rvalue(I_passive, "I_passive",
-              stan::model::index_uni(i), stan::model::index_uni(j)))
-            -
-            stan::model::rvalue(I_active, "I_active",
-              stan::model::index_uni(i), stan::model::index_uni(j)))
-            -
-            stan::model::rvalue(C_active_additional, "C_active_additional",
-              stan::model::index_uni(i), stan::model::index_uni(j)));
-          current_statement__ = 51;
-          if (stan::math::logical_gt(gq_tmp, 0)) {
-            current_statement__ = 49;
-            stan::model::assign(C_passive_additional,
-              stan::math::poisson_rng(gq_tmp, base_rng__),
-              "assigning variable C_passive_additional",
-              stan::model::index_uni(i), stan::model::index_uni(j));
-          } else {
-            current_statement__ = 47;
-            stan::model::assign(C_passive_additional, 0,
-              "assigning variable C_passive_additional",
-              stan::model::index_uni(i), stan::model::index_uni(j));
-          }
-        }
-      }
-      current_statement__ = 15;
-      stan::math::check_greater_or_equal(function__, "C_active_additional",
-        C_active_additional, 0);
-      current_statement__ = 16;
-      stan::math::check_greater_or_equal(function__, "C_passive_additional",
-        C_passive_additional, 0);
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          out__.write(C_active_additional[(sym2__ - 1)][(sym1__ - 1)]);
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          out__.write(C_passive_additional[(sym2__ - 1)][(sym1__ - 1)]);
-        }
-      }
-      out__.write(gq_tmp);
-    } catch (const std::exception& e) {
-      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
-    }
-  }
-  template <typename VecVar, typename VecI,
-            stan::require_vector_t<VecVar>* = nullptr,
-            stan::require_vector_like_vt<std::is_integral, VecI>* = nullptr>
-  inline void
-  unconstrain_array_impl(const VecVar& params_r__, const VecI& params_i__,
-                         VecVar& vars__, std::ostream* pstream__ = nullptr) const {
-    using local_scalar_t__ = double;
-    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
-    stan::io::serializer<local_scalar_t__> out__(vars__);
-    int current_statement__ = 0;
-    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
-    // suppress unused var warning
-    (void) DUMMY_VAR__;
-    try {
-      int pos__ = std::numeric_limits<int>::min();
-      pos__ = 1;
-      std::vector<local_scalar_t__> alpha =
-        std::vector<local_scalar_t__>(alpha_1dim__, DUMMY_VAR__);
-      current_statement__ = 1;
-      stan::model::assign(alpha,
-        in__.read<std::vector<local_scalar_t__>>(alpha_1dim__),
-        "assigning variable alpha");
-      out__.write(alpha);
-      std::vector<local_scalar_t__> mort_coef =
-        std::vector<local_scalar_t__>(mort_coef_1dim__, DUMMY_VAR__);
-      current_statement__ = 2;
-      stan::model::assign(mort_coef,
-        in__.read<std::vector<local_scalar_t__>>(mort_coef_1dim__),
-        "assigning variable mort_coef");
-      out__.write(mort_coef);
-      std::vector<std::vector<local_scalar_t__>> logit_hzd =
-        std::vector<std::vector<local_scalar_t__>>(time_groups,
-          std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__));
-      current_statement__ = 3;
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        current_statement__ = 3;
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          current_statement__ = 3;
-          logit_hzd[(sym2__ - 1)][(sym1__ -
-            1)] = in__.read<local_scalar_t__>();
-        }
-      }
-      out__.write(logit_hzd);
-      local_scalar_t__ active_detection = DUMMY_VAR__;
-      current_statement__ = 4;
-      active_detection = in__.read<local_scalar_t__>();
-      out__.write_free_lub(0, 1, active_detection);
-      local_scalar_t__ passive_asymptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 5;
-      passive_asymptomatic_detection = in__.read<local_scalar_t__>();
-      out__.write_free_lub(0, 1, passive_asymptomatic_detection);
-      local_scalar_t__ passive_symptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 6;
-      passive_symptomatic_detection = in__.read<local_scalar_t__>();
-      out__.write_free_lub(0, 1, passive_symptomatic_detection);
-    } catch (const std::exception& e) {
-      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
-    }
-  }
-  template <typename VecVar, stan::require_vector_t<VecVar>* = nullptr>
-  inline void
-  transform_inits_impl(const stan::io::var_context& context__, VecVar&
-                       vars__, std::ostream* pstream__ = nullptr) const {
-    using local_scalar_t__ = double;
-    stan::io::serializer<local_scalar_t__> out__(vars__);
-    int current_statement__ = 0;
-    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
-    // suppress unused var warning
-    (void) DUMMY_VAR__;
-    try {
-      current_statement__ = 1;
-      context__.validate_dims("parameter initialization", "alpha", "double",
-        std::vector<size_t>{static_cast<size_t>(alpha_1dim__)});
-      current_statement__ = 2;
-      context__.validate_dims("parameter initialization", "mort_coef",
-        "double", std::vector<size_t>{static_cast<size_t>(mort_coef_1dim__)});
-      current_statement__ = 3;
-      context__.validate_dims("parameter initialization", "logit_hzd",
-        "double",
-        std::vector<size_t>{static_cast<size_t>(time_groups),
-          static_cast<size_t>(strata_groups)});
-      current_statement__ = 4;
-      context__.validate_dims("parameter initialization", "active_detection",
-        "double", std::vector<size_t>{});
-      current_statement__ = 5;
-      context__.validate_dims("parameter initialization",
-        "passive_asymptomatic_detection", "double", std::vector<size_t>{});
-      current_statement__ = 6;
-      context__.validate_dims("parameter initialization",
-        "passive_symptomatic_detection", "double", std::vector<size_t>{});
-      int pos__ = std::numeric_limits<int>::min();
-      pos__ = 1;
-      std::vector<local_scalar_t__> alpha =
-        std::vector<local_scalar_t__>(alpha_1dim__, DUMMY_VAR__);
-      current_statement__ = 1;
-      alpha = context__.vals_r("alpha");
-      out__.write(alpha);
-      std::vector<local_scalar_t__> mort_coef =
-        std::vector<local_scalar_t__>(mort_coef_1dim__, DUMMY_VAR__);
-      current_statement__ = 2;
-      mort_coef = context__.vals_r("mort_coef");
-      out__.write(mort_coef);
-      std::vector<std::vector<local_scalar_t__>> logit_hzd =
-        std::vector<std::vector<local_scalar_t__>>(time_groups,
-          std::vector<local_scalar_t__>(strata_groups, DUMMY_VAR__));
-      {
-        std::vector<local_scalar_t__> logit_hzd_flat__;
-        current_statement__ = 3;
-        logit_hzd_flat__ = context__.vals_r("logit_hzd");
-        current_statement__ = 3;
-        pos__ = 1;
-        current_statement__ = 3;
-        for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-          current_statement__ = 3;
-          for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-            current_statement__ = 3;
-            stan::model::assign(logit_hzd, logit_hzd_flat__[(pos__ - 1)],
-              "assigning variable logit_hzd", stan::model::index_uni(sym2__),
-              stan::model::index_uni(sym1__));
-            current_statement__ = 3;
-            pos__ = (pos__ + 1);
-          }
-        }
-      }
-      out__.write(logit_hzd);
-      local_scalar_t__ active_detection = DUMMY_VAR__;
-      current_statement__ = 4;
-      active_detection = context__.vals_r("active_detection")[(1 - 1)];
-      out__.write_free_lub(0, 1, active_detection);
-      local_scalar_t__ passive_asymptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 5;
-      passive_asymptomatic_detection = context__.vals_r("passive_asymptomatic_detection")[(1
-        - 1)];
-      out__.write_free_lub(0, 1, passive_asymptomatic_detection);
-      local_scalar_t__ passive_symptomatic_detection = DUMMY_VAR__;
-      current_statement__ = 6;
-      passive_symptomatic_detection = context__.vals_r("passive_symptomatic_detection")[(1
-        - 1)];
-      out__.write_free_lub(0, 1, passive_symptomatic_detection);
-    } catch (const std::exception& e) {
-      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
-    }
-  }
-  inline void
-  get_param_names(std::vector<std::string>& names__, const bool
-                  emit_transformed_parameters__ = true, const bool
-                  emit_generated_quantities__ = true) const {
-    names__ = std::vector<std::string>{"alpha", "mort_coef", "logit_hzd",
-                "active_detection", "passive_asymptomatic_detection",
-                "passive_symptomatic_detection"};
-    if (emit_transformed_parameters__) {
-      std::vector<std::string>
-        temp{"xi", "mortality", "S", "C", "theta", "passive_denom", "xi_tmp",
-             "mort_tmp"};
-      names__.reserve(names__.size() + temp.size());
-      names__.insert(names__.end(), temp.begin(), temp.end());
-    }
-    if (emit_generated_quantities__) {
-      std::vector<std::string>
-        temp{"C_active_additional", "C_passive_additional", "gq_tmp"};
-      names__.reserve(names__.size() + temp.size());
-      names__.insert(names__.end(), temp.begin(), temp.end());
-    }
-  }
-  inline void
-  get_dims(std::vector<std::vector<size_t>>& dimss__, const bool
-           emit_transformed_parameters__ = true, const bool
-           emit_generated_quantities__ = true) const {
-    dimss__ = std::vector<std::vector<size_t>>{std::vector<size_t>{static_cast<
-                                                                    size_t>(
-                                                                    alpha_1dim__)},
-                std::vector<size_t>{static_cast<size_t>(mort_coef_1dim__)},
-                std::vector<size_t>{static_cast<size_t>(time_groups),
-                  static_cast<size_t>(strata_groups)}, std::vector<size_t>{},
-                std::vector<size_t>{}, std::vector<size_t>{}};
-    if (emit_transformed_parameters__) {
-      std::vector<std::vector<size_t>>
-        temp{std::vector<size_t>{static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(time_groups),
-               static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(time_groups),
-               static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{}, std::vector<size_t>{}};
-      dimss__.reserve(dimss__.size() + temp.size());
-      dimss__.insert(dimss__.end(), temp.begin(), temp.end());
-    }
-    if (emit_generated_quantities__) {
-      std::vector<std::vector<size_t>>
-        temp{std::vector<size_t>{static_cast<size_t>(time_groups),
-               static_cast<size_t>(strata_groups)},
-             std::vector<size_t>{static_cast<size_t>(time_groups),
-               static_cast<size_t>(strata_groups)}, std::vector<size_t>{}};
-      dimss__.reserve(dimss__.size() + temp.size());
-      dimss__.insert(dimss__.end(), temp.begin(), temp.end());
-    }
-  }
-  inline void
-  constrained_param_names(std::vector<std::string>& param_names__, bool
-                          emit_transformed_parameters__ = true, bool
-                          emit_generated_quantities__ = true) const final {
-    for (int sym1__ = 1; sym1__ <= alpha_1dim__; ++sym1__) {
-      param_names__.emplace_back(std::string() + "alpha" + '.' +
-        std::to_string(sym1__));
-    }
-    for (int sym1__ = 1; sym1__ <= mort_coef_1dim__; ++sym1__) {
-      param_names__.emplace_back(std::string() + "mort_coef" + '.' +
-        std::to_string(sym1__));
-    }
-    for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-      for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-        param_names__.emplace_back(std::string() + "logit_hzd" + '.' +
-          std::to_string(sym2__) + '.' + std::to_string(sym1__));
-      }
-    }
-    param_names__.emplace_back(std::string() + "active_detection");
-    param_names__.emplace_back(std::string() +
-      "passive_asymptomatic_detection");
-    param_names__.emplace_back(std::string() +
-      "passive_symptomatic_detection");
-    if (emit_transformed_parameters__) {
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "xi" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "mortality" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "S" + '.' +
-            std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C" + '.' +
-            std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "theta" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "passive_denom" + '.' +
-          std::to_string(sym1__));
-      }
-      param_names__.emplace_back(std::string() + "xi_tmp");
-      param_names__.emplace_back(std::string() + "mort_tmp");
-    }
-    if (emit_generated_quantities__) {
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C_active_additional" +
-            '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C_passive_additional" +
-            '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      param_names__.emplace_back(std::string() + "gq_tmp");
-    }
-  }
-  inline void
-  unconstrained_param_names(std::vector<std::string>& param_names__, bool
-                            emit_transformed_parameters__ = true, bool
-                            emit_generated_quantities__ = true) const final {
-    for (int sym1__ = 1; sym1__ <= alpha_1dim__; ++sym1__) {
-      param_names__.emplace_back(std::string() + "alpha" + '.' +
-        std::to_string(sym1__));
-    }
-    for (int sym1__ = 1; sym1__ <= mort_coef_1dim__; ++sym1__) {
-      param_names__.emplace_back(std::string() + "mort_coef" + '.' +
-        std::to_string(sym1__));
-    }
-    for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-      for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-        param_names__.emplace_back(std::string() + "logit_hzd" + '.' +
-          std::to_string(sym2__) + '.' + std::to_string(sym1__));
-      }
-    }
-    param_names__.emplace_back(std::string() + "active_detection");
-    param_names__.emplace_back(std::string() +
-      "passive_asymptomatic_detection");
-    param_names__.emplace_back(std::string() +
-      "passive_symptomatic_detection");
-    if (emit_transformed_parameters__) {
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "xi" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "mortality" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "S" + '.' +
-            std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C" + '.' +
-            std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "theta" + '.' +
-          std::to_string(sym1__));
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        param_names__.emplace_back(std::string() + "passive_denom" + '.' +
-          std::to_string(sym1__));
-      }
-      param_names__.emplace_back(std::string() + "xi_tmp");
-      param_names__.emplace_back(std::string() + "mort_tmp");
-    }
-    if (emit_generated_quantities__) {
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C_active_additional" +
-            '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      for (int sym1__ = 1; sym1__ <= strata_groups; ++sym1__) {
-        for (int sym2__ = 1; sym2__ <= time_groups; ++sym2__) {
-          param_names__.emplace_back(std::string() + "C_passive_additional" +
-            '.' + std::to_string(sym2__) + '.' + std::to_string(sym1__));
-        }
-      }
-      param_names__.emplace_back(std::string() + "gq_tmp");
-    }
-  }
-  inline std::string get_constrained_sizedtypes() const {
-    return std::string("[{\"name\":\"alpha\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(alpha_1dim__) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"parameters\"},{\"name\":\"mort_coef\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(mort_coef_1dim__) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"parameters\"},{\"name\":\"logit_hzd\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"parameters\"},{\"name\":\"active_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"passive_asymptomatic_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"passive_symptomatic_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"xi\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"mortality\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"S\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"transformed_parameters\"},{\"name\":\"C\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"transformed_parameters\"},{\"name\":\"theta\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"passive_denom\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"xi_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"mort_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"C_active_additional\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"C_passive_additional\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"gq_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]");
-  }
-  inline std::string get_unconstrained_sizedtypes() const {
-    return std::string("[{\"name\":\"alpha\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(alpha_1dim__) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"parameters\"},{\"name\":\"mort_coef\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(mort_coef_1dim__) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"parameters\"},{\"name\":\"logit_hzd\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"parameters\"},{\"name\":\"active_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"passive_asymptomatic_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"passive_symptomatic_detection\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"xi\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"mortality\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"S\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"transformed_parameters\"},{\"name\":\"C\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}}},\"block\":\"transformed_parameters\"},{\"name\":\"theta\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"passive_denom\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"real\"}},\"block\":\"transformed_parameters\"},{\"name\":\"xi_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"mort_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"transformed_parameters\"},{\"name\":\"C_active_additional\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"C_passive_additional\",\"type\":{\"name\":\"array\",\"length\":" + std::to_string(time_groups) + ",\"element_type\":{\"name\":\"array\",\"length\":" + std::to_string(strata_groups) + ",\"element_type\":{\"name\":\"int\"}}},\"block\":\"generated_quantities\"},{\"name\":\"gq_tmp\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]");
-  }
-  // Begin method overload boilerplate
-  template <typename RNG> inline void
-  write_array(RNG& base_rng, Eigen::Matrix<double,-1,1>& params_r,
-              Eigen::Matrix<double,-1,1>& vars, const bool
-              emit_transformed_parameters = true, const bool
-              emit_generated_quantities = true, std::ostream*
-              pstream = nullptr) const {
-    const size_t num_params__ = (((((alpha_1dim__ + mort_coef_1dim__) +
-      (time_groups * strata_groups)) + 1) + 1) + 1);
-    const size_t num_transformed = emit_transformed_parameters *
-      ((((((((strata_groups + strata_groups) + (time_groups * strata_groups))
-      + (time_groups * strata_groups)) + strata_groups) + strata_groups) + 1)
-      + 1));
-    const size_t num_gen_quantities = emit_generated_quantities *
-      ((((time_groups * strata_groups) + (time_groups * strata_groups)) + 1));
-    const size_t num_to_write = num_params__ + num_transformed +
-      num_gen_quantities;
-    std::vector<int> params_i;
-    vars = Eigen::Matrix<double,-1,1>::Constant(num_to_write,
-             std::numeric_limits<double>::quiet_NaN());
-    write_array_impl(base_rng, params_r, params_i, vars,
-      emit_transformed_parameters, emit_generated_quantities, pstream);
-  }
-  template <typename RNG> inline void
-  write_array(RNG& base_rng, std::vector<double>& params_r, std::vector<int>&
-              params_i, std::vector<double>& vars, bool
-              emit_transformed_parameters = true, bool
-              emit_generated_quantities = true, std::ostream*
-              pstream = nullptr) const {
-    const size_t num_params__ = (((((alpha_1dim__ + mort_coef_1dim__) +
-      (time_groups * strata_groups)) + 1) + 1) + 1);
-    const size_t num_transformed = emit_transformed_parameters *
-      ((((((((strata_groups + strata_groups) + (time_groups * strata_groups))
-      + (time_groups * strata_groups)) + strata_groups) + strata_groups) + 1)
-      + 1));
-    const size_t num_gen_quantities = emit_generated_quantities *
-      ((((time_groups * strata_groups) + (time_groups * strata_groups)) + 1));
-    const size_t num_to_write = num_params__ + num_transformed +
-      num_gen_quantities;
-    vars = std::vector<double>(num_to_write,
-             std::numeric_limits<double>::quiet_NaN());
-    write_array_impl(base_rng, params_r, params_i, vars,
-      emit_transformed_parameters, emit_generated_quantities, pstream);
-  }
-  template <bool propto__, bool jacobian__, typename T_> inline T_
-  log_prob(Eigen::Matrix<T_,-1,1>& params_r, std::ostream* pstream = nullptr) const {
-    Eigen::Matrix<int,-1,1> params_i;
-    return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
-  }
-  template <bool propto__, bool jacobian__, typename T_> inline T_
-  log_prob(std::vector<T_>& params_r, std::vector<int>& params_i,
-           std::ostream* pstream = nullptr) const {
-    return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
-  }
-  inline void
-  transform_inits(const stan::io::var_context& context,
-                  Eigen::Matrix<double,-1,1>& params_r, std::ostream*
-                  pstream = nullptr) const final {
-    std::vector<double> params_r_vec(params_r.size());
-    std::vector<int> params_i;
-    transform_inits(context, params_i, params_r_vec, pstream);
-    params_r = Eigen::Map<Eigen::Matrix<double,-1,1>>(params_r_vec.data(),
-                 params_r_vec.size());
-  }
-  inline void
-  transform_inits(const stan::io::var_context& context, std::vector<int>&
-                  params_i, std::vector<double>& vars, std::ostream*
-                  pstream__ = nullptr) const {
-    vars.resize(num_params_r__);
-    transform_inits_impl(context, vars, pstream__);
-  }
-  inline void
-  unconstrain_array(const std::vector<double>& params_constrained,
-                    std::vector<double>& params_unconstrained, std::ostream*
-                    pstream = nullptr) const {
-    const std::vector<int> params_i;
-    params_unconstrained = std::vector<double>(num_params_r__,
-                             std::numeric_limits<double>::quiet_NaN());
-    unconstrain_array_impl(params_constrained, params_i,
-      params_unconstrained, pstream);
-  }
-  inline void
-  unconstrain_array(const Eigen::Matrix<double,-1,1>& params_constrained,
-                    Eigen::Matrix<double,-1,1>& params_unconstrained,
-                    std::ostream* pstream = nullptr) const {
-    const std::vector<int> params_i;
-    params_unconstrained = Eigen::Matrix<double,-1,1>::Constant(num_params_r__,
-                             std::numeric_limits<double>::quiet_NaN());
-    unconstrain_array_impl(params_constrained, params_i,
-      params_unconstrained, pstream);
-  }
-};
-}
-using stan_model = model_severity_estimate_by_surveillance_namespace::model_severity_estimate_by_surveillance;
-#ifndef USING_R
-// Boilerplate
-stan::model::model_base&
-new_model(stan::io::var_context& data_context, unsigned int seed,
-          std::ostream* msg_stream) {
-  stan_model* m = new stan_model(data_context, seed, msg_stream);
-  return *m;
-}
-stan::math::profile_map& get_stan_profile_data() {
-  return model_severity_estimate_by_surveillance_namespace::profiles__;
-}
-#endif
-#endif
diff --git a/tests/testthat/test-estimate_severity.R b/tests/testthat/test-estimate_severity.R
deleted file mode 100644
index 44de73b..0000000
--- a/tests/testthat/test-estimate_severity.R
+++ /dev/null
@@ -1,566 +0,0 @@
-test_that("Input Validation", {
-  # Ideal valid inputs
-  linelist <- data.frame(
-    patient_id = letters,
-    week = rep_len(1L:3L, 26L),
-    sex = rep_len(c("Female", "Male"), 26L),
-    testing_type = rep_len(c("A", "A", "A", "P", "P"), 26L),
-    patient_status = rep_len(c("A", "D", "S", "S"), 26L)
-  )
-  population <- data.frame(
-    sex = c("Female", "Male"),
-    value = c(4000L, 3975L)
-  )
-  # linelist param
-  expect_error(
-    estimate_severity(
-      NULL,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "`linelist` is not 'data.frame' like.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      list(abc = letters, def = LETTERS),
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "`linelist` is not 'data.frame' like.",
-    fixed = TRUE
-  )
-  # population param
-  expect_error(
-    estimate_severity(
-      linelist,
-      NULL,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "`population` is not 'data.frame' like.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      list(abc = letters, def = LETTERS),
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "`population` is not 'data.frame' like.",
-    fixed = TRUE
-  )
-  # surveillance param
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      1L,
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "Assertion on 'surveillance' failed: ",
-      "Must be of type 'string', not 'integer'."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      letters,
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "Assertion on 'surveillance' failed: Must have length 1.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      NA_character_,
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "Assertion on 'surveillance' failed: May not be NA.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "not present in linelist",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "`linelist` is missing required string columns: not present in linelist."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "week",
-      "patient_status",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "The 'week' column of `linelist` is not a character or factor, ",
-      "instead is: integer."
-    ),
-    fixed = TRUE
-  )
-  # outcome param
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      1L,
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "Assertion on 'outcome' failed: Must be of type 'string', not 'integer'."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      letters,
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "Assertion on 'outcome' failed: Must have length 1.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      NA_character_,
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = "Assertion on 'outcome' failed: May not be NA.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "not present in linelist",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "`linelist` is missing required string columns: not present in linelist."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "week",
-      time_period = "week",
-      strata = "sex"
-    ),
-    regexp = paste0(
-      "The 'week' column of `linelist` is not a character or factor, ",
-      "instead is: integer."
-    ),
-    fixed = TRUE
-  )
-  # hazard_std param
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      hazard_std = "abc"
-    ),
-    regexp = paste0(
-      "Assertion on 'hazard_std' failed: ",
-      "Must be of type 'number', not 'character'."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      hazard_std = 1:5
-    ),
-    regexp = "Assertion on 'hazard_std' failed: Must have length 1.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      hazard_std = NA_real_
-    ),
-    regexp = "Assertion on 'hazard_std' failed: May not be NA.",
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      hazard_std = -pi
-    ),
-    regexp = paste0(
-      "Assertion on 'hazard_std' failed: Element 1 is not >= 2.22045e-16."
-    ),
-    fixed = TRUE
-  )
-  # degrees_of_freedom param
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      degrees_of_freedom = "abc"
-    ),
-    regexp = paste0(
-      "Assertion on 'degrees_of_freedom' failed: ",
-      "Must be of type 'integerish', not 'character'."
-    ),
-    fixed = TRUE
-  )
-  expect_error(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      degrees_of_freedom = pi
-    ),
-    regexp = paste0(
-      "Assertion on 'degrees_of_freedom' failed: ",
-      "Must be of type 'integerish', ",
-      "but element 1 is not close to an integer."
-    ),
-    fixed = TRUE
-  )
-})
-
-test_that("Output Validation", {
-  # Slow, skip if not manual run
-  skip_on_cran()
-  # Sample data
-  linelist <- data.frame(
-    patient_id = letters,
-    week = rep_len(1L:3L, 26L),
-    sex = rep_len(c("Female", "Male"), 26L),
-    testing_type = rep_len(c("A", "A", "A", "P", "P"), 26L),
-    patient_status = rep_len(c("A", "D", "S", "S"), 26L)
-  )
-  population <- data.frame(
-    sex = c("Female", "Male"),
-    value = c(4000L, 3975L)
-  )
-  # Categories inferred from linelist/population
-  severity_est <- suppressWarnings(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      # Params passed to rstan::stan
-      chains = 1L,
-      iter = 100L,
-      seed = 1L,
-      cores = 1L,
-      open_progress = FALSE,
-      refresh = 0L
-    )
-  )
-  # Expectations on output structure
-  expect_s4_class(severity_est, "SeverityEstimateFit")
-  expect_equal(
-    slotNames(severity_est),
-    c(
-      "model_fit",
-      "population",
-      "incidence",
-      "time_period",
-      "strata",
-      "surveillance",
-      "outcome"
-    )
-  )
-  expect_s4_class(severity_est@model_fit, "stanfit")
-  expected_population <- array(
-    data = c(4000L, 3975L),
-    dimnames = list("strata" = seq_len(2))
-  )
-  expect_identical(severity_est@population, expected_population)
-  expected_incidence <- array(
-    data = c(
-      2L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      1L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      2L,
-      2L,
-      0L,
-      0L,
-      0L,
-      1L,
-      1L,
-      0L,
-      1L,
-      2L,
-      1L,
-      1L,
-      1L,
-      1L,
-      1L,
-      0L,
-      1L,
-      1L,
-      1L,
-      1L
-    ),
-    dim = c(3L, 2L, 2L, 3L),
-    dimnames = list(
-      "time_period" = seq_len(3L),
-      "strata" = seq_len(2L),
-      "surveillance" = seq_len(2L),
-      "outcome" = seq_len(3L)
-    )
-  )
-  expect_identical(severity_est@incidence, expected_incidence)
-  expected_time_period <- data.frame(week = 1L:3L)
-  expect_identical(severity_est@time_period, expected_time_period)
-  expected_strata <- data.frame(sex = c("Female", "Male"))
-  expect_identical(severity_est@strata, expected_strata)
-  expected_surveillance <- data.frame(
-    testing_type = c("Active", "Passive"),
-    stringsAsFactors = TRUE
-  )
-  expect_identical(severity_est@surveillance, expected_surveillance)
-  expected_outcome <- data.frame(
-    patient_status = c("Asymptomatic", "Death", "Symptomatic"),
-    stringsAsFactors = TRUE
-  )
-  expect_identical(severity_est@outcome, expected_outcome)
-  # Time period/strata categories are user specified
-  severity_est <- suppressMessages(
-    estimate_severity(
-      linelist,
-      population,
-      "testing_type",
-      "patient_status",
-      time_period = "week",
-      strata = "sex",
-      time_period_reference = 0L:4L,
-      strata_reference = c("Female", "Male", "NA"),
-      # Params passed to rstan::stan
-      chains = 1L,
-      iter = 100L,
-      seed = 1L,
-      cores = 1L,
-      open_progress = FALSE,
-      refresh = 0L
-    )
-  )
-  # Expectations on output structure
-  expect_s4_class(severity_est, "SeverityEstimateFit")
-  expect_equal(
-    slotNames(severity_est),
-    c(
-      "model_fit",
-      "population",
-      "incidence",
-      "time_period",
-      "strata",
-      "surveillance",
-      "outcome"
-    )
-  )
-  expect_s4_class(severity_est@model_fit, "stanfit")
-  expected_population <- array(
-    data = c(4000L, 3975L, 0L),
-    dimnames = list("strata" = seq_len(3L))
-  )
-  expect_identical(severity_est@population, expected_population)
-  expected_incidence <- array(
-    data = c(
-      0L,
-      2L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      2L,
-      2L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      2L,
-      1L,
-      0L,
-      0L,
-      1L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      1L,
-      0L,
-      1L,
-      0L,
-      0L,
-      1L,
-      1L,
-      1L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L,
-      0L
-    ),
-    dim = c(5L, 3L, 2L, 3L),
-    dimnames = list(
-      "time_period" = seq_len(5L),
-      "strata" = seq_len(3L),
-      "surveillance" = seq_len(2L),
-      "outcome" = seq_len(3L)
-    )
-  )
-  expect_identical(severity_est@incidence, expected_incidence)
-  expected_time_period <- data.frame(week = 0L:4L)
-  expect_identical(severity_est@time_period, expected_time_period)
-  expected_strata <- data.frame(sex = c("Female", "Male", "NA"))
-  expect_identical(severity_est@strata, expected_strata)
-  expected_surveillance <- data.frame(
-    testing_type = c("Active", "Passive"),
-    stringsAsFactors = TRUE
-  )
-  expect_identical(severity_est@surveillance, expected_surveillance)
-  expected_outcome <- data.frame(
-    patient_status = c("Asymptomatic", "Death", "Symptomatic"),
-    stringsAsFactors = TRUE
-  )
-  expect_identical(severity_est@outcome, expected_outcome)
-})
diff --git a/tests/testthat/test-stan_model.R b/tests/testthat/test-stan_model.R
index 5eb8a72..d39a48f 100644
--- a/tests/testthat/test-stan_model.R
+++ b/tests/testthat/test-stan_model.R
@@ -1,16 +1,3 @@
-test_that("By-surveillance Stan model compiles", {
-  output <- rstan::stanc(
-    file = system.file(
-      "stan",
-      "severity_estimate_by_surveillance.stan",
-      package = "SeverityEstimate"
-    ),
-    model_name = "severity_estimate_by_surveillance",
-    verbose = FALSE
-  )
-  expect_true(output$status)
-})
-
 test_that("`stan_model()` errors for an unknown precompiled model", {
   expect_error(
     stan_model("does_not_exist"),