2
\$\begingroup\$

This script is a further development of Finding crude death rate and age-standardized death rates from data. In reflection to the kind reviews I created a module for the data; saved the data in dictionaries, with the age groups signified as keys; improved variable names; added function definitions and typing; a description at the beginning of the script; used _ separators for 100_000, and 100_000 consistently instead of (10**5); made output more descriptive. In summary, I tried making the code more accessible to readers.

I am unsure about the correctness of the calculations, whether the values are approximately realistic, and why there is a difference of values.

copd_data.py

"""
Data sources:
    - [UN World Population Prospects (2022) — Population by Five-year Age Groups - Both Sexes](https://population.un.org/wpp/Download/Files/1_Indicators%20(Standard)/EXCEL_FILES/2_Population/WPP2022_POP_F02_1_POPULATION_5-YEAR_AGE_GROUPS_BOTH_SEXES.xlsx)
    - [WHO Standard Population — Table 1 in 'Ahmad OB, Boschi-Pinto C, Lopez AD, Murray CJ, Lozano R, Inoue M (2001). Age standardization of rates: a new WHO standard.'](https://www.researchgate.net/publication/238744905_Ahmad_OB_Boschi-Pinto_C_Lopez_AD_et_al_2000_Age_Standardization_of_Rates_A_New_WHO_Standard_GPE_Discussion_Paper_Series_No_31_World_Health)
    - [Table of age-specific death rates of COPD](https://owid.notion.site/Data-analysis-exercise-Our-World-in-Data-Junior-Data-Scientist-application-ab287a3c07264b4d91aadc436021b8c0)
"""

# COPD deaths in US, in 2019, per 100_000
us_copd_death_rates_2019 = {
    "0-4": 0.04,
    "5-9": 0.02,
    "10-14": 0.02,
    "15-19": 0.02,
    "20-24": 0.06,
    "25-29": 0.11,
    "30-34": 0.29,
    "35-39": 0.56,
    "40-44": 1.42,
    "45-49": 4.00,
    "50-54": 14.13,
    "55-59": 37.22,
    "60-64": 66.48,
    "65-69": 108.66,
    "70-74": 213.10,
    "75-79": 333.06,
    "80-84": 491.10,
    "85+": 894.45,
}

# COPD deaths in Uganda, in 2019, per 100_000
uganda_copd_death_rates_2019 = {
    "0-4": 0.40,
    "5-9": 0.17,
    "10-14": 0.07,
    "15-19": 0.23,
    "20-24": 0.38,
    "25-29": 0.40,
    "30-34": 0.75,
    "35-39": 1.11,
    "40-44": 2.04,
    "45-49": 5.51,
    "50-54": 13.26,
    "55-59": 33.25,
    "60-64": 69.62,
    "65-69": 120.78,
    "70-74": 229.88,
    "75-79": 341.06,
    "80-84": 529.31,
    "85+": 710.40,
}

# US 2019 population by age group in thousands
us_population_by_age_group_2019 = {
    "0-4": 19849,
    "5-9": 20697,
    "10-14": 22092,
    "15-19": 21895,
    "20-24": 21872,
    "25-29": 23407,
    "30-34": 22842,
    "35-39": 22297,
    "40-44": 20695,
    "45-49": 21244,
    "50-54": 21346,
    "55-59": 22348,
    "60-64": 20941,
    "65-69": 17501,
    "70-74": 13689,
    "75-79": 9273,
    "80-84": 6119,
    "85+": 6214,
}

# Uganda 2019 population by age group in thousands
uganda_population_by_age_group_2019 = {
    "0-4": 7329,
    "5-9": 6614,
    "10-14": 5899,
    "15-19": 5151,
    "20-24": 4348,
    "25-29": 3500,
    "30-34": 2619,
    "35-39": 1903,
    "40-44": 1504,
    "45-49": 1235,
    "50-54": 953,
    "55-59": 687,
    "60-64": 500,
    "65-69": 353,
    "70-74": 197,
    "75-79": 93,
    "80-84": 44,
    "85+": 20,
}

# World Health Organisation standard population weights, by age group, in percent
who_standard_weights_by_age_group = {
    "0-4": 8.86,
    "5-9": 8.69,
    "10-14": 8.6,
    "15-19": 8.47,
    "20-24": 8.22,
    "25-29": 7.93,
    "30-34": 7.61,
    "35-39": 7.15,
    "40-44": 6.59,
    "45-49": 6.04,
    "50-54": 5.37,
    "55-59": 4.55,
    "60-64": 3.72,
    "65-69": 2.96,
    "70-74": 2.21,
    "75-79": 1.52,
    "80-84": 0.91,
    "85+": 0.63,
}

copd_death_rates.py

"""
The purpuse of this code is to calculate 2019 crude death rates, and age-standardised death rates for chronic obstructive pulmonary disease (COPD), in the US and Uganda.
"""

import numpy as np
from typing import Dict
from copd_data import *


def calculate_total_deaths(
    population_by_age_group: Dict[str, int], death_rates: Dict[str, float]
) -> float:
    """
    Calculate total deaths of total population, from population by age group (in thousands), and death rates by age groups.

    Input type dictionary, outputs np.ndarray.
    """
    pop_arr = np.array(list(population_by_age_group.values())) * 1000
    death_rates_arr = np.array(list(death_rates.values()))
    return np.sum(pop_arr * death_rates_arr) / (100_000)


def crude_death_rate(
    population_by_age_group: Dict[str, int], total_deaths: np.ndarray
) -> float:
    """
    Calculate crude death rate per 100_000, from given population by age group (in thousands), and total deaths in total population.

    Input population as type dictionary, outputs np.ndarray.
    """

    pop_arr = np.array(list(population_by_age_group.values())) * 1000
    return (total_deaths / np.sum(pop_arr)) * (100_000)


def age_standardised_death_rate(
    population_by_age_group: Dict[str, int],
    standard_weights: Dict[str, int],
    death_rates: Dict[str, int],
) -> float:
    """
    Calculate age-standardised death rate, total standard deaths to a standard population:

    sum(population_by_age_group * weights * death_rates) / sum(population_by_age_group * weights)
    """

    pop_arr = np.array(list(population_by_age_group.values())) * 1000
    weights_arr = np.array(list(standard_weights.values())) / 100
    death_rates_arr = np.array(list(death_rates.values()))
    standard_deaths = np.sum((pop_arr / 100_000) * weights_arr * death_rates_arr)
    total_standard_population = np.sum(pop_arr * weights_arr)
    return (standard_deaths / total_standard_population) * (100_000)


us_total_copd_deaths_2019 = calculate_total_deaths(
    us_population_by_age_group_2019, us_copd_death_rates_2019
)
uganda_total_copd_deaths_2019 = calculate_total_deaths(
    uganda_population_by_age_group_2019, uganda_copd_death_rates_2019
)

us_crude_death_rate = crude_death_rate(
    us_population_by_age_group_2019, us_total_copd_deaths_2019
)
uganda_crude_death_rate = crude_death_rate(
    uganda_population_by_age_group_2019, uganda_total_copd_deaths_2019
)

us_age_standardised_copd_death_rate_2019 = age_standardised_death_rate(
    us_population_by_age_group_2019,
    who_standard_weights_by_age_group,
    us_copd_death_rates_2019,
)
uganda_age_standardised_copd_death_rate_2019 = age_standardised_death_rate(
    uganda_population_by_age_group_2019,
    who_standard_weights_by_age_group,
    uganda_copd_death_rates_2019,
)

print(
    f"US 2019 COPD crude death rate, (total deaths / total population) per 100_000: {us_crude_death_rate:.1f}"
)
print(
    f"Uganda 2019 COPD crude death rate, (total deaths / total population) per 100_000: {uganda_crude_death_rate:.1f}"
)
print(
    f"US 2019 COPD age standardized death rate, (total standard deaths / total standard population), per 100_000: {us_age_standardised_copd_death_rate_2019:.1f}"
)
print(
    f"Uganda 2019 COPD age standardized death rate, (total standard deaths / total standard population), per 100_000: {uganda_age_standardised_copd_death_rate_2019:.1f}"
)

the output:

US 2019 COPD crude death rate, (total deaths / total population) per 100_000: 57.2
Uganda 2019 COPD crude death rate, (total deaths / total population) per 100_000: 5.8
US 2019 COPD age standardized death rate, (total standard deaths / total standard population), per 100_000: 16.5
Uganda 2019 COPD age standardized death rate, (total standard deaths / total standard population), per 100_000: 2.2
\$\endgroup\$
0

3 Answers 3

5
\$\begingroup\$

Use data frames

For me, global variables with very long names are a code smell. In this case, uganda_age_standardised_copd_death_rate_2019 and us_age_standardised_copd_death_rate_2019 are not unmanageable. But what happens when you have 100 countries, 20 years of data and five derived metrics?

I am inclined to do this using pandas.DataFrame.pipe(). This allows you to write code like this:

def calculate_death_metrics(df, groups = ["org", "year"]):
    return (
        df
        .groupby(groups)
        .pipe(calculate_total_deaths)
        .merge(df, how = "right", on = groups)
        .groupby(groups)
        .pipe(crude_death_rate)
        .merge(df, how = "right", on = groups)
        .groupby(groups)
        .pipe(age_standardised_death_rate)
    )

Which will produce output that scales better:

calculate_death_metrics(df)
      org  year  total_deaths  crude_death_rate  age_standardised_death_rate
0      US  2019   191354.2049         57.236669                    16.465936
1  Uganda  2019     2500.8292          5.822788                     2.249041

Getting your data into a data frame

Here is a function that will do this:

import pandas as pd

def create_df(data_dict: Dict[str, int], org: str, metric: str, year: int) -> pd.DataFrame:
    df_dict = {k: [org, year, v] for k,v in data_dict.items()}
    return (
        pd.DataFrame(
            df_dict
        )
        .transpose()
        .reset_index()
        .rename(
            columns = {0: "org", 1: "year", 2: metric, "index" : "age_group"}
        )
    )

We can then use this with your input to create a data frame called df:

rates_df = pd.concat([
    create_df(us_copd_death_rates_2019, "US", "copd_death_rate", 2019),
    create_df(uganda_copd_death_rates_2019, "Uganda", "copd_death_rate", 2019)
])
pop_df = pd.concat([
    create_df(us_population_by_age_group_2019 , "US", "pop", 2019),
    create_df(uganda_population_by_age_group_2019 , "Uganda", "pop", 2019)
])
weights_df = (
    create_df(who_standard_weights_by_age_group , "WHO", "weights", 2019)
)
df = rates_df.merge(
    pop_df,
    on = ["age_group", "org", "year"]
).merge(weights_df[["age_group", "year", "weights"]],
    on = ["age_group", "year"])

Rewriting your functions

We can also rewrite each of your functions to be used in this workflow. I haven't bothered with type annotations here but they all take a grouped data frame and return an ungrouped data frame.

def calculate_total_deaths(grp, per = 1):
    return pd.DataFrame([
        [name[0], name[1], np.sum(vals["pop"] * vals["copd_death_rate"] / per)]
                     for name, vals in grp
    ], columns = ["org", "year", "total_deaths"])
def crude_death_rate(grp):
     return pd.DataFrame([
        [
             name[0], 
             name[1], 
             total_deaths := vals["total_deaths"].unique()[0],
             (total_deaths / np.sum(vals["pop"]))
         ]
          for name, vals in grp
    ], columns = ["org", "year", "total_deaths", "crude_death_rate"])
def age_standardised_death_rate(grp):
    return pd.DataFrame([
        [
             name[0], 
             name[1],
             vals["total_deaths"].unique()[0] / 100,
             vals["crude_death_rate"].unique()[0],
             np.sum(vals["pop"] * vals["weights"] * vals["copd_death_rate"]) /
             np.sum(vals["pop"] * vals["weights"])
         ]
          for name, vals in grp
    ], columns = ["org", "year", "total_deaths", "crude_death_rate", "age_standardised_death_rate"])

In your example, you have two countries, one year and three metrics, so only six global variables. But hopefully you see the advantage of data frames once you have more combinations.

\$\endgroup\$
4
\$\begingroup\$

As far as correctness is concerned, if your us_copd_death_rates_2019 dictionary really is deaths per 100_000 broken down by age groups and if us_population_by_age_group_2019 really is the population in thousands by age groups, then your function calculate_total_deaths does calculate correctly the total deaths within a precision that is limited by the precision of the input data. Because of this imprecision, I am not sure what the correct maximum number of significant digits you can use in your output. I also have not spent any great amount of time looking at the other calculations. So my comments will be few:

Documenting the Input Data

I see that you have given the input data good descriptions and meaningful names. Great!

Documenting the Functions

You have provided both docstrings and type hints. Again, great! As long as the information provided in each is not redundant, then having both is a plus. But you do have, for example in function crude_death_rate, a docstring that describes the input and output that is redundant with the type hints. Also, for example, the type hint for the total_deaths argument should be float.

Can We Simplify the Calculations?

Let's look at function calculate_total_deaths. Your calculation is:

    ...

    pop_arr = np.array(list(population_by_age_group.values())) * 1000
    death_rates_arr = np.array(list(death_rates.values()))
    return np.sum(pop_arr * death_rates_arr) / (100_000)

If we want to emphasize clarity over efficiency or getting the most precise answer (having the fewest floating point rounding errors), then you might consider instead:

    ...

    pop_arr = np.array(list(population_by_age_group.values())) * 1000
    death_rates_arr = np.array(list(death_rates.values())) / 100_000
    return np.sum(pop_arr * death_rates_arr)

This makes it clearer that the death_rates_arr has converted input that represents death rates per per 100,000 to death rates per person. Then, clearly when you multiply the total number of people (pop_arr) by the death rate per person, you arrive at the total number of deaths. In general, however, the two calculations can produce slightly different results due to the finite precision of floating point. The second, "clearer" calculation introduces more rounding errors than the first,

If you want to optimize the calculation as far as efficiency and reducing rounding errors are concerned, then you can perform constant folding:

    ...

    pop_arr = np.array(list(population_by_age_group.values()))
    death_rates_arr = np.array(list(death_rates.values()))
    return np.sum(pop_arr * death_rates_arr) / 100

But if we were really concerned with efficiency, you could consider eliminating the use of numpy altogether. numpy makes multiplying arrays together very simple. The overhead, however, is that you have to first convert the values of the population_by_age_group dictionary to a list before you can initialize the death_rates_arr numpy array. We could instead do:

from itertools import starmap
from operator import mul

    ...

    return sum(
        starmap(mul, zip(population_by_age_group.values(), death_rates.values()))
    ) / 10

Function crude_death_rate becomes even simpler:

def crude_death_rate(
    population_by_age_group: Dict[str, int], total_deaths: float
) -> float:
    """
    Calculate crude death rate per 100_000 from a given population
    by age group (in thousands), and total deaths in that population.
    """

    return 100 * total_deaths / sum(population_by_age_group.values())

Note that I have updated the docstring to remove redundant information that was already provided by type hinting and to keep the line lengths to a more reasonable size.

These functions that do not use numpy actually run faster than the original functions, but admittedly speed does not seem to be an issue with these calculations. So the question is, "Have we sacrificed clarity by not using numpy?" I don't think so.

\$\endgroup\$
3
\$\begingroup\$

I am unsure about the correctness of the calculations, whether the values are approximately realistic, and why there is a difference of values.

This is a serious problem; anything else you might do to your code is moot if it doesn't work, or even if you're not sure if it works. So somehow you have to figure out how to test/validate your work.

Some formatting stuff:

  • Line length: I'm not a stickler for the the 80-char limit, but it's never nice to be scrolling sideways to read stuff, and line-wrapping is only a little better. Just use more line breaks, mostly in your comments and your output.
  • Since version 3.9, you don't need to from typing import Dict, you can just annotate stuff as
    : dict[t1, t2]. That said, typing is a big subject, and saying something is a dict is actually quite over-specific. The best annotation for everything you're doing here is typing.Mapping.
  • Since you're using type annotations, don't bother with them in the comments.

Engineering stuff:

  • Don't use types if you're not going to typecheck them! Your IDE may have a built-in type-checker for python; otherwise an easy choice is to pip install mypy, and then you can run (on the command line, not inside your code) mypy filename. This flags a few issues in your code; these aren't "errors", they're just things that the type-checker thinks are wrong. You should certainly "fix" them (add or change annotations); convincing the type-checker the code is correct is a first step toward convincing yourself!
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.