-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlisteriaSummary.Rmd
94 lines (75 loc) · 5.28 KB
/
listeriaSummary.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
---
title: "Modelling the decline in Listeria cases"
subtitle: "Initial concept note - Confidential"
date: "Last updated: 2018-03-29"
author: "Prepared by SACEMA"
output:
pdf_document:
fig_caption: yes
highlight: tango
latex_engine: xelatex
html_document: default
mainfont: Arial
geometry: margin=.75in
fontsize: 12pt
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(tidyverse)
require(bbmle)
source('fitFunctions.R')
source('epiCurve.R')
sampleDelay <- 3
bsIterations <- 200
source('incubationDistSens.R') # Assumes samples are acquired within sampleDelay days of onset
#source('incubationDist.R')
source('epiFunctions.R')
baseline <- mean(cases$count)
recallEffect <- 0.65 # after 1 week
expectedShelfLife <- 2 # weeks
```
## Baseline cases
The expected baseline number of cases was set to the mean weekly number of cases for Week 40 of 2017 through week 8 of 2018 (n = `r round(baseline,1)`). Observed variation around this mean was approximately Poisson distributed, so Poisson noise around the expected value is assumed in the model.
## Decline in risky exposures
An initial assessment suggests that the product recall was `r recallEffect*100`\% effective within one week. We assume a corresponding decline in risky exposures over this timeframe, but this assumption is easily modifiable should better information become available. We assume a further geometric decline through time as the remaining contaminated product is consumed or discarded. We have initially assumed the average remaining time that the contaminated product poses a risk is `r round(expectedShelfLife,1)` weeks. This time may actually be shorter, and the assumption is easily modified. The resulting temporal profile for risky exposures is represetnted below (Figure 1A).
```{r exposure, echo=FALSE}
ep <- exposureProfile(baseline,recallEffect)
minGen <- qpois(0.025,ep$riskExp)
maxGen <- qpois(0.975,ep$riskExp)
```
## Anticipated future cases
The number of future cases is determined by two primary factors: the temporal profile of risky exposures and the incubation period. The incubation period distribution has been characterized, initially, based on data from previous outbreaks in Europe and the US (see Goulet _et al_. 2013 _BMC Infectious Diseases_, Figure 1B). As data on the incubation period becomes available for the current outbreak in South Africa, we will use these data to update estimates that are likely to be more relevant for this outbreak, given the differences in background population and the source of exposure. Note that we cannot confidently capture the tail of the distribution based on the data available.
```{r,echo=FALSE,fig.width=7.5, fig.height=3, fig.cap='(A) Temporal profile of risky expsoures, shown as plausible range for future cases generated each week, based on the stated model assumptions. (B) Comparison of non-parametric (grey) and parametric (green) survival analyses for incubation data. Based on data from Goulet 2013.'}
par(mar=c(5,4,2,2), mfcol=c(1,2))
plot(ep$week+.5,rep(NA,nrow(ep)),ylim=c(0,45),bty='L',ann=F,yaxt='n',xaxt='n',yaxs='i',xpd=T,type='s',lwd=3)
arrows(ep$week+.5,minGen,ep$week+.5,maxGen,le=.02,angle = 90,lwd=3,code = 3,xpd=T,col='darkgrey')
axis(2,seq(0,45,5),cex.axis=.85)
axis(1,seq(1,18,2),cex.axis=.85)
mtext('Week of the year',1,3)
mtext('Future cases generated',2,2.5)
mtext('A',3,line=1,adj=-.3,font=2)
par(mar=c(5,4,2,1))
plot(npmle_perInc,ylim=c(0,1),bty='L',ann=F,yaxt='n',yaxs='i',xpd=NULL)
axis(2,seq(0,1,.4),cex.axis=.85)
axis(2,seq(.2,1,.4),cex.axis=.85)
mtext('Time since infection (days)',1,3)
mtext('Proportion incubating',2,2.5)
mtext('B',3,line=1,adj=-.3,font=2)
curve(1-pgamma(x,meanIncGamma*rateIncGamma,rateIncGamma),add=T,lwd=3,col='darkgreen')
points(meanIncGamma,1.1,pch=21,cex=1.1,xpd=T,col='darkgreen',lwd=3)
arrows(ciIncGamma[1],1.1,ciIncGamma[2],1.1,le=.02,angle = 90,lwd=3,code = 3,col='darkgreen',xpd=T)
mtext('Mean and 95% CI',3,line=.5,adj=1.1,font=2,col='darkgreen',cex=.8,xpd=NULL)
abline(v=3.5+seq(0,70,7),col='red',lty=3,xpd=F)
```
Combining the temporal profile for risky exposures and the estimated distribution of the incubation period, we can begin to assess the anticipated decline in the epidemic curve as a result of increased awareness and product recall. Figure 2 shows plausible ranges for the number of cases under the assumptions of the model.
```{r,echo=FALSE,fig.width=4, fig.height=4, fig.cap='Range of cases expected through time, from Week 8 to Week 18 of 2018.'}
weekRange <- 8:18
par(mar=c(5,4,2,2))
plot(weekRange,rep(NA,11),ylim=c(0,45),bty='L',yaxs='i',xlab='Week of the year',ylab='Anticipated number of cases')
arrows(weekRange,projCurve(minGen)[weekRange],weekRange,projCurve(maxGen)[weekRange],le=.02,angle = 90,lwd=3,code = 3,col='darkblue',xpd=T)
#projCurve(ep$riskExp)
# Can we get to probablility of seeing no cases for two consecutive weeks through time?
```
## Interpretation
**If the number of cases through time is consistently or dramatically above the maximum of the expected range, this indicates that the assumptions of the model are violated in important ways and is probably a sign that there is an additonal source of risky exposure that remains unidentified.**
**Contact:** Please address questions, feedback, and requests related to this document to Prof. Juliet Pulliam, Director of SACEMA, at <[email protected]>.