Title:  Accounting for incomplete testing in estimation of epidemic parameters
Authors:  Rebecca A. Betensky and Yang Feng

R functions 


#Figure 1
doubtime=function(){
postscript(file="doubtimesens9.ps")
par(mfrow=c(3,4),omi=c(.5,.5,.5,.5))
daily=read.csv(file="daily.csv",header=T)
#state=unique(daily$state)
state=c( "CA","FL","GA","IL","LA","MA","MI","NJ","NY","PA","TX","WA")
lstate=length(state)
for(i in 1:lstate){
dat=daily[daily$state==state[i],]
pos=dat$positive
ldays=length(pos)
compare=(pos[1]/pos[(2):10])
ind=1:9
doub0= max(ind[compare[ind]<2])
slp= min(compare[compare>=2])-max(compare[compare <2])
inter=max(compare[compare<2])-slp*doub0
doub00=(2-inter)/slp
ll=0
for(eps in c(.11,.25,.33)){
ll=ll+1
doub=rep(0,9)
doubb=doub
for(j in 1:9){
rat=rep(1,9)
rat[j:9]=1+eps
# j indicates the day at which probability of testing given pos decreased relative to previous day by factor of 1/(1+eps) (i.e., 90%, 80%, 75%): assume constant across all other days
ind=1:9
doub[j]= max(ind[compare[ind]*rat[ind]<2])
cr=compare[ind]*rat[ind]
doub[j]= max(ind[cr[ind]<2])
slp= min(cr [cr>=2])-max(cr[cr <2])
inter=max(cr[cr<2])-slp*doub[j]
doubb[j]=(2-inter)/slp
}
if(ll==1) {plot(0:(-8),doubb[1:9],type="l",xlab="day of decrease in testing",ylab="doubling time",ylim=c(0,8),xlim=c(0,-8),cex=.75)
title(paste(state[i]))
points(0,doub00)}
#points(1,doub0,pch=2)}
if(ll>1) lines(0:(-8),doubb[1:9],lty=ll,col=ll)
}}}

#Figure 2
doubtime2=function(){
postscript(file="doubtimesens29.ps")
par(mfrow=c(3,4),omi=c(.5,.5,.5,.5))
daily=read.csv(file="daily.csv",header=T)
state=unique(daily$state)
state=c( "CA","FL","GA","IL","LA","MA","MI","NJ","NY","PA","TX","WA")
lstate=length(state)
for(i in 1:lstate){
dat=daily[daily$state==state[i],]
pos=dat$positive
pos[is.na(pos)==T]=0
ldays=length(pos)
doub00=rep(0,20)
for(j in 1:20){
maxind=min(j+9,length(pos))
chk=pos[j:maxind]
if(length(chk[chk>0])>2){
compare=(pos[j]/pos[(j+1):maxind])
compare=compare[compare<Inf]
if(min(compare)<2 & max(compare)>2){
ind=1:length(compare)
doub0= max(ind[compare[ind]<2])
slp= min(compare[compare>=2])-max(compare[compare <2])
inter=max(compare[compare<2])-slp*doub0
doub00[j]=(2-inter)/slp
}
if(min(compare)>2 | max(compare)<2) doub00[j]=1
}
}
plot(0:(-19),doub00[1:20],type="l",xlab="day",ylab="doubling time",ylim=c(0,8),xlim=c(0,-19),cex=.75)
title(paste(state[i]))
ll=1
for(eps in c(.11,.25,.33)){
ll=ll+1
doubb=rep(0,20)
for(j in 1:20){
maxind=min(j+9,length(pos))
chk=pos[j:maxind]
if(length(chk[chk>0])>2){
compare=(pos[j]/pos[(j+1):maxind])
ind=1:length(compare)
cr=compare[ind]*(1+eps)
if(min(cr)<2 & max(cr)>2){
doub= max(ind[cr[ind]<2])
slp= min(cr [cr>=2])-max(cr[cr <2])
inter=max(cr[cr<2])-slp*doub
doubb[j]=(2-inter)/slp
}
if(min(cr)>=2 | max(cr)<=2) doubb[j]=1
}}
lines(0:(-19),doubb[1:20],lty=ll,col=ll)
}}}



#Figure 3
curve=function(){
postscript(file="curve9.ps")
par(mfrow=c(3,4),omi=c(.5,.5,.5,.5))
daily=read.csv(file="daily.csv",header=T)
#state=unique(daily$state)
state=c( "CA","FL","GA","IL","LA","MA","MI","NJ","NY","PA","TX","WA")
lstate=length(state)
for(i in 1:lstate){
dat=daily[daily$state==state[i],]
pos=dat$positive
if(pos[10]>0){
ldays=length(pos)
compare=(pos[1:9]/pos[10])
plot(-8:0,(compare[9:1]),type="l", xlab="day",ylab="cases relative to day -8", cex=.4,ylim=c(0,7),xlim=c(-8,0),lwd=2)
title(paste(state[i]))
ll=1
for(eps in c(.11,.25,.33)){
ll=ll+1
lines(-8:0,(compare[9:1]*(1+eps)),lty=ll,col=ll)
}}}}