Illustrative plots for PCR positivity

2024_01_delay_PCR.R

# - - - - - - - - - - - - - - - - - - - - - - - 
# Code by Adam Kucharski to accompany post 'Counting current COVID infections' (2nd Jan 2024)
# Post is available at: https://kucharski.substack.com
# Code in this file is shared under an MIT license (https://opensource.org/license/mit/)
# - - - - - - - - - - - - - - - - - - - - - - - 

# Load libraries
library(dplyr)
library(readr)

# Import PCR positivity estimates from Hellewell et al (BMC Med, 2021)
pcr_curve <- read_csv("https://raw.githubusercontent.com/cmmid/pcr-profile/main/PCR_curve_summary.csv")

# Convert median PCR curves to daily probabilities
pcr_curve <- pcr_curve |>  mutate(p_neg = 1-median)

# Calculate the daily median value, averaging across each day to get a discretised distribution
daily_data <- pcr_curve |> 
  mutate(day = round(days_since_infection)) |> 
  group_by(day) |> 
  summarise(daily_median = mean(median)) # Or use sum(median) if appropriate

# Function to extract daily value
p_by_day <- function(x){daily_data[match(x,daily_data$day),]$daily_median}


# Plot PCR positivity over time
max_days <- 30 # Time period to plot

par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(0:max_days,p_by_day(0:max_days),type="l",xlab="days after infection",ylab="probability will test positive by PCR",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(0,0.8))
points(0:max_days,p_by_day(0:max_days),col="dark orange",pch=19)

dev.copy(png,paste0("delay.png"),units="cm",width=15,height=10,res=150)
dev.off()

# Plot normalised curves, i.e. P(infected on a given day | positive)
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(0:max_days,p_by_day(0:max_days)/sum(p_by_day(0:max_days)),type="l",xlab="days since infection",ylab="probability",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(0,0.12))
points(0:max_days,p_by_day(0:max_days)/sum(p_by_day(0:max_days)),col="dark orange",pch=19)

graphics::arrows(x0=25,y0=0.085,x1=25,y1=0.01,length=0.1,lwd=2)
graphics::text(x=25,y=0.09,labels="less likely",adj=0.5)

graphics::arrows(x0=0.5,y0=0.105,x1=0.5,y1=0.01,length=0.1,lwd=2)
graphics::text(x=1,y=0.11,labels="less likely",adj=0.2)

graphics::arrows(x0=8,y0=0.105,x1=6,y1=0.09,length=0.1,lwd=2)
graphics::text(x=8,y=0.11,labels="more likely",adj=0)

dev.copy(png,paste0("delay_norm.png"),units="cm",width=15,height=10,res=150)
dev.off()


# Infection timing given exponential growth ------------------------------------------------------

# Define time period and growth
xx_days <- 0:max_days
daily_growth <- 0.06
yy <- 100*exp(xx_days*daily_growth)

# Plot growth in new infections over time
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(xx_days,yy,type="l",xlab="days",ylab="new infections (scaled)",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(100,610))
points(xx_days,yy,col="dark orange",pch=19)

dev.copy(png,paste0("infections.png"),units="cm",width=15,height=10,res=150)
dev.off()

# Plot P(infection time | infected) in growing epidemic
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(xx_days,rev(yy)/sum(yy),type="l",xlab="days since infection",ylab="probability",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(0,0.1))
points(xx_days,rev(yy)/sum(yy),col="dark orange",pch=19)
       
dev.copy(png,paste0("infections_time.png"),units="cm",width=15,height=10,res=150)
dev.off()

# Adjust P(infection time | positive) for epidemic dynamics
adjust_yy <- (p_by_day(0:max_days)/sum(p_by_day(0:max_days)))*rev(yy)/sum(yy)

# Plot P(infection time | positive) in growing and flat epidemic
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(xx_days,adjust_yy/sum(adjust_yy),type="l",xlab="days since infection",ylab="probability",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(0,0.14))
lines(0:max_days,p_by_day(0:max_days)/sum(p_by_day(0:max_days)),lwd=2,col="dark orange",lty=2)
points(xx_days,adjust_yy/sum(adjust_yy),col="dark orange",pch=19)

graphics::arrows(x0=8,y0=0.125,x1=6,y1=0.11,length=0.1,lwd=2)
graphics::text(x=7,y=0.13,labels="during growing epidemic",adj=0)

dev.copy(png,paste0("PCR_delay_adjusted.png"),units="cm",width=15,height=12,res=150)
dev.off()


# New daily infections vs PCR positivity ------------------------------------------------------

# Simulate 3 scenarios
xx_days <- 0:max_days
yy <- 0.16*exp(xx_days*0.06)
yy2 <- 1.67*exp(-xx_days*0.0527)
yy3 <- 0.43*sin(40*pi*xx_days/365)+0.45

# Plot new infections over time
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(xx_days,yy,type="l",xlab="days",ylab="new infections (% population)",
     xaxs="i",yaxs="i",bty="l",lwd=2,col="dark orange",ylim=c(0,2.5),xlim=c(0,30.5))
points(xx_days,yy,col="dark orange",pch=19)

lines(xx_days,yy2,col="dark blue")
points(xx_days,yy2,col="dark blue",pch=19)

lines(xx_days,yy3,col="dark green")
points(xx_days,yy3,col="dark green",pch=19)

# Calculate % PCR positive on day = max_days
print(sum(rev(yy)*p_by_day(0:max_days)))
print(sum(rev(yy2)*p_by_day(0:max_days)))
print(sum(rev(yy3)*p_by_day(0:max_days)))

dev.copy(png,paste0("infections_to_prevalence.png"),units="cm",width=15,height=10,res=150)
dev.off()