Fetching the data
In this blog post I want to look at visualization, and basic transformations that serve visualization with pandas. I have multiple goals with this post:
- Demonstrate some practical use cases for simple data pocessing and data visualization with pandas.
- Show how the stats for COVID-19 mortality, so far, compare to those of the Spanish flu.
The public data set used for this analysis is the deaths csv file from the Novel Corona Virus 2019 Dataset.
imports
We start our analysis by importing the nessasary Python libraries into Jupyter Notebook.
#The pandas dataframe library
import pandas as pd
#We use numpy because sometimes exponential data need a logaritmic scale
import numpy as np
#matplotlib for plotting our data
import matplotlib.pyplot as plt
# urllib for fetching web pages from wikipedia
import urllib.request as rq
# BeautifulSoup for parsing wikipedia HTML pages from wikipedia.
from bs4 import BeautifulSoup as bs
%matplotlib inline
A first transformation
The data in the csv file has two columns "Country/Region" and "Province/State" that together act as a unique key identifying either a whole country or a specific region within that country. The mortality data for those countries and regions is contained in many per day columns.
For visualization we want to transform this data. This transformation could be done in multiple ways. We could use pandas its transpose function for data frames, and then in other steps remove coordinates, combine the two part index into a single one and shift timelines to all start on the first COVID-19 death. Instead we choose to itterate the region rows and the date columns from the 20th of January till the 30th of march, in the original data frame and use that data to create a new transposed dataframe.
We put out transformation code into the constructor of a simple Python object and add some methods we'll be using later on.
class DataSet(object):
def __init__(self):
deaths = pd.read_csv("time_series_covid_19_deaths.csv")
deaths.fillna('', inplace=True)
trydates = []
for day in range(22,32):
trydates.append("1/" + str(day) + "/20")
for day in range(1,30):
trydates.append("2/" + str(day) + "/20")
for day in range(1,31):
trydates.append("3/" + str(day) + "/20")
recdata_mort = dict()
deathset = set()
for index, row in deaths.iterrows():
province = row["Province/State"]
country = row["Country/Region"]
deathset.add(country + ":" + province)
for region in deathset:
cy, province = region.split(":")
cydeaths = deaths[deaths["Country/Region"] == cy]
cydeaths = cydeaths[cydeaths["Province/State"] == province]
results_mort = []
lastdth = 0
for day in trydates:
try:
dth = cydeaths[day].tolist()[0]
except:
dth = lastdth
if dth > 0:
results_mort.append(dth)
lastdth = dth
recdata_mort[region] = results_mort[:]
dfdata_mort = list()
for index in range(0,len(trydates)):
obj_mort = dict()
for key in deathset:
try:
obj_mort[key] = recdata_mort[key][index]
except:
pass
dfdata_mort.append(obj_mort)
self.timelines_mort = pd.DataFrame(dfdata_mort)
def deaths(self):
return self.timelines_mort
def setset(self, popset):
cyset = set()
for index, row in popset.iterrows():
if ":" not in row["cy"]:
cyset.add(row["cy"] + ":")
else:
cyset.add(row["cy"])
return set(self.timelines_mort.columns) - cyset
def dropother(self, popset):
cyset = set()
for index, row in popset.iterrows():
if ":" not in row["cy"]:
cyset.add(row["cy"] + ":")
else:
cyset.add(row["cy"])
dropset = set(self.timelines_mort.columns) - cyset
rframe = self.timelines_mort.copy()
for dropcy in dropset:
print("dropping", dropcy)
del rframe[dropcy]
return rframe
The setset and dropother method methods both take a data frame with per region population size info. The setset method returns a list of mortality data regions that are missing from the provided population size data frame. The dropother method returns a stripped down version of the transformed mortality data frame with just those regions and countries that are also in the population size data frame. We will revisit these later.
The first thing we do now is instantiate a DataSet object from the class above.
df = DataSet()
Getting a feel for the data
Let start by simply plotting the data as is. Given the huge diferences in population sizes and curve starts, the initial view on the data won't be all that usefull, but it gives us a rough idea what we are looking at.
deathcomul = df.deaths()
deathcomul.plot(figsize=(10,10), legend=False)
<matplotlib.axes._subplots.AxesSubplot at 0x7f225f2170b8>
Because much of the data curves are hidden from view by the scale of the mortality numbers in a few countries and/or regions, and because we know part of many of the curves theoretically start of with exponential growth, we have a look at the log of the same data.
ldeathcomul = np.log(df.deaths() +1)
/home/rob/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in log
"""Entry point for launching an IPython kernel.
ldeathcomul.plot(figsize=(10,10), legend=False)
<matplotlib.axes._subplots.AxesSubplot at 0x7f225c856d68>
Collecting population sizes.
Wikipedia has a number of pages with tables of per country population data. Here is a wikipedia page with countries and their population sizes. We fetch this file and see if we can read it into a data frame. We need to do some HTML traversal using the Python BeautifulSoup library.
Note that some countries require mapping from Wikipedia names to the names used in our COVID-19 data file.
wikipedia = 'https://en.wikipedia.org/wiki/'
page1 = rq.urlopen(wikipedia + 'List_of_countries_and_dependencies_by_population')
s1 = bs(page1, 'html.parser')
t1 = s1.find_all("table", attrs={"class": "wikitable sortable"})[0]
altnames = dict()
altnames["Myanmar"] = "Burma"
altnames["Cape Verde"] = "Cabo Verde"
altnames["Ivory Coast"] = "Cote d'Ivoire"
altnames["Congo"] = "Congo (Brazzaville)"
altnames["DR Congo"] = "Congo (Kinshasa)"
altnames["Czech Republic"] = "Czechia"
altnames["United States"] = "US"
altnames["East Timor"] = "Timor-Leste"
altnames["Taiwan"] = "Taiwan*"
altnames["South Korea"] = "Korea, South"
objlist = list()
for row in t1.find_all('tr')[1:-1]:
tds = row.find_all('td')
links = tds[1].find_all('a')
count = len(links)
candidate = links[0].getText()
if count > 1:
alt = links[1].getText()
if "[" not in alt:
if not candidate:
candidate = alt
else:
candidate = alt + ":" + candidate
if candidate in altnames.keys():
print("Altname for", candidate, ":", altnames[candidate])
candidate = altnames[candidate]
if ":" not in candidate:
candidate += ":"
obj = dict()
obj["cy"] = candidate
obj["population"] = int("".join(tds[2].getText().split(",")))
if obj["cy"]:
objlist.append(obj)
populations = pd.DataFrame(objlist)
populations.head(6)
Altname for United States : US
Altname for DR Congo : Congo (Kinshasa)
Altname for Myanmar : Burma
Altname for South Korea : Korea, South
Altname for Ivory Coast : Cote d'Ivoire
Altname for Taiwan : Taiwan*
Altname for Czech Republic : Czechia
Altname for Congo : Congo (Brazzaville)
Altname for East Timor : Timor-Leste
Altname for Cape Verde : Cabo Verde
| cy | population | |
|---|---|---|
| 0 | China: | 1402001120 |
| 1 | India: | 1360456490 |
| 2 | US: | 329544974 |
| 3 | Indonesia: | 266911900 |
| 4 | Pakistan: | 219126520 |
| 5 | Brazil: | 211327887 |
Now we use the setset method we discussed at the beginning of this post, to see what countries and regions are still missing.
df.setset(populations)
{'Australia:Australian Capital Territory',
'Australia:New South Wales',
'Australia:Queensland',
'Australia:Victoria',
'Australia:Western Australia',
'Canada:Alberta',
'Canada:British Columbia',
'Canada:Diamond Princess',
'Canada:Manitoba',
'Canada:Newfoundland and Labrador',
'Canada:Ontario',
'Canada:Quebec',
'Canada:Saskatchewan',
'China:Anhui',
'China:Beijing',
'China:Chongqing',
'China:Fujian',
'China:Gansu',
'China:Guangdong',
'China:Guangxi',
'China:Guizhou',
'China:Hainan',
'China:Hebei',
'China:Heilongjiang',
'China:Henan',
'China:Hubei',
'China:Hunan',
'China:Inner Mongolia',
'China:Jiangxi',
'China:Jilin',
'China:Liaoning',
'China:Shaanxi',
'China:Shandong',
'China:Shanghai',
'China:Sichuan',
'China:Tianjin',
'China:Xinjiang',
'China:Yunnan',
'China:Zhejiang',
'Diamond Princess:',
'France:Guadeloupe',
'France:Martinique',
'France:St Martin',
'Netherlands:Curacao',
'United Kingdom:Cayman Islands',
'United Kingdom:Channel Islands',
'West Bank and Gaza:'}
We see two big ones, China and Australia. Wikipedia has pages with tables for those as well.
page2 = rq.urlopen(wikipedia + 'List_of_Chinese_administrative_divisions_by_population')
s2 = bs(page2, 'html.parser')
t2 = s2.find_all("table", attrs={"class": "wikitable sortable"})[1]
for row in t2.find_all('tr')[2:-1]:
tds = row.find_all('td')
links = tds[0].find_all('a')
province = links[0].getText()
if len(links) > 2:
province = links[1].getText()
population = int("".join(tds[1].getText().split(",")))
obj = dict()
obj["cy"] = "China:" + province
obj["population"] = population
objlist.append(obj)
populations = pd.DataFrame(objlist)
page3 = rq.urlopen(wikipedia + 'States_and_territories_of_Australia')
s3 = bs(page3, 'html.parser')
t3 = s3.find_all("table", attrs={"class": "wikitable sortable"})[0]
for row in t3.find_all('tr')[1:]:
tds = row.find_all('td')
province = tds[1].find_all('a')[0].getText()
population = tds[5].getText()
population = int("".join(tds[5].getText().split(",")))
obj = dict()
obj["cy"] = "Australia:" + province
obj["population"] = population
objlist.append(obj)
populations = pd.DataFrame(objlist)
Let's look at what's left.
df.setset(populations)
{'Australia:Australian Capital Territory',
'Canada:Alberta',
'Canada:British Columbia',
'Canada:Diamond Princess',
'Canada:Manitoba',
'Canada:Newfoundland and Labrador',
'Canada:Ontario',
'Canada:Quebec',
'Canada:Saskatchewan',
'Diamond Princess:',
'France:Guadeloupe',
'France:Martinique',
'France:St Martin',
'Netherlands:Curacao',
'United Kingdom:Cayman Islands',
'United Kingdom:Channel Islands',
'West Bank and Gaza:'}
We could repeat the same for Canada, and maybe look what we can do about the rest, for now we choosate to leave it at the regions we have collected so far.
Converting comultative mortality to mortality rate.
Now that we have a decent population sizes data frame, it is time to request the subset of our mortality data that matches the population size data frame countries and regions. But there is one extra thing to consider other than fetching the propper subset. The data from our COVID-19 data set are commultative mortalities.
While commultative numbers have value, usually rate-numbers are more relevant to most people. Pandas has a diff method for data frames. To make the visualisation less noicy, we suppy diff with a number for its periods parameter. That means we get a value for the new deaths in the last seven days from day seven till the last day of the data set.
ldeadrate2 = df.dropother(populations).diff(7)
dropping United Kingdom:Cayman Islands
dropping Netherlands:Curacao
dropping Canada:Diamond Princess
dropping Diamond Princess:
dropping United Kingdom:Channel Islands
dropping Canada:Manitoba
dropping France:Guadeloupe
dropping Canada:Newfoundland and Labrador
dropping France:St Martin
dropping France:Martinique
dropping Australia:Australian Capital Territory
dropping Canada:Alberta
dropping Canada:British Columbia
dropping Canada:Saskatchewan
dropping West Bank and Gaza:
dropping Canada:Ontario
dropping Canada:Quebec
Normalizing
Now that we have our motrality-rate data and our population data, it is time to normalize our data towards units used in a graph that we would like to compare our data to: The Spanish-flu mortality graph for England.
We are going to transform our mortality rate data from population-wide deaths per week to a yearly mortality rate standarized rate per thousand people. This means we need to multipy our data points by 52 weeks times a thousand people and then divide that nmber by the countries or regions total population size. Let's see what our data looks like then.
for cy in ldeadrate2.columns:
if cy in list(populations["cy"]):
pop = populations[populations["cy"] == cy]["population"].iloc[0]
ldeadrate2[cy] *= 52000/pop
ldr = ldeadrate2[cy]
else:
del ldeadrate2[cy]
ldeadrate2.plot(figsize=(10,10), legend=False)
<matplotlib.axes._subplots.AxesSubplot at 0x7f22599ba5f8>
Let see if we can add some labels to the data to see what regions we are looking at that rice above the rest.
m = ldeadrate2.max(axis=0)
m[m > 0.3].sort_values(ascending=False)
San Marino: 23.232263
Spain: 5.967253
Italy: 4.759492
Andorra: 4.694170
Netherlands: 1.939906
Belgium: 1.917661
France: 1.677619
Switzerland: 1.447380
Luxembourg: 1.185872
China:Hubei 0.846696
United Kingdom: 0.839852
Sweden: 0.608896
Iran: 0.599720
Portugal: 0.592024
Denmark: 0.526898
Austria: 0.508166
Ireland: 0.507162
US: 0.382018
Cyprus: 0.356205
Germany: 0.326449
dtype: float64
topmortality = ldeadrate2[["Andorra:",
"Austria:",
"Belgium:",
"China:Hubei",
"Denmark:",
"France:",
"Iran:",
"Ireland:",
"Italy:",
"Luxembourg:",
"Netherlands:",
"Portugal:",
"San Marino:",
"Spain:",
"Sweden:",
"Switzerland:",
"US:",
"United Kingdom:"]]
topmortality.plot(figsize=(10,10))
<matplotlib.axes._subplots.AxesSubplot at 0x7f22599adb00>
The mortality rate in San Marino is massive. At its peak about twice the average all cause mortality rate for europe. I won't speculate on the meaning or potential spuriousness of the San Marino data, or the dact that unlike most other curves, San Marino seems to be past its peak already.
Before we look at the runners up, let us change the scale a little bit so the x-axis spans a whole year and the y axis goes up to 30.
topmortality.plot(figsize=(10,5),xlim=(0,365), ylim=(0,30))
<matplotlib.axes._subplots.AxesSubplot at 0x7f2252d4e860>
So why bother with whe x and the y axis, well, doing so allows us to compare the current curves with a curve from the infamous Spanish flu.
Looking at it at this level, I don't think this data tells us COVID-19 cant end up as bad or worse as the Spanish Flu. But then, the data also doesn't tell us this whole thing couldn't be over and done with in weeks after a single wave.
But lets look now at a few of the runners up in our data.
topmortality = ldeadrate2[["Andorra:",
"Belgium:",
"China:Hubei",
"France:",
"Italy:",
"Luxembourg:",
"Netherlands:",
"Spain:",
"Switzerland:",]]
topmortality.plot(figsize=(10,10))
<matplotlib.axes._subplots.AxesSubplot at 0x7f2252893470>
The green line is the China:Hubei region, apart from San Marino the only curve in our data past a (first) peak. Very little signs of any of the other curves even approaching the top of the (first) wave. Each of these curves might either soon move towards a peak, or might continue up to San Marino levels or second wave Spanish flu levels. The mortality data on its own doesn't reveal much of the final peak of any of these curves yet. Realize that the current mortality rate for Andora, Spain and Italy has already surpassed that of cardiovascular discease right now.
Lets look at the curves a bit further down.
topmortality = ldeadrate2[["Austria:",
"China:Hubei",
"Denmark:",
"Iran:",
"Ireland:",
"Portugal:",
"Sweden:",
"US:",
"United Kingdom:"]]
topmortality.plot(figsize=(10,10))
<matplotlib.axes._subplots.AxesSubplot at 0x7f22599ba6a0>
Most of the curves look as gloomy as our previous bunch, just a few days behind, but here we do see two curves that hint at slowdown or even that the (first) wave might approach or be at its peak. Iran and Denmark, with Denmark being more likely to be spurious.
Conclusions
I hope this posts shows a number of things, both about the use of pandas, and the COVID-19 mortality data. I hope this data shows that sometimes all that is needed to make data more accessible is data preprocessing and visualization, and that Pandas and numpy are powerfull tools even when you aren't going to be doing any kind of statistical analysis on the data.
I also hope the resulting visualization of COVID-19 data shows that today, just by looking at the graphs, that in most countries and regions, we are in a section of the (first) mortality wave where there is a high amount of uncertainty of both the eventual height end width of the per-country waves.