SR11_Assessing_risk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistics review 11: Assessing risk\n",
    "\n",
    "R code accompanying [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC522855/pdf/cc2908.pdf)\n",
    "\n",
    "## Key learning points\n",
    "\n",
    "- Relative Risk\n",
    "- Odds Ratio\n",
    "- Measuring the impact of exposure to a risk factor\n",
    "- Measures of the success of a treatment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "suppressPackageStartupMessages(library(tidyverse))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "That study investigated the association between surfactant protein B and acute respiratory distress syndrome (ARDS). Patients were classified according to their thymine/cytosine (C/T) gene coding, and patients with the C allele present (genotype CC or CT) were compared with those with genotype TT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>ARDS</th><th scope=col>NoARDS</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>CC/CT</th><td>11 </td><td>208</td></tr>\n",
       "\t<tr><th scope=row>TT</th><td> 1 </td><td>182</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & ARDS & NoARDS\\\\\n",
       "\\hline\n",
       "\tCC/CT & 11  & 208\\\\\n",
       "\tTT &  1  & 182\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | ARDS | NoARDS | \n",
       "|---|---|\n",
       "| CC/CT | 11  | 208 | \n",
       "| TT |  1  | 182 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "      ARDS NoARDS\n",
       "CC/CT 11   208   \n",
       "TT     1   182   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ARDS <- c(11,1)\n",
    "NoARDS <- c(208, 182)\n",
    "df <- data.frame(row.names=c(\"CC/CT\", \"TT\"), ARDS=ARDS, NoARDS=NoARDS)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relative risk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>ARDS</th><th scope=col>NoARDS</th><th scope=col>Total</th><th scope=col>risk</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>11         </td><td>208        </td><td>219        </td><td>0.050228311</td></tr>\n",
       "\t<tr><td> 1         </td><td>182        </td><td>183        </td><td>0.005464481</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       " ARDS & NoARDS & Total & risk\\\\\n",
       "\\hline\n",
       "\t 11          & 208         & 219         & 0.050228311\\\\\n",
       "\t  1          & 182         & 183         & 0.005464481\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "ARDS | NoARDS | Total | risk | \n",
       "|---|---|\n",
       "| 11          | 208         | 219         | 0.050228311 | \n",
       "|  1          | 182         | 183         | 0.005464481 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  ARDS NoARDS Total risk       \n",
       "1 11   208    219   0.050228311\n",
       "2  1   182    183   0.005464481"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df1 <- df %>% mutate(Total = ARDS+NoARDS) %>% mutate(risk = ARDS/Total)\n",
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.19"
      ],
      "text/latex": [
       "9.19"
      ],
      "text/markdown": [
       "9.19"
      ],
      "text/plain": [
       "[1] 9.19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rr <- df1$risk[1]/df1$risk[2]\n",
    "round(rr, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rr <- function(tbl) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    round((a/(a+b))/(c/(c+d)), 2)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "rr_ci <- function(tbl, alpha=0.05) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    log.se <- sqrt(1/a - 1/(a+b) + 1/c - 1/(c+d))\n",
    "    log.rr <- log((a/(a+b))/(c/(c+d)))\n",
    "    k <- qnorm(1-alpha/2)\n",
    "    log.ci <- c(log.rr - k*log.se, log.rr + k*log.se)\n",
    "    round(exp(log.ci), 2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.19"
      ],
      "text/latex": [
       "9.19"
      ],
      "text/markdown": [
       "9.19"
      ],
      "text/plain": [
       "[1] 9.19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rr(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>1.2</li>\n",
       "\t<li>70.53</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1.2\n",
       "\\item 70.53\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1.2\n",
       "2. 70.53\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]  1.20 70.53"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rr_ci(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Odds ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>ARDS</th><th scope=col>NoARDS</th><th scope=col>odds</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>11         </td><td>208        </td><td>0.052884615</td></tr>\n",
       "\t<tr><td> 1         </td><td>182        </td><td>0.005494505</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " ARDS & NoARDS & odds\\\\\n",
       "\\hline\n",
       "\t 11          & 208         & 0.052884615\\\\\n",
       "\t  1          & 182         & 0.005494505\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "ARDS | NoARDS | odds | \n",
       "|---|---|\n",
       "| 11          | 208         | 0.052884615 | \n",
       "|  1          | 182         | 0.005494505 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  ARDS NoARDS odds       \n",
       "1 11   208    0.052884615\n",
       "2  1   182    0.005494505"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df2 <- df %>% mutate(odds = ARDS/NoARDS)\n",
    "df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.625"
      ],
      "text/latex": [
       "9.625"
      ],
      "text/markdown": [
       "9.625"
      ],
      "text/plain": [
       "[1] 9.625"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "or <- df2$odds[1]/df2$odds[2]\n",
    "or"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "or <- function(tbl) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    round((a/b)/(c/d), 2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "or_ci <- function(tbl, alpha=0.05) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    log.se <- sqrt(1/a + 1/b + 1/c - 1/d)\n",
    "    log.rr <- log((a/b)/(c/d))\n",
    "    k <- qnorm(1-alpha/2)\n",
    "    log.ci <- c(log.rr - k*log.se, log.rr + k*log.se)\n",
    "    round(exp(log.ci), 2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.62"
      ],
      "text/latex": [
       "9.62"
      ],
      "text/markdown": [
       "9.62"
      ],
      "text/plain": [
       "[1] 9.62"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "or(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>1.24</li>\n",
       "\t<li>74.5</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1.24\n",
       "\\item 74.5\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1.24\n",
       "2. 74.5\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]  1.24 74.50"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "or_ci(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advantages of odds ratio\n",
    "\n",
    "- Can be estimated in case-control study\n",
    "- OR is a symmetric ratio in that the OR for the disease given the risk factor is the same as the OR for the risk factor given the disease.\n",
    "- form part of the output when carrying out logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Attributable risk\n",
    "\n",
    "The proportion of cases in a population that could be prevented if the risk factor were to be eliminated. The AR is the difference between the actual number of cases in a sample and the number of cases that would be expected if exposure to the risk factor were eliminated, expressed as a proportion of the former."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ar <- function(tbl) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    n <- a + b + c + d\n",
    "    ar <- ((a+c) - (n*c)/(c+d))/(a+c)   \n",
    "    round(100*ar, 2)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ar_ci <- function(tbl, alpha=0.05) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "    n <- a + b + c + d\n",
    "    k <- qnorm(1-alpha/2)\n",
    "    u <- (k*(a+c)*(c+d))/(a*d-b*c) * sqrt(((a*d*(n-c) + c^2*b)/(n*c*(a+c)*(c+d))))\n",
    "    \n",
    "    hi <- ((a*d - b*c)*exp(u))/(n*c + (a*d-b*c)*exp(u))\n",
    "    lo <- ((a*d - b*c)*exp(-u))/(n*c + (a*d-b*c)*exp(-u)) \n",
    "    round(100*c(lo, hi), 2)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "81.69"
      ],
      "text/latex": [
       "81.69"
      ],
      "text/markdown": [
       "81.69"
      ],
      "text/plain": [
       "[1] 81.69"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ar(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>31.16</li>\n",
       "\t<li>97.78</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 31.16\n",
       "\\item 97.78\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 31.16\n",
       "2. 97.78\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 31.16 97.78"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ar_ci(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Risk measurements in clinical trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>survived</th><th scope=col>died</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>early</th><td>79</td><td>38</td></tr>\n",
       "\t<tr><th scope=row>standard</th><td>60</td><td>59</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & survived & died\\\\\n",
       "\\hline\n",
       "\tearly & 79 & 38\\\\\n",
       "\tstandard & 60 & 59\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | survived | died | \n",
       "|---|---|\n",
       "| early | 79 | 38 | \n",
       "| standard | 60 | 59 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "         survived died\n",
       "early    79       38  \n",
       "standard 60       59  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df3 <- data.frame(survived=c(79, 60), died=c(38, 59), \n",
    "                  row.names=c(\"early\", \"standard\"))\n",
    "df3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.34"
      ],
      "text/latex": [
       "1.34"
      ],
      "text/markdown": [
       "1.34"
      ],
      "text/plain": [
       "[1] 1.34"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rr(df3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.04"
      ],
      "text/latex": [
       "2.04"
      ],
      "text/markdown": [
       "2.04"
      ],
      "text/plain": [
       "[1] 2.04"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "or(df3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "14.39"
      ],
      "text/latex": [
       "14.39"
      ],
      "text/markdown": [
       "14.39"
      ],
      "text/plain": [
       "[1] 14.39"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ar(df3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Risk difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arr <- function(tbl) {\n",
    "    a <- tbl[1,1]\n",
    "    b <- tbl[1,2]\n",
    "    c <- tbl[2,1]\n",
    "    d <- tbl[2,2]\n",
    "     \n",
    "    r <- (d/(c+d) - b/(a+b))\n",
    "    round(100*r, 2)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "17.1"
      ],
      "text/latex": [
       "17.1"
      ],
      "text/markdown": [
       "17.1"
      ],
      "text/plain": [
       "[1] 17.1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "arr(df3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number needed to treat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nnt <- function(df) {\n",
    "    round(100/arr(df), 0)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "6"
      ],
      "text/latex": [
       "6"
      ],
      "text/markdown": [
       "6"
      ],
      "text/plain": [
       "[1] 6"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nnt(df3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1**. Write a function to calculate the confidence intervals of the absolute risk reduction (ARR) given a $2 \\times 2$ table of outcomes. What is the 90% CI for the data given in `df3`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2** Write a function to calculate the confidence intervals for the number needed to treat (NNT) 2×2\n",
    "  table of outcomes. What is the 90% CI for the data given in df3?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}