Skip to content

Instantly share code, notes, and snippets.

@adamkucharski
Created June 19, 2024 20:26
Show Gist options
  • Save adamkucharski/1542a39fec9bab7dfe80800f25069956 to your computer and use it in GitHub Desktop.
Save adamkucharski/1542a39fec9bab7dfe80800f25069956 to your computer and use it in GitHub Desktop.
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment