Ratio vs subtraction method for seropositivity

serology_roc.R

# Define parameters
bkg_val <- 15 # background value
xx <- 1e3 # Number of values to simulate from positive and negative distribution
xx_pos <- xx; xx_neg <- xx

set.seed(1)

# Set up plot
par(mfcol=c(2,2),mgp=c(2,0.7,0),mar = c(3,3,1,3))

for(kk in 1:2){
  
  # Define positive/negative distribution scenarios
  if(kk==1){
    neg_dist <- rnorm(xx_neg,mean=15,sd=5) %>% pmax(.,0.1)
    pos_dist <- rnorm(xx_pos,mean=40,sd=15) %>% pmax(.,0.1) %>% pmin(.,200) #rlnorm(xx,mean=log(40),sd=log(1.5)) %>% pmax(.,0.1) %>% pmin(.,200)
  }else{
    neg_dist <- rnorm(xx_neg,mean=30,sd=10) %>% pmax(.,0.1)
    pos_dist <- rnorm(xx_pos,mean=60,sd=10) %>% pmax(.,0.1) %>% pmin(.,200)
  }
  
  # Plot positive negative distributions
  break_n <- seq(0,200,4)
  hist(pos_dist,xlim=c(0,100),col=rgb(0,0.5,0,0.4),ylim=c(0,0.1),prob=T,breaks=break_n,main="",xlab="assay value")
  hist(neg_dist,add=T,col=rgb(1,0,0,0.4),prob=T,breaks=break_n)
  lines(c(bkg_val,bkg_val),c(0,1e3),lty=2,lwd=2)
  graphics::text(x=bkg_val+1,y=0.1,labels="background value",adj=0)
  
  # Calculate false positive/negative for different cutoffs
  inc1 <- 0.01; xx1 <- seq(0,5,inc1) # Range of cutoffs for ratio
  inc2 <- 0.2; xx2 <- seq(-50,50,inc2) # Range of cutoffs for subtraction
  
  # Proportion of positives above cutoff
  prop_cut1 <- sapply(xx1,function(x){sum((pos_dist/bkg_val)>x)/xx})
  prop_cut2 <- sapply(xx2,function(x){sum((pos_dist-bkg_val)>x)/xx})
  
  # Proportion of negatives below cutoff
  n_cut1 <- sapply(xx1,function(x){sum((neg_dist/bkg_val)<x)/xx})
  n_cut2 <- sapply(xx2,function(x){sum((neg_dist-bkg_val)<x)/xx})
  
  # Plot ROC curves
  plot(1-n_cut1,prop_cut1,xlim=c(0,1),ylim=c(0,1),xlab="false positive rate",ylab="true positive rate",type="l",lwd=2)
  lines(1-n_cut2,prop_cut2,col="blue",lwd=2)
  graphics::text(x=0.5,0.8,labels="ratio",col="black",adj=0)
  graphics::text(x=0.5,0.7,labels="subtraction",col="blue",adj=0)

}

dev.copy(pdf,"MIA_simulation.pdf",width=6,height=6)
dev.off()