Skip to content

Instantly share code, notes, and snippets.

Forked from notbanker/
Created Oct 23, 2016
What would you like to do?
Bayesian calculation for twin Down Syndrome probabilities (Flask app)
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'
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="" type="text/javascript">
def introduction():
return """
This calculator might be of interest to you if:
<li> you are expecting twins and
<li> received a "positive" test result for Down Syndrome.
But I am not a doctor (not the right kind, anyway) so you should seek medical advice.
<b>Computing the probability of both your twins having Down Syndrome</b> (idealized model)
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 """
<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).
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.)
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).
Once you have adjusted these inputs, please <input type="submit" value="submit" />.
""" % (url_for('twin_calculator_output'),)
def twin_calculator_input():
return mathjax() + introduction() + input_form()
def priorOdds( age ):
# Taken from
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>
Based on an age of %s the prior odds of Down Syndrome are %s:1.
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.
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.
<table border="1" cellpadding="3" cellspacing="0" width="300">
<th> Neither has Down Syndrome</th>
<td> %s percent</td>
<th> One has Down Syndrome</th>
<td> %s percent</td>
<th> Both have Down Syndrome</th>
<td> 1 in %s </td>
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.
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.
Please provide feedback at github where the
<a href="">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) )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment