serology_roc_vary_background.R

# Define parameters
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
  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)

  # Define background distribution
  if(kk==1){
    bkg_val <- rnorm(xx_neg,mean=15,sd=5) %>% pmax(.,0.1) # background value
  }else{
    bkg_val <- rnorm(xx_neg,mean=50,sd=10) %>% pmax(.,0.1) # background value
    
  }
  
  # 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)
  hist(bkg_val,add=T,col=rgb(0,0,0,0.5),prob=T,breaks=break_n)
  graphics::text(x=median(bkg_val)+3,y=0.1,labels="background values",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()