Effect size and why it matters.

In this notebook I analyze diabetes data downloaded from NCD RisC. The dataset contains age standardized diabetes prevalence percentages for every country of the world and spans the years from 1980 to 2014.

Here I demonstrate that statistical hypothesis testing has limitations that unless taken into consideration can lead to erroneous conlusions even at a compellingly low p value of a hypothesis test.

The notebook and part of the analysis was inspired by the work of Gail M. Sullivan and Richard Feinn published in the Journal of Medical Education in September 2012 under the title "Using effect size - or why the P values is not enough" (and references therein).

# data libs
import numpy as np
import pandas as pd
import scipy.stats as sts

#vizualizations
import plotly.express as px
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import seaborn as sns
from plotly.subplots import make_subplots

# machine learning / forecasting

# catch warnings for better readability
import warnings
warnings.filterwarnings('ignore')

# read in data
data = pd.read_csv('NCD_RisC_Lancet_2016_DM_age_standardised_countries.csv')

data.head()

	Country/Region/World	ISO	Sex	Year	Age-standardised diabetes prevalence	Lower 95% uncertainty interval	Upper 95% uncertainty interval
0	Afghanistan	AFG	Men	1980	0.044712	0.015339	0.094918
1	Afghanistan	AFG	Men	1981	0.046114	0.016883	0.093777
2	Afghanistan	AFG	Men	1982	0.047601	0.018745	0.094018
3	Afghanistan	AFG	Men	1983	0.049173	0.020375	0.093950
4	Afghanistan	AFG	Men	1984	0.050834	0.022269	0.093679

Some feature names are a bit long and cumbersome to work with. Let's rename them to something easier to work with. Also, note that the prevalence is reported in floating point number ranging from 0 to 1. I will turn those to percentages.

data.rename(columns = {'Country/Region/World' : 'Country',
                    'Age-standardised diabetes prevalence':'AgStdPrev',
                    'Lower 95% uncertainty interval':'ci95low',
                    'Upper 95% uncertainty interval':'ci95high'}, inplace=True)

#turn the values into percentages and round them
for c in data.columns[4:]:
    data[c] = data[c]*100

#round to 3rd decimal place
data = data.round(3)

data.head()

	Country	ISO	Sex	Year	AgStdPrev	ci95low	ci95high
0	Afghanistan	AFG	Men	1980	4.471	1.534	9.492
1	Afghanistan	AFG	Men	1981	4.611	1.688	9.378
2	Afghanistan	AFG	Men	1982	4.760	1.875	9.402
3	Afghanistan	AFG	Men	1983	4.917	2.037	9.395
4	Afghanistan	AFG	Men	1984	5.083	2.227	9.368

It might be useful at some point to get the entire confidence interval. Therefore, let's define its' range.

data['ci95'] = data['ci95high'] - data['ci95low']

The data should be clean. But just in case, let's check.

data.isna().sum()

Country      0
ISO          0
Sex          0
Year         0
AgStdPrev    0
ci95low      0
ci95high     0
ci95         0
dtype: int64

We don't have any indication about the trend of the data over time. So the reasonable thing to do is to get the distribution of prevalence values. We do so for men and women separately and compare.

plt.figure(figsize = (12, 7))
sns.histplot(data = data, x = 'AgStdPrev', hue='Sex')

<AxesSubplot:xlabel='AgStdPrev', ylabel='Count'>

The distibutions are reminiscent of Gaussians with a slight right skew. But we'd like to confirm with rigorous statistical testing. The following function performs 4 such tests and summarizes results in the output table. We apply the function on data for men and women.

def normality_tests(x, alpha=0.05):
    """
    Apply multiple tests of normality on vector x at significance alpha
    x: array_like
    alpha: float
    
    Returns
    res_df: Pandas dataframe with test statistic, corresponding p value, 
        and pass/fail of null hypothesis
    """
    # define the tests
    s_shap, p_shap = sts.shapiro(x)
    s_dag, p_dag = sts.normaltest(x)
    s_jb, p_jb = sts.jarque_bera(x)
    s_k, p_k = sts.kstest(x, 'norm')
    
    # define the name of the tests
    test_names = ["Shapiro-Wilk", "D'Agostino", "Jarque-Bera", "Kolmogorov-Smirnov"]
    # create lists of the extracted statistics and p_values
    statistics = [s_shap, s_dag, s_jb, s_k]
    p_values = [p_shap, p_dag, p_jb, p_k]
    # create the dataframe of results
    res_df = pd.DataFrame(dict(test = test_names, statistic = statistics, p_val = p_values))
    res_df['Ho'] = res_df['p_val'].map(lambda y: 'Reject' if y < alpha else 'Accept')
    return res_df

normality_tests(data[data['Sex'] == 'Men']['AgStdPrev'])

	test	statistic	p_val	Ho
0	Shapiro-Wilk	0.762138	0.0	Reject
1	D'Agostino	3948.194457	0.0	Reject
2	Jarque-Bera	33643.630432	0.0	Reject
3	Kolmogorov-Smirnov	0.980218	0.0	Reject

normality_tests(data[data['Sex'] == 'Women']['AgStdPrev'])

	test	statistic	p_val	Ho
0	Shapiro-Wilk	0.798992	0.0	Reject
1	D'Agostino	3268.323622	0.0	Reject
2	Jarque-Bera	18263.775251	0.0	Reject
3	Kolmogorov-Smirnov	0.986663	0.0	Reject

All tests consistently reject the null hypothesis, which is that the distribution is a Gaussian. One final check is QQ plots. We use the probplot function from the scipy library to produce those.

sts.probplot(data[data['Sex'] == 'Men']['AgStdPrev'], dist = 'norm', plot = plt);

sts.probplot(data[data['Sex'] == 'Women']['AgStdPrev'], dist = 'norm', plot = plt);

The conclusion is that these are not Gaussian distributions.

The extended overlap seems to leave no doubt that drawing a sample would be almost impossible to distinguish if it came from a woman or a man. Nevertheless, we go ahead and run a hypothesis test. We will use a non-parametric test suitable for independent samples, the Mann-Whitney U test. Mann-Whitney does not depend on the distribution of values (i.e. if it is normal or not etc). The null hypothesis is that the median of the distributions are equivalent and the alternative is that one is higher than the other.

Note that the two samples are quite large with a length of 7000 each!

stat_mwu, pval_mwu = sts.mannwhitneyu(data[data['Sex'] == 'Women']['AgStdPrev'], 
                 data[data['Sex'] == 'Men']['AgStdPrev'],
                use_continuity=True)
print(f"The p-value of the Mann-Whitney non parametric test is {np.round(pval_mwu, 3)}")

The p-value of the Mann-Whitney non parametric test is 0.036

The hypothesis testing rejects the Ho! Our intuition, based on the distribution plots, was telling us the exact opposite!

Note that since the sample size is quite large we could run a 2 (independent) sample ttest. That would only result in a lower p value.

This is were effect size comes at play. Effect size is a statistical measure of the magnitude of the difference of the central tendencies of the distributions between 2 groups.

To begin we combine the data of the two arrays in a new dataframe, named comparison_df, and then calculate means and standard deviations for men (M) and women (W).

comparison_df = pd.DataFrame(data = {"M":np.array(data[data['Sex'] == 'Men']['AgStdPrev']),
                                    "W":np.array(data[data['Sex'] == 'Women']['AgStdPrev'])})

means = comparison_df.mean()
stds = comparison_df.std()

Effect Size Calculations

We are now ready to calculate effect size. We concentrate on the standardized mean difference between two groups. We use two metrics Cohen's d and Hedges' g defined as following.

$Cohen's d = \frac{M_{1}-M_{2}}{SD_{pooled}}$, $SD_{pooled} = \sqrt{\frac{SD_{1}^{2} + SD_{2}^{2}}{2}}$

$Hedges' g = \frac{M_{1}-M_{2}}{SD_{pooled}^{*}}$, $SD_{pooled}^{*} = \sqrt{\frac{(n_{1}-1)SD_{1}^{2} + (n_{2}-1)SD_{2}^{2}}{n_{1} + n_{2} - 2}}$

$M$ is the mean, $SD$ refers to the standard deviation of the group, and $n$ is the population size.

We will also use Pearon's r coefficent as a measure of association between the two groups. For all metrics the following table is a guide on the result.

Value	Effect size
0.2	low
0.5	medium
0.8	strong

# calcuate pooled standard deviations
s_pooled = np.sqrt(stds.loc['M']**2 + stds.loc['W']**2)*0.5

l = len(comparison_df) - 1  #the lengths for men and women are equal therefore use on value
s_pooled_star = np.sqrt((l*stds.loc['M']**2 + l*stds.loc['W']**2)/((2*l) - 2))

cohend = np.round((means.loc['W'] - means.loc['M'])/s_pooled, 3)
hedgeg = np.around((means.loc['W'] - means.loc['M'])/s_pooled_star, 3)
pr, pval = sts.pearsonr(comparison_df['W'], comparison_df['M'])

effects_df = pd.DataFrame(dict(Statistical_index = ["Cohen's d", "Hedge's g", "Pearson r assoc."],
                              index_value = [cohend, hedgeg, pr]))
effects_df.round(3)

	Statistical_index	index_value
0	Cohen's d	0.075
1	Hedge's g	0.053
2	Pearson r assoc.	0.951

The conclusion is that the effect size analysis reveals a negligible substansive significance between the distributions of the two groups. This was expected from the histogram plots but was not supported by the hypothesis testing.

Therefore, there are 2 key take aways:

1. Statistical plots, such as histograms, are a very useful way of identifying differences in distributions.
2. Although a P value can be misleading, the effect size is not.

Map prevalence in time

Since there is no statistical significance of the values between men and women I define a new mean of the age standardized prevalence and use plotly's choropleth function to give a visual of the time dependence of that mean. The visual was saved as a separate html file.

From the time dependent map we see that diabetes had a significant increase in the North African countries, the Arab peninsula.