Statistics review 13: Receiver operating characteristic
=======================================================

R code accompanying
`paper <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1065080/pdf/cc3000.pdf>`__

Key learning points
-------------------

-  Methods for assessing the performance of a diagnostic test
-  Sensitivity, specificity and likelihood ratio of a test
-  Receiver operating characteristic curve and the area under the curve

.. code:: r

    suppressPackageStartupMessages(library(tidyverse))

.. code:: r

    options(repr.plot.width=4, repr.plot.height=3)

Data
----

.. code:: r

    died <- c(81, 45)
    survived <- c(591, 674)
    
    df1 <- data.frame(died=died, survived=survived,
                      row.names=c("Lactate >1.5 mmol/l", "Lactate ≤1.5 mmol/l"))
    df1



.. raw:: html

    <table>
    <thead><tr><th></th><th scope=col>died</th><th scope=col>survived</th></tr></thead>
    <tbody>
    	<tr><th scope=row>Lactate &gt;1.5 mmol/l</th><td>81 </td><td>591</td></tr>
    	<tr><th scope=row>Lactate ≤1.5 mmol/l</th><td>45 </td><td>674</td></tr>
    </tbody>
    </table>



Data set with different prevalence
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    died <- c(386, 214)
    survived <- c(370, 421)
    
    df2 <- data.frame(died=died, survived=survived,
                      row.names=c("Lactate >1.5 mmol/l", "Lactate ≤1.5 mmol/l"))
    df2



.. raw:: html

    <table>
    <thead><tr><th></th><th scope=col>died</th><th scope=col>survived</th></tr></thead>
    <tbody>
    	<tr><th scope=row>Lactate &gt;1.5 mmol/l</th><td>386</td><td>370</td></tr>
    	<tr><th scope=row>Lactate ≤1.5 mmol/l</th><td>214</td><td>421</td></tr>
    </tbody>
    </table>



Sensitivity and specificity
---------------------------

.. code:: r

    prop.ci <- function(p, n, alpha=0.05) {
        se <- sqrt(p*(1-p)/n)
        k <- qnorm(1 - alpha/2)
        round(c(p - k*se, p + k*se), 2)
    }

.. code:: r

    sens <- function(tbl) {
        a <- tbl[1,1]
        b <- tbl[1,2]
        c <- tbl[2,1]
        d <- tbl[2,2]
        n <- a + c   
        p <- a/n
       
        round(p, 2)
    }

.. code:: r

    sens(df1)



.. raw:: html

    0.64


.. code:: r

    p <- sens(df1)
    n <- sum(df1)
    prop.ci(p, n)



.. raw:: html

    <ol class=list-inline>
    	<li>0.61</li>
    	<li>0.67</li>
    </ol>



.. code:: r

    spec <- function(tbl) {
        a <- tbl[1,1]
        b <- tbl[1,2]
        c <- tbl[2,1]
        d <- tbl[2,2]
        n <- b + d
        
        p <- d/n
        round(p, 2)
    }

.. code:: r

    spec(df1)



.. raw:: html

    0.53


Positive and negative predictive values
---------------------------------------

.. code:: r

    ppv <- function(tbl) {
        a <- tbl[1,1]
        b <- tbl[1,2]
        c <- tbl[2,1]
        d <- tbl[2,2]
        n <- a + b
        p <- a/n
        
        round(p, 2)
    }

.. code:: r

    ppv(df1)



.. raw:: html

    0.12


.. code:: r

    npv <- function(tbl, alpha=0.05) {
        a <- tbl[1,1]
        b <- tbl[1,2]
        c <- tbl[2,1]
        d <- tbl[2,2]
        n <- c + d
        
        p <- d/n    
        round(p, 2)
    }

.. code:: r

    npv(df1)



.. raw:: html

    0.94


Effect of prevalence
--------------------

.. code:: r

    c(sens(df1)[1], sens(df2)[1])



.. raw:: html

    <ol class=list-inline>
    	<li>0.64</li>
    	<li>0.64</li>
    </ol>



.. code:: r

    c(spec(df1)[1], spec(df2)[1])



.. raw:: html

    <ol class=list-inline>
    	<li>0.53</li>
    	<li>0.53</li>
    </ol>



.. code:: r

    c(ppv(df1)[1], ppv(df2)[1])



.. raw:: html

    <ol class=list-inline>
    	<li>0.12</li>
    	<li>0.51</li>
    </ol>



.. code:: r

    c(npv(df1)[1], npv(df2)[1])



.. raw:: html

    <ol class=list-inline>
    	<li>0.94</li>
    	<li>0.66</li>
    </ol>



Likelihood ratios
-----------------

.. code:: r

    lr <- function(tbl) {
        round(sens(tbl)[1]/(1 - spec(tbl)[1]), 2)
    }

.. code:: r

    lr(df1)



.. raw:: html

    1.36


Receiver operating characteristic curve (ROC)
---------------------------------------------

.. code:: r

    threshold <- c(0,1,1.5,2,3,5,25)
    died <- c(126,114,81,58,37,19,0)
    survived <- c(1265,996,591,329,131,27,0)
    
    df3 <- data.frame(threshold=threshold, died=died, survived=survived)
    df3



.. raw:: html

    <table>
    <thead><tr><th scope=col>threshold</th><th scope=col>died</th><th scope=col>survived</th></tr></thead>
    <tbody>
    	<tr><td> 0.0</td><td>126 </td><td>1265</td></tr>
    	<tr><td> 1.0</td><td>114 </td><td> 996</td></tr>
    	<tr><td> 1.5</td><td> 81 </td><td> 591</td></tr>
    	<tr><td> 2.0</td><td> 58 </td><td> 329</td></tr>
    	<tr><td> 3.0</td><td> 37 </td><td> 131</td></tr>
    	<tr><td> 5.0</td><td> 19 </td><td>  27</td></tr>
    	<tr><td>25.0</td><td>  0 </td><td>   0</td></tr>
    </tbody>
    </table>



.. code:: r

    n.died <- 126
    n.survived <- 1265

.. code:: r

    sens.1 <-function(died, ndied) {
        a <- died
        c <- ndied-died
        a/(a + c)
    }

.. code:: r

    spec.1 <-function(survived, nsurvived) {
        b <- survived
        d <- nsurvived - survived
        b/(b+d)
    }

.. code:: r

    df4 <- df3 %>% mutate(sens=sens.1(died, n.died), one.minus.spec=spec.1(survived, n.survived)) 
    df4



.. raw:: html

    <table>
    <thead><tr><th scope=col>threshold</th><th scope=col>died</th><th scope=col>survived</th><th scope=col>sens</th><th scope=col>one.minus.spec</th></tr></thead>
    <tbody>
    	<tr><td> 0.0      </td><td>126       </td><td>1265      </td><td>1.0000000 </td><td>1.00000000</td></tr>
    	<tr><td> 1.0      </td><td>114       </td><td> 996      </td><td>0.9047619 </td><td>0.78735178</td></tr>
    	<tr><td> 1.5      </td><td> 81       </td><td> 591      </td><td>0.6428571 </td><td>0.46719368</td></tr>
    	<tr><td> 2.0      </td><td> 58       </td><td> 329      </td><td>0.4603175 </td><td>0.26007905</td></tr>
    	<tr><td> 3.0      </td><td> 37       </td><td> 131      </td><td>0.2936508 </td><td>0.10355731</td></tr>
    	<tr><td> 5.0      </td><td> 19       </td><td>  27      </td><td>0.1507937 </td><td>0.02134387</td></tr>
    	<tr><td>25.0      </td><td>  0       </td><td>   0      </td><td>0.0000000 </td><td>0.00000000</td></tr>
    </tbody>
    </table>



.. code:: r

    ggplot(df4, aes(x=one.minus.spec, y=sens)) + geom_point() + geom_line() +
    geom_abline(intercept = 0, slope = 1, color='grey', linetype='dashed') +
    coord_fixed()





.. image:: SR13_Receiver_operating_characteristic_files/SR13_Receiver_operating_characteristic_33_1.png


Area under the curve (AUCROC)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    suppressPackageStartupMessages(require(pracma))

.. code:: r

    round(-1 * trapz(x=df4$one.minus.spec, y=df4$sens), 2)



.. raw:: html

    0.64