Essays, Posts & Presentations

John Snow & the Birth of Epidemiology: Data Analysis & Visualization

Roman Parker, Mads Bahrami

Author

Roman Parker, Mads Bahrami

Title

John Snow & the Birth of Epidemiology: Data Analysis & Visualization

Description

Analysis and visualization of cholera case data from the 1854 London outbreak using modern statistical techniques available in the Wolfram Language.

John Snow & the Birth of Epidemiology: Data Analysis & Visualization

August 3, 2021

Roman Parker, Community Advancement, and Mads Bahrami, Community & Content Developer, Community Advancement

In 1854, there was a major cholera outbreak in Soho, a neighborhood in London that Judith Summers described as full of “cow-sheds, animal droppings, slaughterhouses, grease-boiling dens and primitive, decaying sewers.” At the time, the cause of the outbreak was unknown because germ theory was still being developed and disease transmission was not well understood. Miasma theory was the dominant hypothesis, and it proposed that diseases, including cholera and the plague, were spread by foul gasses emitted from decomposing organic matter.

John Snow, a physician known for his work with anesthetics, managed to trace the outbreak back to a water pump on Broad Street contaminated by a nearby cesspit. He did this by mapping the locations of cholera cases and then comparing infections recorded near the Broad Street pump with cases in other neighborhoods. This was a major development in the formulation of germ theory and helped convince many scientists that cholera actually spread through contaminated water instead of noxious air.

We still have Snow’s data on the infected cases and their locations in the Soho neighborhood. In this computational essay, we will explore this data with modern statistical techniques available in the Wolfram Language and the

SpatialPointData

functionality introduced recently.

First, we’ll define the London cholera outbreak as a spatial point configuration:

In[]:=

sp=SpatialPointData

geo-positions



data

region

;

The data contains 322 buildings (locations) in the region that, according to Snow, contained infected individuals (i.e. cases). Additionally, the data includes the number of cases at each location (Cases), how far each location was from the Broad Street pump (DistanceToContaminatedPump) and the nearest uncontaminated pump (DistanceToNonContaminatedPump). The data also contains information on whether the Broad Street pump or a non-contaminated pump was the closest to the location (ContaminatedPumpClosest). The data does not directly use latitude-longitude coordinates, so the pumps’ coordinates must be scaled appropriately.

Now we can define the location of the contaminated Broad Street pump and other non-contaminated pumps in a manner compatible with the data:

In[]:=

contaminatedPump=GeoPosition[{51.51320947964693`,-0.13658234129516839`}];nonContaminatedPumps=

LocationsOfNonContaminatedPumps

;

With the Broad Street pump in red and non-contaminated pumps in green, infection locations are plotted on a map of the area:

In[]:=

RowGeoGraphics[{EdgeForm[Black],sp["ObservationRegion"]},GeoRange

1.5

,ImageSizeMedium,PlotLabel"Infected neighborhood"]," ",Showcases=GeoGraphics[Point@sp["Points"],GeoBackground"VectorMonochrome"],

PumpsMap

,PlotLabel"Houses with some cases; markers are pump locations",ImageSizeMedium

Out[]=

Each black dot is a house with one or more cases, the red marker is the location of the contaminated pump and the green markers are non-contaminated pumps.

We can drill deeper into the data to show the number of cases at each location:

In[]:=

ShowGeoBubbleChart[sp["Points"]sp["Annotations"]["Cases"],BubbleSizes{0.01,0.07},GeoBackground"VectorMonochrome"],

PumpsMap

,PlotLabel"Geo bubble (bubble size corresponds to number of cases per house)"

Out[]=

The data also contains the following annotations:

In[]:=

Column[sp[{"Annotations","Keys"}]]

Out[]=

Cases

ContaminatedPumpClosest

DistanceToContaminatedPump

DistanceToNonContaminatedPump

ContaminatedPumpClosest is Boolean, where “true” means the contaminated pump was the closest pump, and “false” means the closest pump was non-contaminated.

In this dataset, each household (represented by a point) contains at least one case of cholera. We can create a histogram that plots the number of these cases per location to show the proportion:

In[]:=

Histogramsp["Annotations"]["Cases"],

Options



Out[]=

Fit the case counts to a curve and plot them:

In[]:=

edist=EstimatedDistribution[sp["Annotations"]["Cases"],LogSeriesDistribution[a]]

Out[]=

LogSeriesDistribution[0.661006]

Visualize fitting with the histogram:

In[]:=

Show[Histogram[sp["Annotations"]["Cases"],{.5,20.5,1},"PDF",Frame->True],DiscretePlot[PDF[edist,x],{x,1,10,1}]]

Out[]=

Test the data against the obtained distribution:

In[]:=

DistributionFitTest[sp["Annotations"]["Cases"],edist,"ShortTestConclusion"]

Out[]=

Do not reject

The data follows a negative exponential progression. This could be due to the fact that the population per building distribution approximately follows this distribution, assuming the same infection rate per building, on average.

As seen earlier, this dataset also includes the distance to the nearest non-contaminated pump and to the Broad Street pump, which we will explore now.

Create a distribution chart of the non-contaminated and contaminated pump distances:

In[]:=

With{data=Query[{"DistanceToNonContaminatedPump","DistanceToContaminatedPump"}][sp["Annotations"]]},DistributionChartdata,ChartLegendsKeys[data],



Out[]=

	DistanceToNonContaminatedPump
	DistanceToContaminatedPump

This chart shows that because there are more non-contaminated pumps, the maximum distance to them is lower. Cases are, on average, closer to the contaminated pump, and the minimum value of DistanceToContaminatedPump is lower than that of DistanceToNonContaminatedPump.

Show how the cases are distributed in space:

In[]:=

ShowGeoSmoothHistogram[WeightedData[sp["Points"][[1]],sp["Annotations"]["Cases"]],Mesh15,ColorFunction"Rainbow",GeoBackground"VectorMonochrome"],

PumpsMap

,PlotLabel->"Weighted (by Cases) Density Histogram"

Out[]=

This smooth density histogram counts the cases at each location and creates a “density map” of the amount of cases in that area. This shows us that the largest case “hotspot” takes the shape of an elongated oval approximately centered on the Broad Street pump and oriented at about a 45-degree angle from the streets. Additionally, almost all areas with high density are near the contaminated pump. This provides additional evidence that the Broad Street pump is the cause of the cases.

Plot how many points are at each distance from the contaminated pump and the closest non-contaminated pump:

In[]:=

Module{data1=GroupBy[Transpose[{sp["Annotations"]["ContaminatedPumpClosest"],sp["Annotations"]["DistanceToContaminatedPump"]}],FirstRest,Flatten],data2=GroupBy[Transpose[{sp["Annotations"]["ContaminatedPumpClosest"],sp["Annotations"]["DistanceToNonContaminatedPump"]}],FirstRest,Flatten]},LegendedRowBoxWhiskerChartValues@data1,

," ",BoxWhiskerChartValues@data2,

,

Placed[

]



Out[]=

Contaminated Pump Closest

Non-contaminated Pump Closest

Note that in these figures, the interquartile ranges do not overlap. Strangely enough, it appears that some cases used the contaminated Broad Street pump in spite of having a shorter distance to a non-contaminated one. Perhaps the roadways made it easier to get to the Broad Street pump or people were in the area for another reason?

Now we will attempt to find the source of the cholera outbreak. First, we will determine the spatial median. Then, we will calculate how important the relative distance to the Broad Street pump is compared to a non-contaminated pump in affecting cases.

We will do most of this analysis using both unweighted data and data weighted by the number of cases at the location and with different spatial measures. It is very likely that weighting by cases gives a more accurate picture of case hotspots. It is less likely, however, to detect a phenomenon that happens rarely. While unweighted data will detect phenomena like that, it’s easily skewed by outliers.

We can calculate the

SpatialMedian

of the data, weighting it by the number of cases at the location:

In[]:=

weightedSpatialMedian=SpatialMedian[WeightedData[sp["Points"],sp["Annotations"]["Cases"]]];

When we weight the data, the spatial median will be more related to the locations of the cases. This is because areas with very few cases won’t influence the data as much, and case hotspots will provide more of an influence.

When the data is left unweighted, the spatial median is more related to the range and locations of the data because a location with only one case on the edge will provide as much influence as a location with 20 cases toward the center.

The spatial median location minimizes the distance to each case location. Therefore, if the cases were primarily stemming from a single source, it would be likely that the source is near the spatial median. This makes it useful to identify an area to examine more closely for possible sources. Note: unless otherwise specified, the default distance function is

EuclideanDistance

(or

GeoDistance

for geo spatial data).

Find the unweighted spatial median:

In[]:=

unweightedSpatialMedian=SpatialMedian[sp];

Plot the location of the spatial medians (using EuclideanDistance) with the cases:

In[]:=

ShowGeoListPlot[List/@{weightedSpatialMedian,unweightedSpatialMedian},PlotLegends{"Spatial median of weighted data","Spatial median of unweighted data"},GeoBackground"VectorMonochrome"],cases,

PumpsMap



Out[]=

	Spatial median of weighted data
	Spatial median of unweighted data

Let’s use different distance functions to calculate spatial medians:

In[]:=

spList=GeoPosition[SpatialMedian[sp["Points"]〚1〛,DistanceFunction#]]&/@

DistanceFunctions

;

Plot the location of the unweighted spatial medians using different distance functions with the pump location and the cases (note: some of the spatial median locations overlap):

In[]:=

RowShowGeoListPlot[List/@spList,PlotLegendsNone,GeoBackground"VectorMonochrome"],cases,

PumpsMap

,ImageSize->Medium," ",ShowGeoListPlotList/@spList,PlotLegends

DistanceFunctions

,GeoBackground"VectorMonochrome",GeoGraphicsGeoMarkercontaminatedPump,"Color"Directive

,Opacity[.6],EdgeForm[Gray],ImageSize->Medium

Out[]=

	BrayCurtisDistance
	CanberraDistance
	ChessboardDistance
	EuclideanDistance
	ManhattanDistance
	NormalizedSquaredEuclideanDistance
	SquaredEuclideanDistance
	WarpingDistance

The spatial medians, whether weighted, unweighted or with different distance functions, are near the location of the Broad Street pump and each other (see the previous visualizations). This tells us that the number of cases at a location is approximately rotationally symmetric around the spatial median. Otherwise, the unweighted spatial median would be quite different than the spatial median.

Additionally, looking at the following visualization, it appears that the case locations are more common toward the center of the data. This provides some amount of anecdotal evidence that there is a single data source.

Show a smooth density histogram of the data weighted by the number of cases at that location:

In[]:=

ShowGeoSmoothHistogram[WeightedData[sp["Points"][[1]],sp["Annotations"]["Cases"]],Mesh15,ColorFunction"Rainbow",GeoBackground"VectorMonochrome"],

PumpsMap

,PlotLabel->"Weighted (by Cases) Density Histogram"

Out[]=

Show a smooth density histogram of the case locations, ignoring the number of cases:

In[]:=

ShowGeoSmoothHistogram[sp["Points"],Mesh15,ColorFunction"Rainbow",GeoBackground"VectorMonochrome"],

PumpsMap

,PlotLabel->"Unweighted Density Histogram"

Out[]=

Next, we will calculate how much of an effect being near the Broad Street pump instead of another pump had on the number of cases at a location. It may provide additional evidence for the Broad Street pump being responsible.

We calculate a measure that is simply the ratio between DistanceToNonContaminatedPump to DistanceToContaminatedPump. We would expect that as our number of cases increases, at least up to a point, the value of the ratio should also increase because the ratio will increase as a point is closer to the Broad Street pump:

In[]:=

ratio=Rescale[Divide@@Lookup[{"DistanceToNonContaminatedPump","DistanceToContaminatedPump"}][sp["Annotations"]]];

We can create a

BoxWhiskerChart

of cases and ratios to see if it shows a similar phenomenon or if it is unique to this one method of visualization.

Show box-whisker charts of the ratio for each number of cases:

In[]:=

With[{data=KeySort@GroupBy[Transpose[{ratio,sp["Annotations"]["Cases"]}],LastMost,Flatten]},BoxWhiskerChart[data,"Outliers",ChartLabelsKeys[data],PlotRange{0,.1}]]

Out[]=

Cases and ratios seem to rise together up until five cases, and although the rise in the median still seems to be linear, the 75th percentile appears to grow exponentially, which is interesting. Additionally, the lowest value and 25th percentile continue to rise when there are five cases. This suggests that our hypothesis from earlier that the low values are potentially outliers may be correct. We still do not have, however, the evidence to prove this, and there very well could be some secondary factor here impacting case numbers.

Next, we will test to see if the ratio and cases are independent from each other:

In[]:=

IndependenceTest[ratio,sp["Annotations"]["Cases"],"TestConclusion"]

Out[]=

The null hypothesis that the populations are independent is rejected at the 5 percent level based on the Spearman Rank test.

Showing a lack of independence provides another piece of evidence that the ratio, and therefore the distance to the Broad Street pump and other pumps, is related to cases. Between all these visualizations and data analysis, we can say that the Broad Street pump seems more likely to be related to the cholera outbreak, while the other pumps are significantly less likely to be contributors to it.

Snow’s map persuaded the Board of Guardians of St. James’s Parish, in whose parish the pump fell, to remove the handle from the pump. While the overall effectiveness of this tactic remains in doubt as the cholera outbreak was already declining, the impact of Snow’s disease transmission discovery went far beyond the 1854 cholera outbreak.

In fact, Snow could be considered “the father of epidemiology” because, as the US Centers for Disease Control and Prevention (CDC) states on its website, “his work illustrates the classic sequence from descriptive epidemiology to hypothesis generation to hypothesis testing (analytic epidemiology) to application....”

Today, if you visit the site of the offending Broad Street pump, you’ll find a double commemoration of Snow’s work: a memorial in the form of a water pump and his very own “watering hole,” the John Snow pub (as seen thanks to Google Street View).

Acknowledgments

We would like to thank Gosia Konwerska (Wolfram Research) for providing the data and helping to improve the essay and Jon McLoone (Wolfram Research) for his help on transforming data into actual geo positions and providing historical and contemporary context about John Snow’s contributions to epidemiology.

Cite this as: Roman Parker, Mads Bahrami, "John Snow & the Birth of Epidemiology: Data Analysis & Visualization" from the Notebook Archive (2022), https://notebookarchive.org/2022-01-9pj6xxp