Bayesian calculation for twin Down Syndrome probabilities (Flask app)

app.py

from flask import Flask, request, url_for
from scipy.optimize import bisect
import math
import numpy as np

app = Flask(__name__)
app.secret_key = 'Whatever'

@app.route('/')
def index():
    return """A Bayesian <a href="%s">calculator</a>
              for Down Syndrome odds in twin pregnancies.""" % (url_for('twin_calculator_input'),)

def mathjax():
    return """<script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript">
        </script>"""

def introduction():
    return """
       This calculator might be of interest to you if:
       <ol>
       <li> you are expecting twins and
       <li> received a "positive" test result for Down Syndrome.
       </ol>
       But I am not a doctor (not the right kind, anyway) so you should seek medical advice.
      <p>
    <b>Computing the probability of both your twins having Down Syndrome</b> (idealized model)
         <p>
        Because the blood tests cannot distinguish between the two foetuses, the results are
        difficult to interpret. Yet they tend to be communicated in the same manner as results for single
        foetus pregnancies - which I personally find confusing. What is intended here is a calculation to estimate the probability of giving birth to
      zero, one or two Down Syndrome children in a twin pregnancy after a positive blood test for Down Syndrome has
        been received. There are many caveats, but let's proceed...
      """

def input_form():
    return """
        <p>
        <b> User inputs </b>
        <form method = "POST" action="%s">
            The calculation requires your age: <input type = "number" value = "35" min = "20" max = "49" step="1" name="age" /> (please adjust).
            <p>
             We will also need a probability of identical twins equal to
                <input type = "number" value = "0.1" min = "0.01" step="0.01" name="prob_identical" /> (please adjust).
            This should be a reasonable estimate of the probability that your twins are identical, not fraternal.
         In other words, based on previous examinations (and not taking into account the blood test, which would be circular), roughly what probability
         did your ostetrician assign to your twins being identical? Please enter a number between 0 and 1.
         (The default setting of 0.1 is very roughly the population incidence of identical versus fraternal
         twins. However it is usually possible for a medical professional to form an opinion in your individual case.)
         <p>
         We will assume Down Syndrome odds of 1 in <input type = "number" value = "150" min = "1" step="1" name="odds_down" /> (please adjust).
         This is the number that was provided in the integrated test. (If your test was "positive" it generally
         implies the odds are 200/1 or less, but the terminology "positive" seems largely meaningless without the number itself).
         <p>
         Once you have adjusted these inputs, please <input type="submit" value="submit" />.
         </form>
         """ % (url_for('twin_calculator_output'),)


@app.route('/twins_calculator_input')
def twin_calculator_input():
    return mathjax() + introduction() + input_form()

def priorOdds( age ):
    # Taken from http://downsyndrome.about.com/od/diagnosingdownsyndrome/a/Matagechart.htm
    table = { 20:1667, 21:1429,22:1429, 23:1429, 24:1250,
              25:1250, 26:1176,27:1111, 28:1053, 29:1000,
              30: 952, 31:909, 32:769,  33:625,  34: 500,
              35: 385, 36:294, 37:227,  38:175,  39:137,
              40: 106, 41:82,  43:50,   44:38,   45:30,
              46: 23,  47:18,  48:18,   49:11 }
    return table[age]

muNormal = 0.0
muDown   = 1.0

def implyMarker( odds_ratio ):
    """ Given odds of n:1 , interpret as a marker measurement """

    def markerRatio(x ):
         return math.exp( -(x-muDown)*(x-muDown)/2. ) / math.exp( -(x-muNormal)*(x-muNormal)/2. )
    def markerRatioError( x ):
         return markerRatio(x) - odds_ratio

    a  = muNormal - 10.0
    b  = muDown + 10.0

    marker = bisect( markerRatioError, a, b)
    return marker

@app.route('/twins_calculator_output', methods = ['POST'])
def twin_calculator_output():

    # Interpret the odds ratio supplied as a lab result
    age            = int(request.form["age"])
    odds_down      = float(request.form["odds_down"])
    odds_prior     = priorOdds( age )
    odds_ratio     = odds_prior / (2*odds_down)
    prob_identical = float(request.form["prob_identical"])
    marker         = implyMarker( odds_ratio )

    # Posterior probabilities for the four possibilities
    prob_down      = 1/(2*odds_down)
    pDD            = prob_identical*prob_down + (1-prob_identical)*prob_down*prob_down
    pDU            = (1-prob_identical)*prob_down*(1-prob_down)
    pUD            = (1-prob_identical)*prob_down*(1-prob_down)
    pUU            = (1-prob_identical)*(1-prob_down)*(1-prob_down) + prob_identical*(1-prob_down)
    probs          = [ pDD, pDU, pUD, pUU ]
    markers        = [ (muDown,muDown), (muDown,muNormal), (muNormal,muDown), (muNormal,muNormal) ]
    posteriors     = list()
    for (marker1, marker2),prior in zip( markers, probs ):
        meanMarker = 0.5*marker1 + 0.5*marker2
        prob       = math.exp( -( marker-(meanMarker )*( marker-meanMarker )/2. ) )
        posterior  = prior*prob
        posteriors.append( posterior )

    p0 = posteriors[3] / np.sum( posteriors )
    p1 = (posteriors[1]+posteriors[2]) / np.sum( posteriors )
    p2 = posteriors[0] / np.sum( posteriors )

    return """ <b> Results </b>
               <p>
               Based on an age of %s the prior odds of Down Syndrome are %s:1.
               <p>
               So the updated odds of %s:1 might (I presume) be interpreted as odds of
               roughly %s:1 for each foetus and thus
               as an odds ratio of %s and, furthermore, as
               a marker measurement of %s in standardized terms.
               <p>
               However, assuming a fifty/fifty mix in the blood test (i.e. half coming from each foetus, which
               is admittedly quite an assumption) the
               twin probabilities may actually be as follows.
               <p>
               <table border="1" cellpadding="3" cellspacing="0" width="300">
               <thead>
                  <tr>
                     <th>Outcome</th>
                     <th>Probability</th>
                  </tr>
               </thead>
               <tbody>
                  <tr>
                     <th> Neither has Down Syndrome</th>
                     <td> %s percent</td>
                  </tr>
                  <tr>
                     <th> One has Down Syndrome</th>
                     <td> %s percent</td>
                  </tr>
                  <tr>
                     <th> Both have Down Syndrome</th>
                     <td> 1 in %s </td>
                  </tr>
               </tbody>
            </table>
            <p>
             These probabilites should be taken with a grain of salt. There are many assumptions, including
             a prior probability of identical twins equal to %s percent. Perhaps the biggest simplification, however,
             is that the test result is from a single marker. In reality it is from multiple markers.
               <p>
               As a final caveat, I may well have misinterpreted the odds ratio supplied to patients
               who have twins. Here I have assumed that the same scale is used for single pregnancies in converting
               from marker levels to an odds ratio. I'd be happy to be corrected on that one.
               <p>
               Please provide feedback at github where the
               <a href="https://gist.github.com/notbanker/223652684d8c0b270b40">code</a> for this little app is posted.

              """ % ( age, odds_prior, odds_down, 2*odds_down, int(10*odds_ratio)/10., int(100*marker)/100.,
                      int(p0*1000)/10., int(p1*1000)/10., int(1/p2), int(prob_identical*100) )