COVID-19/All-cause deaths/Scripts
The following Python scripts were used to produce charts. They are small and very uncomplicated, very straightforward.
plotHmd.py[edit | edit source]
import sys, csv, datetime
fileName = sys.argv[1] # e.g. stmf.csv
countryCode = sys.argv[2] # e.g. USA
figureFieldName = sys.argv[3] # e.g. D0_14 (deaths for 0-14y), DTotal (total deaths)
smoothingCount = int(sys.argv[4]) if len(sys.argv) >= 4 + 1 else 1 # e.g. 3
data = []
with open(fileName) as file1:
file1.readline() # Two intro lines without data
file1.readline()
for line in csv.DictReader(file1):
if line and line["CountryCode"] == countryCode and line["Sex"] == "b":
deaths = int(float(line[figureFieldName]))
date = line["Year"] + " " + line["Week"] + " 0"
date = datetime.datetime.strptime(date, "%Y %U %w")
data.append( (date, deaths) )
data.pop() # Drop last two weeks for too big registration delay effect
data.pop()
def movingAverage(values, itemCount):
average = 0
averages = []
idx = -1
for val in values:
idx += 1
average += val - values[idx - itemCount] if idx >= itemCount else val
if idx >= itemCount - 1:
averages.append(average / float(itemCount))
else:
averages.append(None)
return averages
values = [v for k, v in data]
if smoothingCount > 1:
values = movingAverage(values, smoothingCount)
maxValue = max([float(v) for v in values if v is not None])
valueFormatString = "%.1f" if maxValue < 100 else "%.0f"
sys.stdout.write("|x = " + ", ".join([k.strftime("%Y-%m-%d") for k, v in data]) + "\n")
sys.stdout.write("|y = " + ", ".join([valueFormatString % v if v is not None else "" for v in values]) + "\n")
Usage:
- plotHmd.py stmf.csv USA DTotal
- plotHmd.py stmf.csv USA D0_14
- plotHmd.py stmf.csv USA D0_14 3
To obtain stmf.csv, go to mpidr.shinyapps.io/stmortality and at the bottom of the left pane click on the icon to the right of "CSV".
Country codes: AUS (Australia), AUT (Austria), BEL (Belgium), BGR (Bulgaria), CAN (Canada), CHE (Switzerland), CHL (Chile), CZE (Czechia), DEUTNP (Germany), DNK (Denmark), ESP (Spain), EST (Estonia), FIN (Finland), FRATNP (France), GBRTENW (England and Wales), GBR_NIR (Northern Ireland), GBR_SCO (Scotland), GRC (Greece), HRV (Croatia), HUN (Hungary), ISL (Iceland), ISR (Israel), ITA (Italy), KOR (Republic of Korea), LTU (Lithuania), LUX (Luxembourg), LVA (Latvia), NLD (Netherlands), NOR (Norway), NZL_NP (New Zealand), POL (Poland), PRT (Portugal), RUS (Russia), SVK (Slovakia), SVN (Slovenia), SWE (Sweden), TWN (Taiwan), USA (U.S.). See also W:ISO 3166-1 alpha-3; note that some codes in stmf.csv have custom suffix such as the custom codes for parts of the United Kingdom (GBRTENW, etc.) or FRATNP instead of FRA. See also HMD-countries-codes.pdf, mortality.org; however, this list contains multiple codes absent from stmf.csv, such as FRACNP and UKR.
Field codes: D0_14 (deaths for 0-14y), D15_64 (deaths for 15-64y), D65_74 (deaths for 65-74y), D75_84 (deaths for 75-84y), D85p (deaths for 85+y), DTotal (deaths total).
The script drops last two data points to prevent the worst effect of registration delay; for some countries, the last two weeks were obviously very badly affected by registration delay.
plotHmdPerYear.py[edit | edit source]
import sys, csv, datetime
fileName = sys.argv[1] # e.g. stmf.csv
countryCode = sys.argv[2] # e.g. USA
figureFieldName = sys.argv[3] # e.g. D0_14 (deaths for 0-14y), DTotal (total deaths)
data = []
file1 = open(fileName)
file1.readline() # Two intro lines without data
file1.readline()
for line in csv.DictReader(file1):
if line and line["CountryCode"] == countryCode and line["Sex"] == "b":
deaths = int(float(line[figureFieldName]))
data.append( (int(line["Year"]), int(line["Week"]), deaths) )
data.pop() # Drop last two weeks for too big registration delay effect
data.pop()
years = sorted(list({year for year, week, deaths in data}))
maxYear = max(years)
maxYearWeeks = [week for year, week, deaths in data if year == maxYear]
maxWeek = maxYearWeeks[-1]
deathsUpToMaxWeek = []
for year in years:
deathsUpToMaxWeek1 = 0
for year1, week, deaths in data:
if year1 == year and week <= maxWeek:
deathsUpToMaxWeek1 += deaths
deathsUpToMaxWeek.append(deathsUpToMaxWeek1)
yearsOut = ", ".join([str(year) for year in years])
deathsOut = ", ".join([str(deaths) for deaths in deathsUpToMaxWeek])
sys.stdout.write("Last week in %i: %i\n" % (maxYear, maxWeek))
sys.stdout.write("|x = " + yearsOut + "\n")
sys.stdout.write("|y = " + deathsOut + "\n")
Usage: similar to plotHmd.py.
plotHmdPerSeason.py[edit | edit source]
import sys, csv, datetime
# Plot deaths per season: week 40 of year before to week x of the year
fileName = sys.argv[1] # e.g. stmf.csv
countryCode = sys.argv[2] # e.g. USA
figureFieldName = sys.argv[3] # e.g. D0_14 (deaths for 0-14y), DTotal (total deaths)
seasonStartWeek = 40
data = []
file1 = open(fileName)
file1.readline() # Two intro lines without data
file1.readline()
for line in csv.DictReader(file1):
if line and line["CountryCode"] == countryCode and line["Sex"] == "b":
deaths = int(float(line[figureFieldName]))
data.append( (int(line["Year"]), int(line["Week"]), deaths) )
data.pop() # Drop last two weeks for too big registration delay effect
data.pop()
years = sorted(list({year for year, week, deaths in data}))
maxYear = max(years)
maxYearWeeks = [week for year, week, deaths in data if year == maxYear]
maxWeek = maxYearWeeks[-1]
if maxWeek >= seasonStartWeek:
maxWeek = seasonStartWeek - 1
deathsInSeason = []
for year in years[1:]:
deathsInSeason1 = 0
for year1, week, deaths in data:
if year1 == year and week <= maxWeek:
deathsInSeason1 += deaths
if year1 == (year - 1) and week >= seasonStartWeek:
deathsInSeason1 += deaths
deathsInSeason.append(deathsInSeason1)
yearsOut = ", ".join([str(year) for year in years[1:]])
deathsOut = ", ".join([str(deaths) for deaths in deathsInSeason])
write = sys.stdout.write
write("Last week in %i: %i\n" % (maxYear, maxWeek))
write("All-cause deaths in weeks %i+ year before and weeks 1-%i of the year, year by year:\n" % (seasonStartWeek, maxWeek))
write("|x = " + yearsOut + "\n")
write("|y = " + deathsOut + "\n")
Usage: similar to plotHmd.py.
plotHmdExcessDeathPercPerYear.py[edit | edit source]
import sys, csv
fileName = sys.argv[1] # e.g. stmf.csv
countryCode = sys.argv[2] # e.g. USA
figureFieldName = sys.argv[3] # e.g. D0_14 (deaths for 0-14y), DTotal (total deaths)
baseYearRangeLen = int(sys.argv[4]) # e.g. 5
data = []
file1 = open(fileName)
file1.readline() # Two intro lines without data
file1.readline()
for line in csv.DictReader(file1):
if line and line["CountryCode"] == countryCode and line["Sex"] == "b":
deaths = int(float(line[figureFieldName]))
data.append( (int(line["Year"]), int(line["Week"]), deaths) )
data.pop() # Drop last two weeks for too big registration delay effect
data.pop()
years = sorted(list({year for year, week, deaths in data}))
years = years[:-1] # Drop the last year since it is usually incomplete
minYear = min(years)
maxYear = max(years)
calculableYears = [year for year in years if year >= minYear + baseYearRangeLen]
deathsPerYear = {}
for year, week, deaths in data:
if year <= maxYear:
if not year in deathsPerYear:
deathsPerYear[year] = 0
deathsPerYear[year] += deaths
excessPercentPerYear = {}
for year in calculableYears:
baseYearRangeMin = None
for yearOffset in range(baseYearRangeLen, 0, -1):
if baseYearRangeMin is None:
baseYearRangeMin = deathsPerYear[year - yearOffset]
else:
baseYearRangeMin = min(baseYearRangeMin, deathsPerYear[year - yearOffset])
excessPercentPerYear[year] = deathsPerYear[year] / float(baseYearRangeMin) - 1
yearsOut = ", ".join([str(year) for year in calculableYears])
excessPercOut = ", ".join([("%.3f" % excessPercentPerYear[year]) for year in calculableYears])
sys.stdout.write("|x = " + yearsOut + "\n")
sys.stdout.write("|y = " + excessPercOut + "\n")
Usage (similar to plotHmd.py).
- plotHmdExcessDeathPercPerYear.py stmf.csv BGR DTotal 5
Rationale for using minimum of a range of years as the baseline: visual inspection of the weekly charts shows that significant year-specific variation occurs in the upward direction but not in the downward direction. Thus, the year for which as little year-specific variation occurred as possible is taken to be the baseline; if we took the average, we would include previous year-specific upward variation into the baseline.
plotUsCdc.py[edit | edit source]
import sys, csv, datetime
fileName = sys.argv[1] # e.g. "Excess_Deaths_Associated_with_COVID-19.csv"
jurisdiction = sys.argv[2] # e.g. "New York City" or "Alabama"
data = []
file1 = open(fileName)
firstThree = file1.read(3) # Drop BOM if present
if firstThree == "Wee": # No BOM
file1.seek(0)
for line in csv.DictReader(file1):
if line["State"] == jurisdiction and line["Type"] == "Unweighted":
date = datetime.datetime.strptime(line["Week Ending Date"], "%Y-%m-%d")
deaths = line["Observed Number"].replace(",", "")
data.append( (date, deaths) )
data.sort(key=lambda x: x[0])
datesOut = ", ".join([k.strftime("%Y-%m-%d") for k, v in data])
deathsOut = ", ".join([v for k, v in data])
sys.stdout.write("|x =" + datesOut + "\n")
sys.stdout.write("|y =" + deathsOut + "\n")
Usage:
- plotUsCdc.py Excess_Deaths_Associated_with_COVID-19.csv "New York City"
To obtain Excess_Deaths_Associated_with_COVID-19.csv:
- 1) Visit CDC[1].
- 2) In section "Download Data:", click on "National and State Estimates of Excess Deaths".
- 3) Save file "Excess_Deaths_Associated_with_COVID-19.csv", which contains data for all jurisdictions.
plotWmd.py[edit | edit source]
import sys, csv, datetime, argparse
parser = argparse.ArgumentParser(description='Plot Wmd for wiki.')
parser.add_argument("fileName") # e.g. world_mortality.csv from akarlinsky, hithub.com
parser.add_argument("countryCode") # e.g. USA
parser.add_argument("smoothingCount", nargs="?", type=int, default=1) # e.g. 3
parser.add_argument("--dlti", dest="dropLastTwoItems", action="store_true",
help="Drop last two items in the data to prevet registration delay")
parser.add_argument("--ep", dest="excessDeathRangeLen", type=int,
help="Excess death percentage: length of ref. year range")
args = parser.parse_args()
def fillDataFromWmdFile(fileName, countryCode, dropLastTwoItems):
data = []
with open(fileName) as file1:
for line in csv.DictReader(file1):
if line["iso3c"] == countryCode:
if line["time_unit"] == "weekly":
week = line["time"]
date = line["year"] + " " + week + " 0"
date = datetime.datetime.strptime(date, "%Y %U %w")
elif line["time_unit"] == "monthly":
month = line["time"]
date = line["year"] + " " + month + " 1"
date = datetime.datetime.strptime(date, "%Y %m %d")
else:
sys.stderr.write("Unexpected time unit. Aborting.")
sys.exit(1)
deaths = int(float(line["deaths"]))
data.append( (date, deaths) )
if dropLastTwoItems:
data.pop() # Drop last two weeks for too big registration delay effect
data.pop()
return data
def movingAverage(values, itemCount):
average = 0
averages = []
idx = -1
for val in values:
idx += 1
average += val - values[idx - itemCount] if idx >= itemCount else val
if idx >= itemCount - 1:
averages.append(average / float(itemCount))
else:
averages.append(None)
return averages
def outputWeeklyOrMonthlyTimeSeriesForWikiChart(data, smoothingCount):
values = [v for k, v in data]
if smoothingCount > 1:
values = movingAverage(values, smoothingCount)
maxValue = max([float(v) for v in values if v is not None])
valueFormatString = "%.1f" if maxValue < 100 else "%.0f"
sys.stdout.write("|x = " + ", ".join([k.strftime("%Y-%m-%d") for k, v in data]) + "\n")
sys.stdout.write("|y = " + ", ".join([valueFormatString % v if v is not None else "" for v in values]) + "\n")
def outputYearlyExcessDeathPercentForWikiChart(data, baseYearRangeLen): # data is a list of (date, death) pairs
years = sorted(list({date.year for date, deaths in data}))
years = years[:-1] # Drop last year as incomplete
minYear = min(years)
maxYear = max(years)
calculableYears = [year for year in years if year >= minYear + baseYearRangeLen]
deathsPerYear = {}
for date, deaths in data:
year = date.year
if year <= maxYear:
if not year in deathsPerYear:
deathsPerYear[year] = 0
deathsPerYear[year] += deaths
excessPercentPerYear = {}
for year in calculableYears:
baseYearRangeMin = None
for yearOffset in range(baseYearRangeLen, 0, -1):
if baseYearRangeMin is None:
baseYearRangeMin = deathsPerYear[year - yearOffset]
else:
baseYearRangeMin = min(baseYearRangeMin, deathsPerYear[year - yearOffset])
excessPercentPerYear[year] = deathsPerYear[year] / float(baseYearRangeMin) - 1
yearsOut = ", ".join([str(year) for year in calculableYears])
excessPercOut = ", ".join([("%.3f" % excessPercentPerYear[year]) for year in calculableYears])
sys.stdout.write("|x = " + yearsOut + "\n")
sys.stdout.write("|y = " + excessPercOut + "\n")
data = fillDataFromWmdFile(args.fileName, args.countryCode, args.dropLastTwoItems)
if args.excessDeathRangeLen is None:
outputWeeklyOrMonthlyTimeSeriesForWikiChart(data, args.smoothingCount)
else:
outputYearlyExcessDeathPercentForWikiChart(data, args.excessDeathRangeLen)
Usage:
- plotWmd.py world_mortality.csv PER 3
- plotWmd.py world_mortality.csv PER --ep 3
Moving average via awk[edit | edit source]
You can calculate the 7-day moving average using awk on Windows:
- echo 1, 0, 4, 5, 18, 15, 28, 26, 64, 77, 101 | awk -F, -vn=7 "{for(i=1;i<=NF; i++) {s+=i>n?$i-$(i-n):$i; if(i>=n){printf \"%.0f, \", s/n}else{printf \", \"}}}"
You can put the result into clipboard:
- echo 1, 0, 4, 5, 18, 15, 28, 26, 64, 77, 101 | awk -F, -vn=7 "{for(i=1;i<=NF; i++) {s+=i>n?$i-$(i-n):$i; if(i>=n){printf \"%.0f, \", s/n}else{printf \", \"}}}" | clip
You can do the calculation on Linux:
- echo 1, 0, 4, 5, 18, 15, 28, 26, 64, 77, 101 | awk -F, -vn=7 '{for(i=1;i<=NF; i++) {s+=i>n?$i-$(i-n):$i; if(i>=n){printf "%.0f, ", s/n}else{printf ", "}}}'
If you are on Linux or a modern Mac, you already have awk. For Windows, you can install awk from ezwinports or GnuWin32 project.