title: Simple way to plot mean/average ROC curve in R
body: <p>I am doing molecular docking research by testing 4 different docking techniques. I have 4 data tables: each has columns of Ligand Name, Binding Score, and Whether the compound is active (active = 1, not active = 0). I have plotted ROC curves for each technique on a graph, but I cannot figure out how to plot a "mean", or average ROC curve for these. </p>
  
  <p>Here is a picture of my ROC Curve
  <img src="https://i.stack.imgur.com/gZXpz.jpg" alt="plot">.</p>
  
  <p>I do not know if this is possible with the data charts I have, or if I need to create a new file with ROC prediction data. If so, how can I do this?</p>
  
  <p>I am new to R, so a simpler explanation would be greatly appreciated!</p>
  
  <pre><code>#load ROCR
  library(ROCR)
  library(pROC)
  
  #load ligands and decoys
  lig &lt;- unique(read.table("C:/ROC-Curve/Input-Files/actives")[,1]);
  dec &lt;- unique(read.table("C:/ROC-Curve/Input-Files/decoys")[,1]);
  
  #load data file from docking
  m10n &lt;- read.table("C:/ROC-Curve/Input-Files/med10-nmr",header=T);
  m10x &lt;- read.table("C:/ROC-Curve/Input-Files/med10-xray",header=T);
  m30n &lt;- read.table("C:/ROC-Curve/Input-Files/med30-nmr-results",header=T);
  m30x &lt;- read.table("C:/ROC-Curve/Input-Files/med30-xray-results",header=T);
  
  #change colnames
  colnames(m10n)[1]="LigandName";
  colnames(m10n)[2]="BindingScore"
  
  colnames(m10x)[1]="LigandName";
  colnames(m10x)[2]="BindingScore"
  
  colnames(m30n)[1]="LigandName";
  colnames(m30n)[2]="BindingScore"
  
  colnames(m30x)[1]="LigandName";
  colnames(m30x)[2]="BindingScore"
  
  #add column with ligand/decoy info
  m10n$IsActive &lt;- as.numeric(m10n$LigandName %in% lig)
  m10n$IsActive = ifelse(m10n$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), 
  m10n$IsActive+1,m10n$IsActive*0) 
  
  m10x$IsActive &lt;- as.numeric(m10x$LigandName %in% lig)
  m10x$IsActive = ifelse(m10x$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m10x$IsActive+1,m10n$IsActive*0)
  
  m30n$IsActive &lt;- as.numeric(m30n$LigandName %in% lig)
  m30n$IsActive = ifelse(m30n$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m30n$IsActive+1,m10n$IsActive*0)
  
  m30x$IsActive &lt;- as.numeric(m30x$LigandName %in% lig)
  m30x$IsActive = ifelse(m30x$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m30x$IsActive+1,m10n$IsActive*0)
  
  #define ROC parameters 
  #here pred is selected to compare between ligands using AutoDock
  PREDm10n &lt;- prediction(m10n$BindingScore*-1, m10n$IsActive {class(x) = "Prediction"; return(x)})
  PERFm10n &lt;- performance(PREDm10n, 'tpr','fpr')
  
  PREDm10x &lt;- prediction(m10x$BindingScore*-1, m10x$IsActive)
  PERFm10x &lt;- performance(PREDm10x, 'tpr','fpr')
  
  PREDm30n &lt;- prediction(m30n$BindingScore*-1, m30n$IsActive)
  PERFm30n &lt;- performance(PREDm30n, 'tpr','fpr')
  
  PREDm30x &lt;- prediction(m30x$BindingScore*-1, m30x$IsActive)
  PERFm30x &lt;- performance(PREDm30x, 'tpr','fpr')
  
  jpeg("Medium Runs.jpg")
  #plot m10n
  m10nsemilog=PERFm10n@x.values[[1]]
  m10nsemilog[m10nsemilog &lt; 0.0005]=0.0005
  plot(m10nsemilog,PERFm10n@y.values[[1]],type="l",lwd=2,xlab="False Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="blue",main="Medium Runs")
  axis(1, c(0,0.001,0.01,0.1,1))
  x&lt;-seq(0,1,0.001)
  points(x,x,col="gray",type="l")
  
  #plotm10x
  m10xsemilog=PERFm10x@x.values[[1]]
  m10xsemilog[m10xsemilog &lt; 0.0005]=0.0005
  lines(m10xsemilog,PERFm10x@y.values[[1]],lwd=2,type="l",xlab="False 
  Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", 
  col="red",main="Short Runs")
  axis(1, c(0,0.001,0.01,0.1,1))
  x&lt;-seq(0,1,0.001)
  points(x,x,col="gray",type="l")
  
  #plotm30n
  m30nsemilog=PERFm30n@x.values[[1]]
  m30nsemilog[m30nsemilog &lt; 0.0005]=0.0005
  lines(m30nsemilog,PERFm30n@y.values[[1]],lwd=2,type="l",xlab="False Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="green",main="Short Runs")
  axis(1, c(0,0.001,0.01,0.1,1))
  x&lt;-seq(0,1,0.001)
  points(x,x,col="gray",type="l")
  
  #plotm30x
  m30xsemilog=PERFm30x@x.values[[1]]
  m30xsemilog[m30xsemilog &lt; 0.0005]=0.0005
  lines(m30xsemilog,PERFm30x@y.values[[1]],lwd=2,type="l",xlab="False 
  Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="purple",main="Short Runs")
  axis(1, c(0,0.001,0.01,0.1,1))
  x&lt;-seq(0,1,0.001)
  points(x,x,col="gray",type="l")
  
  #make legend
  legend("topleft", c("medium10-nmr (AUC = 0.586)", "medium10-xray (AUC = 0.564)", "medium30-nmr (AUC = 0.577)", "medium30-xray (AUC= 0.543)", "random"), lty=1, 
     col = c("red", "blue", "green", "purple", "gray"), bty="n")
  
  dev.off()
  </code></pre>