Source code for rsvsim.parameters

'''
Set the parameters for rsvsim.
'''

import numpy as np
import sciris as sc
from .settings import options as cvo # For setting global options
from . import misc as cvm
from . import defaults as cvd

__all__ = ['make_pars', 'reset_layer_pars', 'get_prognoses', 'get_genotype_choices', 'get_vaccine_choices']


[docs]def make_pars(set_prognoses=False, prog_by_age=True, **kwargs): ''' Create the parameters for the simulation. Typically, this function is used internally rather than called by the user; e.g. typical use would be to do sim = cv.Sim() and then inspect sim.pars, rather than calling this function directly. Args: set_prognoses (bool): whether or not to create prognoses (else, added when the population is created) prog_by_age (bool): whether or not to use age-based severity, mortality etc. kwargs (dict): any additional kwargs are interpreted as parameter names Returns: pars (dict): the parameters of the simulation ''' pars = {} # Population parameters pars['pop_size'] = 20e3 # Number of agents, i.e., people susceptible to RSV pars['pop_type'] = 'random' # What type of population data to use -- 'random' (fastest), 'synthpops' (best), 'hybrid' (compromise) pars['location'] = None # What location to load data from -- default Kenya # Simulation parameters pars['start_day'] = '2021-01-01' # Start day of the simulation pars['end_day'] = None # End day of the simulation pars['n_days'] = 365 # Number of days to run, if end_day isn't specified pars['timestep'] = 7 # Number of days to simulate per timestep, default to weekly pars['tspy'] = 52 # Number of timesteps per year, default to weekly pars['rand_seed'] = 1 # Random seed, if None, don't reset pars['verbose'] = cvo.verbose # Whether or not to display information during the run -- options are 0 (silent), 0.1 (some; default), 1 (default), 2 (everything) # Rescaling parameters pars['pop_scale'] = 1 # Factor by which to scale the population -- e.g. pop_scale=10 with pop_size=100e3 means a population of 1 million pars['scaled_pop'] = None # The total scaled population, i.e. the number of agents times the scale factor pars['rescale'] = True # Enable dynamic rescaling of the population -- starts with pop_scale=1 and scales up dynamically as the epidemic grows pars['rescale_threshold'] = 0.05 # Fraction susceptible population that will trigger rescaling if rescaling pars['rescale_factor'] = 1.2 # Factor by which the population is rescaled on each step pars['frac_susceptible'] = dict(groupA=dict(frac=1, date_recovered=-105), groupB=dict(frac=1, date_recovered=-105)) # What proportion of the population is susceptible to infection pars['prior_wave_dist'] = None # Network parameters, including vital dynamics and schools pars['sex_ratio'] = 0.5 pars['age_pregnancy'] = [20,45] # min and max age of pregnancy in years pars['gestation'] = [35,41] # min and max gestational age in weeks pars['dur_post_partum'] = 20 # weeks post-partum (not eligible to get pregnant) pars['max_age'] = 100 pars['pop_data_file'] = 'Kenya.json' pars['contacts'] = None # The number of contacts per layer; set by reset_layer_pars() below pars['dynam_layer'] = None # Which layers are dynamic; set by reset_layer_pars() below pars['beta_layer'] = None # Transmissibility per layer; set by reset_layer_pars() below # Basic disease transmission parameters pars['beta_dist'] = dict(dist='neg_binomial', par1=1.0, par2=0.45, step=0.01) # Distribution to draw individual level transmissibility; dispersion from https://www.researchsquare.com/article/rs-29548/v1 pars['viral_dist'] = dict(frac_time=0.3, load_ratio=2, high_cap=4) # The time varying viral load (transmissibility); estimated from Lescure 2020, Lancet, https://doi.org/10.1016/S1473-3099(20)30200-0 pars['beta'] = 0.0972 # Beta per symptomatic contact; absolute value, calibrated pars['seasonality'] = True pars['beta_seasonality'] = 0.65 pars['phase_shift'] = 5 pars['asymp_factor'] = 0.63 # Multiply beta by this factor for asymptomatic cases; no statistically significant difference in transmissibility: https://www.sciencedirect.com/science/article/pii/S1201971220302502 # Parameters that control settings and defaults for multi-genotype runs pars['n_genotypes'] = 2 # The number of genotypes circulating in the population. By default set to 2, group A and group B # Parameters used to calculate immunity pars['imm_init'] = 1 # Peak efficacy afforded by types of immunity pars['imm_decay'] = {'mab': dict(form='growth_decay', growth_time=14, decay_rate=0.04), 'ab': dict(form='growth_decay', growth_time=14, decay_rate=0.01)} pars['imm_kin'] = None # Constructed during sim initialization using the imm_decay parameters pars['immunity'] = None # Matrix of immunity and cross-immunity factors, set by init_immunity() in immunity.py pars['cord_transfer'] = 1.1 # Degree of immunity passed from mom to baby at time of birth pars['dur_redux'] = 0.9 # Relative reduction in duration of infection for each additional infection (gets applied to rel_trans) # genotype-specific disease transmission parameters. By default, these are set up for a single genotype, but can all be modified for multiple genotypes pars['rel_beta'] = 1.0 # Relative transmissibility varies by genotype pars['rel_imm_genotype'] = 1.0 # Relative own-immunity varies by genotype # Duration parameters: time for disease progression pars['dur'] = {} pars['dur']['exp2inf'] = dict(dist='lognormal_int', par1=4.98, par2=0.9) # Duration from exposed to infectious; see Lauer et al., https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7081172/, appendix table S2, subtracting inf2sym duration pars['dur']['inf2sym'] = dict(dist='lognormal_int', par1=1.1, par2=0.9) # Duration from infectious to symptomatic; see Linton et al., https://doi.org/10.3390/jcm9020538, from Table 2, 5.6 day incubation period - 4.5 day exp2inf from Lauer et al. pars['dur']['sym2sev'] = dict(dist='lognormal_int', par1=6.16, par2=4.9) # Duration from symptomatic to severe symptoms; see Linton et al., https://doi.org/10.3390/jcm9020538, from Table 2, 6.6 day onset to hospital admission (deceased); see also Wang et al., https://jamanetwork.com/journals/jama/fullarticle/2761044, 7 days (Table 1) pars['dur']['sev2crit'] = dict(dist='lognormal_int', par1=1.5, par2=2.0) # Duration from severe symptoms to requiring ICU; average of 1.9 and 1.0; see Chen et al., https://www.sciencedirect.com/science/article/pii/S0163445320301195, 8.5 days total - 6.6 days sym2sev = 1.9 days; see also Wang et al., https://jamanetwork.com/journals/jama/fullarticle/2761044, Table 3, 1 day, IQR 0-3 days; std=2.0 is an estimate # Duration parameters: time for disease recovery pars['dur']['asym2rec'] = dict(dist='lognormal_int', par1=3.0, par2=2.0) # Duration for asymptomatic people to recover; see Wölfel et al., https://www.nature.com/articles/s41586-020-2196-x pars['dur']['mild2rec'] = dict(dist='lognormal_int', par1=4.16, par2=2.0) # Duration for people with mild symptoms to recover; see Wölfel et al., https://www.nature.com/articles/s41586-020-2196-x pars['dur']['sev2rec'] = dict(dist='lognormal_int', par1=10.1, par2=6.3) # Duration for people with severe symptoms to recover, 24.7 days total; see Verity et al., https://www.thelancet.com/journals/laninf/article/PIIS1473-3099(20)30243-7/fulltext; 18.1 days = 24.7 onset-to-recovery - 6.6 sym2sev; 6.3 = 0.35 coefficient of variation * 18.1; see also https://doi.org/10.1017/S0950268820001259 (22 days) and https://doi.org/10.3390/ijerph17207560 (3-10 days) pars['dur']['crit2rec'] = dict(dist='lognormal_int', par1=18.1, par2=6.3) # Duration for people with critical symptoms to recover; as above (Verity et al.) pars['dur']['crit2die'] = dict(dist='lognormal_int', par1=10.7, par2=4.8) # Duration from critical symptoms to death, 18.8 days total; see Verity et al., https://www.thelancet.com/journals/laninf/article/PIIS1473-3099(20)30243-7/fulltext; 10.7 = 18.8 onset-to-death - 6.6 sym2sev - 1.5 sev2crit; 4.8 = 0.45 coefficient of variation * 10.7 # Severity parameters: probabilities of symptom progression pars['rel_symp_prob'] = 1.0 # Scale factor for proportion of symptomatic cases pars['rel_severe_prob'] = 1.0 # Scale factor for proportion of symptomatic cases that become severe pars['rel_crit_prob'] = 1.0 # Scale factor for proportion of severe cases that become critical pars['rel_death_prob'] = 1.0 # Scale factor for proportion of critical cases that result in death pars['prog_by_age'] = prog_by_age # Whether to set disease progression based on the person's age pars['prognoses'] = None # The actual arrays of prognoses by age; this is populated later # Events and interventions pars['interventions'] = [] # The interventions present in this simulation; populated by the user pars['analyzers'] = [] # Custom analysis functions; populated by the user pars['timelimit'] = None # Time limit for the simulation (seconds) pars['stopping_func'] = None # A function to call to stop the sim partway through # Health system parameters pars['n_beds_hosp'] = None # The number of hospital (adult acute care) beds available for severely ill patients (default is no constraint) pars['n_beds_icu'] = None # The number of ICU beds available for critically ill patients (default is no constraint) pars['no_hosp_factor'] = 2.0 # Multiplier for how much more likely severely ill people are to become critical if no hospital beds are available pars['no_icu_factor'] = 2.0 # Multiplier for how much more likely critically ill people are to die if no ICU beds are available # Handle vaccine and genotype parameters pars['vaccine_pars'] = {} # Vaccines that are being used; populated during initialization pars['vaccine_map'] = {} #Reverse mapping from number to vaccine key pars['genotypes'] = [] # genotypes that are circulating in sim; populated by the user, see immunity.py pars['genotype_map'] = dict() # Populated by the user, see immunity.py pars['genotype_pars'] = dict() # Populated by the user, see immunity.py # Update with any supplied parameter values and generate things that need to be generated pars.update(kwargs) reset_layer_pars(pars) if set_prognoses: # If not set here, gets set when the population is initialized pars['prognoses'] = get_prognoses() # Default to age-specific prognoses return pars
# Define which parameters need to be specified as a dictionary by layer -- define here so it's available at the module level for sim.py layer_pars = ['beta_layer', 'contacts', 'dynam_layer']
[docs]def reset_layer_pars(pars, layer_keys=None, force=False): ''' Helper function to set layer-specific parameters. If layer keys are not provided, then set them based on the population type. This function is not usually called directly by the user, although it can sometimes be used to fix layer key mismatches (i.e. if the contact layers in the population do not match the parameters). More commonly, however, mismatches need to be fixed explicitly. Args: pars (dict): the parameters dictionary layer_keys (list): the layer keys of the population, if available force (bool): reset the parameters even if they already exist ''' # Specify defaults for random -- layer 'a' for 'all' layer_defaults = {} layer_defaults['random'] = dict( beta_layer = dict(a=1.0), # Default beta contacts = dict(a=20), # Default number of contacts dynam_layer = dict(a=0), # Do not use dynamic layers by default ) # Specify defaults for hybrid -- household, school, and community layers (h, s, c) layer_defaults['hybrid'] = dict( beta_layer = dict(h=1.0, s=1.5, c=0.75), # Per-population beta weights; relative; in part based on Table S14 of https://science.sciencemag.org/content/sci/suppl/2020/04/28/science.abb8001.DC1/abb8001_Zhang_SM.pdf contacts = dict(h=2.0, s=20, c=20), # Number of contacts per person per day, estimated dynam_layer = dict(h=0, s=0, c=0), # Which layers are dynamic -- none by default ) # Specify defaults for SynthPops layer_defaults['synthpops'] = sc.dcp(layer_defaults['hybrid']) # Choose the parameter defaults based on the population type, and get the layer keys try: defaults = layer_defaults[pars['pop_type']] except Exception as E: errormsg = f'Cannot load defaults for population type "{pars["pop_type"]}": must be hybrid, random, or synthpops' raise ValueError(errormsg) from E default_layer_keys = list(defaults['beta_layer'].keys()) # All layers should be the same, but use beta_layer for convenience # Actually set the parameters for pkey in layer_pars: par = {} # Initialize this parameter default_val = layer_defaults['random'][pkey]['a'] # Get the default value for this parameter # If forcing, we overwrite any existing parameter values if force: par_dict = defaults[pkey] # Just use defaults else: par_dict = sc.mergedicts(defaults[pkey], pars.get(pkey, None)) # Use user-supplied parameters if available, else default # Figure out what the layer keys for this parameter are (may be different between parameters) if layer_keys: par_layer_keys = layer_keys # Use supplied layer keys else: par_layer_keys = list(sc.odict.fromkeys(default_layer_keys + list(par_dict.keys()))) # If not supplied, use the defaults, plus any extra from the par_dict; adapted from https://www.askpython.com/python/remove-duplicate-elements-from-list-python # Construct this parameter, layer by layer for lkey in par_layer_keys: # Loop over layers par[lkey] = par_dict.get(lkey, default_val) # Get the value for this layer if available, else use the default for random pars[pkey] = par # Save this parameter to the dictionary return
[docs]def get_prognoses(): ''' Return the default parameter values for prognoses The prognosis probabilities are conditional given the previous disease state. Returns: prog_pars (dict): the dictionary of prognosis probabilities ''' prognoses = dict( age_cutoffs=np.array([0, 1, 5, 15, 55]), # Age cutoffs (lower limits) sus_ORs=np.array([1.50, 1.00, 1.00, 1.00, 1.50]), # Odds ratios for relative susceptibility trans_ORs=np.array([1.00, 1.00, 1.00, 1.00, 1.00]), # Odds ratios for relative transmissibility comorbidities=np.array([1.00, 1.00, 1.00, 1.00, 1.00]), # Comorbidities by age -- set to 1 by default symp_probs=np.array([0.91, 0.84, 0.49, 0.25, 0.25]), # Overall probability of developing symptoms severe_probs=np.array([0.050, 0.00050, 0.00050, 0.00050, 0.0050]), # Overall probability of developing severe symptoms crit_probs=np.array([0.0003, 0.00003, 0.00003, 0.00003, 0.00003]), # Overall probability of developing critical symptoms death_probs=np.array([0.00003, 0.00003, 0.00003, 0.00003, 0.00003]), # Overall probability of dying ) prognoses = relative_prognoses(prognoses) # Convert to conditional probabilities # Check that lengths match expected_len = len(prognoses['age_cutoffs']) for key,val in prognoses.items(): this_len = len(prognoses[key]) if this_len != expected_len: # pragma: no cover errormsg = f'Lengths mismatch in prognoses: {expected_len} age bins specified, but key "{key}" has {this_len} entries' raise ValueError(errormsg) return prognoses
def relative_prognoses(prognoses): ''' Convenience function to revert absolute prognoses into relative (conditional) ones. Internally, rsvsim uses relative prognoses. ''' out = sc.dcp(prognoses) out['death_probs'] /= out['crit_probs'] # Conditional probability of dying, given critical symptoms out['crit_probs'] /= out['severe_probs'] # Conditional probability of symptoms becoming critical, given severe out['severe_probs'] /= out['symp_probs'] # Conditional probability of symptoms becoming severe, given symptomatic return out def absolute_prognoses(prognoses): ''' Convenience function to revert relative (conditional) prognoses into absolute ones. Used to convert internally used relative prognoses into more readable absolute ones. **Example**:: sim = cv.Sim() abs_progs = cv.parameters.absolute_prognoses(sim['prognoses']) ''' out = sc.dcp(prognoses) out['severe_probs'] *= out['symp_probs'] # Absolute probability of severe symptoms out['crit_probs'] *= out['severe_probs'] # Absolute probability of critical symptoms out['death_probs'] *= out['crit_probs'] # Absolute probability of dying return out #%% genotype, vaccine, and immunity parameters and functions
[docs]def get_genotype_choices(): ''' Define valid pre-defined genotype names ''' # List of choices currently available: new ones can be added to the list along with their aliases choices = { 'groupA': ['groupa', 'A'], 'groupB': ['groupb', 'B'], } mapping = {name:key for key,synonyms in choices.items() for name in synonyms} # Flip from key:value to value:key return choices, mapping
[docs]def get_vaccine_choices(): ''' Define valid pre-defined vaccine names ''' # List of choices currently available: new ones can be added to the list along with their aliases choices = { 'default': ['default', None], } mapping = {name:key for key,synonyms in choices.items() for name in synonyms} # Flip from key:value to value:key return choices, mapping
def get_genotype_pars(default=False): ''' Define the default parameters for the different genotypes ''' pars = dict( groupA = dict( rel_beta = 1.0, # Default values rel_symp_prob = 1.0, # Default values rel_severe_prob = 1.0, # Default values rel_crit_prob = 1.0, # Default values rel_death_prob = 1.0, # Default values ), groupB=dict( rel_beta=1.0, # Default values rel_symp_prob=1.0, # Default values rel_severe_prob=1.0, # Default values rel_crit_prob=1.0, # Default values rel_death_prob=1.0, # Default values ), ) if default: return pars['groupA'] else: return pars def get_cross_immunity(default=False): ''' Get the cross immunity between each genotype in a sim ''' pars = dict( groupA = dict( groupA = 1.0, # Default for own-immunity groupB = 0.5, # Assumption ), groupB = dict( groupA = 0.5, # Assumption groupB = 1.0, # Default for own-immunity ), ) if default: return pars['groupA'] else: return pars def get_vaccine_genotype_pars(default=False): ''' Define the effectiveness of each vaccine against each genotype ''' pars = dict( default = dict( groupA = 1.0, groupB = 1.0, ), ) if default: return pars['default'] else: return pars def get_vaccine_dose_pars(default=False): ''' Define the parameters for each vaccine ''' # Default vaccine NAb efficacy is nearly identical to nab_eff for infection -- only alpha_inf differs. Values derived from fit to data vaccine_nab_eff = dict( alpha_inf = 1.08, beta_inf = 0.967, alpha_symp_inf = -0.739, beta_symp_inf = 0.038, alpha_sev_symp = -0.014, beta_sev_symp = 0.079, ) pars = dict( default = dict( nab_eff = sc.dcp(vaccine_nab_eff), # Efficacy of NAbs of this vaccine nab_init = dict(dist='normal', par1=0, par2=2), # Initial distribution of NAbs nab_boost = 2, # Factor by which a dose increases existing NABs doses = 1, # Number of doses for this vaccine interval = None, # Interval between doses ), ) if default: return pars['default'] else: return pars