61 COVID Intervention Lab

In this lab you will work with the diagnosis and case isolation features of the COVID model from the previous tutorial, and use them to ask how much transmission a screening-based intervention can actually prevent. The tutorial built the machinery; this lab turns it into experiments. The specific learning objectives are to:

Read the diagnosis module and the discordant edgelist code closely enough to know exactly where an intervention acts on transmission;
Run intervention scenarios that vary the diagnosis rates and the case isolation intensity, and interpret the results against the natural history assumptions;
Define concrete next steps for extending diagnosis-based interventions, including contact tracing.

How to use your time

Steps 1 through 3 are the core of the lab and are what everyone should finish in the session. The two extensions that follow are optional. They are specified in enough detail that you can start one here and finish it at home after the course, which is how we expect most people will use them.

61.1 Setup

Once you are ready, start out by clearing your R object environment, to make sure that you do not have any objects lingering from the tutorial. This can be accomplished with:

Code

rm(list = ls())

61.2 Lab Steps

Step 1: Get the model into your workspace, then read the code (about five minutes)

The network model is unchanged from the tutorial, so there is nothing to re-specify and no reason to re-run the diagnostics. The netdx output is on the tutorial page if you want to look at it again.

Run mod13-COVID2.R from the top through the control.net() block, skipping the netdx() diagnostics block and stopping before the final netsim() call. The script sources mod13-COVID2-fx.R itself, so you do not need to source it separately. That leaves est, param, init, and control in your workspace, which is everything you need for the rest of the lab.

Now spend a few minutes reading rather than running. Two places in mod13-COVID2-fx.R are worth the time:

dx_covid(): find the two eligibility vectors, idsElig.sympt and idsElig.other, and see how dx.rate.sympt and dx.rate.other are applied to them. Note that dxStatus is written back with three levels (0 never screened, 1 screened negative, 2 screened positive), and that a false negative from imperfect pcr.sens leaves a truly infected person at dxStatus of 1.
infect2(): find the two lines del$status <- status[del$inf] and del$dxStatus <- dxStatus[del$inf]. These pull the disease state and the diagnosis state of the infected member of each discordant pair onto the edgelist. Everything the intervention does downstream depends on those two columns.

Discuss any conceptual questions about the module design with your group before you start running the model.

Step 2: Step through the infection module in browser mode

Uncomment the browser() line near the top of infect2() in mod13-COVID2-fx.R, re-source the file, and run netsim with a single simulation and a small number of steps so you land in the debugger quickly.

The seeds start as exposed rather than infectious, so for the first several steps there are no infectious nodes, discord_edgelist() is never reached, and del does not exist yet. Type c to continue until the inspection commands below actually return rows. Once they do, inspect the edgelist directly:

head(del)
table(del$status)
table(del$dxStatus)

Work out what a row of del represents, and confirm for yourself that del$status and del$dxStatus describe the infected node, not the susceptible one. Then walk forward through the lines that set del$transProb, del$actRate, and del$finalProb, and watch how the asymptomatic relative risk and the case isolation multiplier change individual rows. Case isolation is the single line

del$actRate[del$dxStatus == 2] <- del$actRate[del$dxStatus == 2] * act.rate.dx.inter.rr

which fires only when at >= act.rate.dx.inter.time. With the tutorial’s parameters that time is Inf, so the line never executes and the baseline model has no isolation at all.

Comment browser() back out before you move on.

Step 3: Run intervention scenarios

Now run a few versions of the model, changing one thing at a time in param.net() and comparing the epidemic output. At minimum, turn case isolation on by setting act.rate.dx.inter.time to a finite value, so the intervention starts partway through the simulation, and compare that run to the tutorial baseline.

Useful things to vary:

act.rate.dx.inter.time: when isolation starts.
act.rate.dx.inter.rr: how strong isolation is once it starts (the tutorial’s 0.05 is a 95% reduction in acts, which is close to perfect isolation; try something weaker and more plausible).
dx.rate.sympt and dx.rate.other: how fast symptomatic and non-symptomatic people are screened.
pcr.sens: test sensitivity, which controls how many true infections are missed and left at a dxStatus of 1.
allow.rescreen: whether someone already screened can be screened again. It is FALSE in the tutorial, which restricts non-symptomatic screening to people at dxStatus of 0, so screening is one-shot per person and raising dx.rate.other mostly exhausts the eligible pool faster. Setting it TRUE is what turns dx.rate.other into an ongoing screening program.
prop.clinical: the natural history side, which determines how many infections ever become symptomatic and therefore ever become eligible for the higher symptomatic screening rate.

Compare runs on cumulative incidence (sum(df$se.flow, na.rm = TRUE) on the mean data frame) as well as on the prevalence curves, and look at the diagnosis counters nDx, nDx.pos, and nDx.pos.fn to see what the screening program is doing mechanically.

The question to hold onto while you do this: how much does disease prevention depend on the natural history of disease, such as the proportion of infections in the subclinical pathway, versus the diagnostic assumptions, such as PCR sensitivity and the screening rates?

61.3 Optional Extensions

These are the take-home half of the lab. Each is specified in enough detail to implement on your own after the course.

Extension 1: Make diagnosis and case isolation age-dependent

Older people are both more likely to be screened in many real programs and more likely to be advised to isolate. Target either or both processes by age.

You already have a worked template for this pattern in the tutorial. The prop.clinical parameter is not a scalar; it is a vector of eight probabilities by age decade, c(0.40, 0.25, 0.37, 0.42, 0.51, 0.59, 0.72, 0.76), and progress2() indexes it with

age.group <- pmin((floor(age[ids.newInf] / 10)) + 1, 8)
prop.clin.vec <- prop.clinical[age.group]

Read those lines carefully before you write anything. They are the same “turn a scalar parameter into a per-node probability vector” move you need here.

The diagnosis half, in dx_covid(). Currently dx.rate.sympt and dx.rate.other are scalars fed straight to rbinom(nElig, 1, rate). To make them age-dependent, pass each as an age-decade vector in param.net(), pull age off dat with get_attr(), and build a per-node probability vector before the rbinom call. The pitfall is positional alignment: the probability vector must be in the same order and of the same length as the eligible ids it is drawn against, so index by age[idsElig.sympt] for the symptomatic draw and by age[idsElig.other] for the other draw. A vector built off the full population, or off the wrong eligibility set, will recycle silently and assign the wrong person’s screening rate without throwing an error. The stop() guard that progress2() uses on prop.clin.vec is a cheap way to catch that.

The isolation half is a different pattern, because isolation is not applied to nodes; it is applied to rows of the discordant edgelist inside infect2(). There is no age column on del, so you have to put one there: del$age <- age[del$inf] after pulling age with get_attr(), in the same style as the existing del$status and del$dxStatus lines. Then make act.rate.dx.inter.rr an age-decade vector and index it by the infected node’s age when you apply the multiplier. Decide deliberately whose age should govern the isolation of a partnership, the infected node’s or the susceptible node’s, and be able to defend the choice.

Extension 2: Contact tracing

Contact tracing means identifying the recent contacts of a diagnosed case and then isolating, screening, or otherwise intervening on them. It is the natural next intervention after diagnosis and case isolation, and unlike everything else in this model it requires the simulation to remember partnerships that have already dissolved.

Rather than deriving it from scratch, start from the maintained EpiModel Gallery example SEIR with Contact Tracing for an Acute, Immunizing Infection. It already solves most of the problem:

a COVID-like SEIR with a presymptomatic infectious split (ip then is), so transmission happens before symptoms and there is something for tracing to catch;
a trace module that finds newly diagnosed indices, calls get_partners() on the cumulative edgelist with a lookback window, and translates the returned unique ids back to positional ids with get_posit_ids();
a Bernoulli reach probability, so only some traced partners are actually reached;
quarantine implemented as an act rate multiplier inside the infection module, the same mechanism as the case isolation you studied in Step 2;
four comparison scenarios that separate tracing speed from tracing coverage.

Clone the whole Gallery repository and run from the repository root, not from inside the example folder:

# from the EpiModel-Gallery repo root
source("examples/seir-contact-tracing/model.R")

Running from the root is required, not a preference: model.R loads its module functions with source("examples/seir-contact-tracing/module-fx.R"), a path written relative to the repository root.

One prerequisite before you port any of this into the course model

The tracing pattern depends on the cumulative edgelist, the running history of every partnership that has been active in the simulation. Our control.net() call does not set cumulative.edgelist, and it defaults to FALSE. If you copy a trace module into this model without setting cumulative.edgelist = TRUE, get_partners() will return nothing, no error will be raised, and the intervention will appear to do nothing at all. Set cumulative.edgelist = TRUE, and set truncate.el.cuml to the number of steps you intend to trace back over, which bounds how much history is held in memory. Do not leave truncate.el.cuml at its default of 0, which retains only currently active partnerships: that is a second silent failure, since tracing would then only ever find partners the index still has.

Once tracing runs, the interesting question is the same one as Step 3, asked of a different intervention: does tracing speed or tracing coverage matter more, and why does the answer depend on the length of the presymptomatic window?

61.4 Lab Questions

Questions

In Step 2, what does one row of the discordant edgelist represent, and why do del$status and del$dxStatus describe the infected node rather than the susceptible one? What would break if they described the susceptible node instead?
Case isolation is implemented as a reduction in act.rate on the rows where del$dxStatus == 2, not as a removal of the partnership. What behavior does that represent, what does it fail to represent, and what would you have to change to model a diagnosed person dropping partnerships entirely?
From your Step 3 scenarios: did prevention depend more on the natural history assumptions (the proportion of infections in the subclinical pathway) or on the diagnostic assumptions (PCR sensitivity and the screening rates)? Explain the result in terms of who is eligible for symptomatic screening in the first place, and what that implies for a screening program facing a pathogen with substantial asymptomatic transmission.

61.5 Next Steps

This is the last lab of the course, so it is a good moment to point the model at your own questions. What intervention scenarios have we not built out for COVID that are of interest to you? Vaccination, targeted screening of high-degree individuals, and time-varying behavioral change are all reachable from where this model already is. What are the natural next steps for building one of those out in code, and which module would each one touch? Sketch it out and we will help.

61.6 Worked solution

A worked solution is in the course repository as mod13-COVID2-Lab-Solution.R. It runs the core steps, including the intervention scenario comparison in Step 3, and implements the age-dependent diagnosis and isolation extension with the alignment issue handled explicitly. Please work through the lab yourself before you consult it. You will get much more out of the questions if you have already formed your own answers.

# COVID Intervention Lab ```{r setup, include = FALSE} library(knitr) opts_chunk$set(comment = NA) par(mar = c(3, 3, 1, 1), mgp = c(2, 1, 0)) ``` In this lab you will work with the diagnosis and case isolation features of the COVID model from the previous tutorial, and use them to ask how much transmission a screening-based intervention can actually prevent. The tutorial built the machinery; this lab turns it into experiments. The specific learning objectives are to: 1. Read the diagnosis module and the discordant edgelist code closely enough to know exactly where an intervention acts on transmission; 2. Run intervention scenarios that vary the diagnosis rates and the case isolation intensity, and interpret the results against the natural history assumptions; 3. Define concrete next steps for extending diagnosis-based interventions, including contact tracing. ::: {.callout-note} ## How to use your time Steps 1 through 3 are the core of the lab and are what everyone should finish in the session. The two extensions that follow are optional. They are specified in enough detail that you can start one here and finish it at home after the course, which is how we expect most people will use them. ::: ## Setup Once you are ready, start out by clearing your R object environment, to make sure that you do not have any objects lingering from the tutorial. This can be accomplished with: ```{r} rm(list = ls()) ``` ## Lab Steps ::: {.callout-note .instructions} ## Step 1: Get the model into your workspace, then read the code (about five minutes) The network model is unchanged from the tutorial, so there is nothing to re-specify and no reason to re-run the diagnostics. The `netdx` output is on the [tutorial page](mod13-COVID2.html) if you want to look at it again. Run [`mod13-COVID2.R`](mod13-COVID2.R) from the top through the `control.net()` block, skipping the `netdx()` diagnostics block and stopping before the final `netsim()` call. The script sources [`mod13-COVID2-fx.R`](mod13-COVID2-fx.R) itself, so you do not need to source it separately. That leaves `est`, `param`, `init`, and `control` in your workspace, which is everything you need for the rest of the lab. Now spend a few minutes reading rather than running. Two places in [`mod13-COVID2-fx.R`](mod13-COVID2-fx.R) are worth the time: - `dx_covid()`: find the two eligibility vectors, `idsElig.sympt` and `idsElig.other`, and see how `dx.rate.sympt` and `dx.rate.other` are applied to them. Note that `dxStatus` is written back with three levels (0 never screened, 1 screened negative, 2 screened positive), and that a false negative from imperfect `pcr.sens` leaves a truly infected person at `dxStatus` of 1. - `infect2()`: find the two lines `del$status <- status[del$inf]` and `del$dxStatus <- dxStatus[del$inf]`. These pull the disease state and the diagnosis state of the *infected* member of each discordant pair onto the edgelist. Everything the intervention does downstream depends on those two columns. Discuss any conceptual questions about the module design with your group before you start running the model. ::: ::: {.callout-note .instructions} ## Step 2: Step through the infection module in browser mode Uncomment the `browser()` line near the top of `infect2()` in [`mod13-COVID2-fx.R`](mod13-COVID2-fx.R), re-source the file, and run `netsim` with a single simulation and a small number of steps so you land in the debugger quickly. The seeds start as exposed rather than infectious, so for the first several steps there are no infectious nodes, `discord_edgelist()` is never reached, and `del` does not exist yet. Type `c` to continue until the inspection commands below actually return rows. Once they do, inspect the edgelist directly: ```r head(del) table(del$status) table(del$dxStatus) ``` Work out what a row of `del` represents, and confirm for yourself that `del$status` and `del$dxStatus` describe the infected node, not the susceptible one. Then walk forward through the lines that set `del$transProb`, `del$actRate`, and `del$finalProb`, and watch how the asymptomatic relative risk and the case isolation multiplier change individual rows. Case isolation is the single line ```r del$actRate[del$dxStatus == 2] <- del$actRate[del$dxStatus == 2] * act.rate.dx.inter.rr ``` which fires only when `at >= act.rate.dx.inter.time`. With the tutorial's parameters that time is `Inf`, so the line never executes and the baseline model has no isolation at all. Comment `browser()` back out before you move on. ::: ::: {.callout-note .instructions} ## Step 3: Run intervention scenarios Now run a few versions of the model, changing one thing at a time in `param.net()` and comparing the epidemic output. At minimum, turn case isolation on by setting `act.rate.dx.inter.time` to a finite value, so the intervention starts partway through the simulation, and compare that run to the tutorial baseline. Useful things to vary: - `act.rate.dx.inter.time`: when isolation starts. - `act.rate.dx.inter.rr`: how strong isolation is once it starts (the tutorial's 0.05 is a 95% reduction in acts, which is close to perfect isolation; try something weaker and more plausible). - `dx.rate.sympt` and `dx.rate.other`: how fast symptomatic and non-symptomatic people are screened. - `pcr.sens`: test sensitivity, which controls how many true infections are missed and left at a `dxStatus` of 1. - `allow.rescreen`: whether someone already screened can be screened again. It is `FALSE` in the tutorial, which restricts non-symptomatic screening to people at `dxStatus` of 0, so screening is one-shot per person and raising `dx.rate.other` mostly exhausts the eligible pool faster. Setting it `TRUE` is what turns `dx.rate.other` into an ongoing screening program. - `prop.clinical`: the natural history side, which determines how many infections ever become symptomatic and therefore ever become eligible for the higher symptomatic screening rate. Compare runs on cumulative incidence (`sum(df$se.flow, na.rm = TRUE)` on the mean data frame) as well as on the prevalence curves, and look at the diagnosis counters `nDx`, `nDx.pos`, and `nDx.pos.fn` to see what the screening program is doing mechanically. The question to hold onto while you do this: how much does disease prevention depend on the natural history of disease, such as the proportion of infections in the subclinical pathway, versus the diagnostic assumptions, such as PCR sensitivity and the screening rates? ::: ## Optional Extensions These are the take-home half of the lab. Each is specified in enough detail to implement on your own after the course. ::: {.callout-note .instructions} ## Extension 1: Make diagnosis and case isolation age-dependent Older people are both more likely to be screened in many real programs and more likely to be advised to isolate. Target either or both processes by age. You already have a worked template for this pattern in the tutorial. The `prop.clinical` parameter is not a scalar; it is a vector of eight probabilities by age decade, `c(0.40, 0.25, 0.37, 0.42, 0.51, 0.59, 0.72, 0.76)`, and `progress2()` indexes it with ```r age.group <- pmin((floor(age[ids.newInf] / 10)) + 1, 8) prop.clin.vec <- prop.clinical[age.group] ``` Read those lines carefully before you write anything. They are the same "turn a scalar parameter into a per-node probability vector" move you need here. **The diagnosis half**, in `dx_covid()`. Currently `dx.rate.sympt` and `dx.rate.other` are scalars fed straight to `rbinom(nElig, 1, rate)`. To make them age-dependent, pass each as an age-decade vector in `param.net()`, pull `age` off `dat` with `get_attr()`, and build a per-node probability vector before the `rbinom` call. The pitfall is positional alignment: the probability vector must be in the same order and of the same length as the eligible ids it is drawn against, so index by `age[idsElig.sympt]` for the symptomatic draw and by `age[idsElig.other]` for the other draw. A vector built off the full population, or off the wrong eligibility set, will recycle silently and assign the wrong person's screening rate without throwing an error. The `stop()` guard that `progress2()` uses on `prop.clin.vec` is a cheap way to catch that. **The isolation half** is a different pattern, because isolation is not applied to nodes; it is applied to rows of the discordant edgelist inside `infect2()`. There is no `age` column on `del`, so you have to put one there: `del$age <- age[del$inf]` after pulling `age` with `get_attr()`, in the same style as the existing `del$status` and `del$dxStatus` lines. Then make `act.rate.dx.inter.rr` an age-decade vector and index it by the infected node's age when you apply the multiplier. Decide deliberately whose age should govern the isolation of a partnership, the infected node's or the susceptible node's, and be able to defend the choice. ::: ::: {.callout-note .instructions} ## Extension 2: Contact tracing Contact tracing means identifying the recent contacts of a diagnosed case and then isolating, screening, or otherwise intervening on them. It is the natural next intervention after diagnosis and case isolation, and unlike everything else in this model it requires the simulation to remember partnerships that have already dissolved. Rather than deriving it from scratch, start from the maintained EpiModel Gallery example [**SEIR with Contact Tracing for an Acute, Immunizing Infection**](https://epimodel.github.io/EpiModel-Gallery/examples/seir-contact-tracing/){target="_blank" rel="noopener"}. It already solves most of the problem: - a COVID-like SEIR with a presymptomatic infectious split (`ip` then `is`), so transmission happens before symptoms and there is something for tracing to catch; - a `trace` module that finds newly diagnosed indices, calls `get_partners()` on the cumulative edgelist with a lookback window, and translates the returned unique ids back to positional ids with `get_posit_ids()`; - a Bernoulli reach probability, so only some traced partners are actually reached; - quarantine implemented as an act rate multiplier inside the infection module, the same mechanism as the case isolation you studied in Step 2; - four comparison scenarios that separate tracing speed from tracing coverage. **Clone the whole [Gallery repository](https://github.com/EpiModel/EpiModel-Gallery){target="_blank" rel="noopener"} and run from the repository root**, not from inside the example folder: ```r # from the EpiModel-Gallery repo root source("examples/seir-contact-tracing/model.R") ``` Running from the root is required, not a preference: `model.R` loads its module functions with `source("examples/seir-contact-tracing/module-fx.R")`, a path written relative to the repository root. ::: {.callout-note .concept} ## One prerequisite before you port any of this into the course model The tracing pattern depends on the **cumulative edgelist**, the running history of every partnership that has been active in the simulation. Our `control.net()` call does not set `cumulative.edgelist`, and it defaults to `FALSE`. If you copy a trace module into this model without setting `cumulative.edgelist = TRUE`, `get_partners()` will return nothing, no error will be raised, and the intervention will appear to do nothing at all. Set `cumulative.edgelist = TRUE`, and set `truncate.el.cuml` to the number of steps you intend to trace back over, which bounds how much history is held in memory. Do not leave `truncate.el.cuml` at its default of 0, which retains only currently active partnerships: that is a second silent failure, since tracing would then only ever find partners the index still has. ::: Once tracing runs, the interesting question is the same one as Step 3, asked of a different intervention: does tracing speed or tracing coverage matter more, and why does the answer depend on the length of the presymptomatic window? ::: ## Lab Questions ::: {.callout-note .question} ## Questions 1. In Step 2, what does one row of the discordant edgelist represent, and why do `del$status` and `del$dxStatus` describe the infected node rather than the susceptible one? What would break if they described the susceptible node instead? 2. Case isolation is implemented as a reduction in `act.rate` on the rows where `del$dxStatus == 2`, not as a removal of the partnership. What behavior does that represent, what does it fail to represent, and what would you have to change to model a diagnosed person dropping partnerships entirely? 3. From your Step 3 scenarios: did prevention depend more on the natural history assumptions (the proportion of infections in the subclinical pathway) or on the diagnostic assumptions (PCR sensitivity and the screening rates)? Explain the result in terms of who is eligible for symptomatic screening in the first place, and what that implies for a screening program facing a pathogen with substantial asymptomatic transmission. ::: ## Next Steps This is the last lab of the course, so it is a good moment to point the model at your own questions. What intervention scenarios have we not built out for COVID that are of interest to you? Vaccination, targeted screening of high-degree individuals, and time-varying behavioral change are all reachable from where this model already is. What are the natural next steps for building one of those out in code, and which module would each one touch? Sketch it out and we will help. ## Worked solution A worked solution is in the course repository as [`mod13-COVID2-Lab-Solution.R`](mod13-COVID2-Lab-Solution.R). It runs the core steps, including the intervention scenario comparison in Step 3, and implements the age-dependent diagnosis and isolation extension with the alignment issue handled explicitly. Please work through the lab yourself before you consult it. You will get much more out of the questions if you have already formed your own answers.