Statistical analysis of data on natural selection

I tested whether humerus length in sparrows is related to surviving a severe winter storm. The dataset has humerus measurements for birds that survived and birds that died. The analysis covers point estimates, distribution estimates, parametric model fitting, simulation, confidence intervals, and t-tests.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as st

pd.set_option('display.float_format', '{:g}'.format)
sns.set()

dataset

K = 16
L = 3
M = ((K + L) * 47) % 11 + 1
pd.DataFrame([K, L, M], index=["K", "L", "M"], columns=[""])

introduction

First I load the data and split the humerus lengths into the two groups. Then I describe the data and estimate the mean, variance, and median for each group separately.

data = pd.read_csv("case0201.csv")

S = data.loc[data["Status"] == "Survived", "Humerus"]
P = data.loc[data["Status"] == "Perished", "Humerus"]

objective

After a severe winter storm, some of the collected sparrows survived and some died. The question is whether the ones that died had a shorter humerus.

data description

The data covers 59 adult male sparrows, of which 24 died and 35 survived. Humerus length is recorded in thousandths of an inch.

point estimates

I estimate the mean with the sample mean

The variance with the sample variance

And the standard deviation with the sample standard deviation

S_n = S.shape[0]
S_mean = np.mean(S)
S_var = np.var(S, ddof=1)
S_std = np.sqrt(S_var)
S_med = np.median(S)

P_n = P.shape[0]
P_mean = np.mean(P)
P_var = np.var(P, ddof=1)
P_std = np.sqrt(P_var)
P_med = np.median(P)


pd.DataFrame({
    "Frequency": [S_n, P_n],
    "Sample mean": [S_mean, P_mean],
    "Sample variance": [S_var, P_var],
    "Sample standard deviation": [S_std, P_std],
    "Median": [S_med, P_med],
}, index=["Survived", "Perished"]).round(4).T

distribution function estimate

For each group, I estimate the density with a histogram and the distribution function with the empirical distribution function.

histogram

I estimate the density with a histogram using \(k\) bins over the range of the data.

empirical distribution function

The empirical distribution function gives the estimate for the distribution function.

def process_bins(vals, bins, normalize=True):
    density_est = []
    my_bins = []
    
    for i in range(len(vals)):
        if i == len(vals) - 1:
            bin = f"[{bins[i]:.0f}, {bins[i + 1]:.0f}]"
        else:
            bin = f"[{bins[i]:.0f}, {bins[i + 1]:.0f})"
        my_bins.append(bin)

    for val in vals:
        if normalize:
            val = val / np.sum(vals)
        density_est.append(val)

    return density_est, my_bins

def histogram(data, bin_count):
    fig, axes = plt.subplots(1, 2, figsize=(10, 5), layout="constrained")
    
    counts, bins = np.histogram(data, bins=bin_count)
    bins = bins.round()
    
    ax = axes[0]
    sns.kdeplot(data, ax=ax, linewidth=1, color="black", alpha=0.3)
    sns.histplot(data, bins=bins, stat="density", ax=ax)
    sns.rugplot(data, ax=ax, color="black")
    ax.set_title("Histogram")
    
    ax = axes[1]
    sns.ecdfplot(data, stat="proportion", ax=ax, linewidth=1.5)
    sns.kdeplot(data, cumulative=True, ax=ax, linewidth=1, color="black", alpha=0.3)
    sns.rugplot(data, ax=ax, color="black", height=0.015)
    ax.set_title("Empirical distribution function")

    plt.show()

    density_est, my_bins =  process_bins(counts, bins)
    density_table = pd.DataFrame([my_bins, np.round(density_est, 4)], index=["Bin", "Estimation \(F(x)\)"], columns=[""]*len(my_bins))

    ecdf = st.ecdf(data)
    distrib_est, my_bins =  process_bins(ecdf.cdf.probabilities[:-1], ecdf.cdf.quantiles, normalize=False)
    distrib_table = pd.DataFrame([my_bins, np.round(distrib_est, 4)], index=["Bin", "Estimation \(F(x)\)"], columns=[""]*len(my_bins))
    
    return density_table, distrib_table

survived

density_table, distrib_table = histogram(S, 9)

density estimation

From the histogram, the density estimate is

density_table

estimation of the distribution function

From the empirical distribution function, the distribution function estimate is

display(distrib_table.iloc[:, :9])
display(distrib_table.iloc[:, 9:18])
display(distrib_table.iloc[:, 18:])

perished

density_table, distrib_table = histogram(P, 6)

density estimation

From the histogram, the density estimate is

density_table

estimation of the distribution function

From the empirical distribution function, the distribution function estimate is

display(distrib_table.iloc[:, :7])
display(distrib_table.iloc[:, 7:14])
display(distrib_table.iloc[:, 14:21])

statistical model

For each group, I fit normal, exponential, and uniform distributions by estimating their parameters. Then I overlay the fitted densities on the histograms and see which one looks closest.

parameter estimation

I estimate the exponential and uniform parameters using the method of moments.

The moments are estimated by sample moments.

S_m1 = 1 / S_n * np.sum(S)
S_m2 = 1 / S_n * np.dot(S, S)

P_m1 = 1 / P_n * np.sum(P)
P_m2 = 1 / P_n * np.dot(P, P)

pd.DataFrame({
    "Survived": [S_m1, S_m2],
    "Perished": [P_m1, P_m2],
}, index=["Sample moment \(m_1\)", "Sample moment \(m_2\)"]).round(4)

Now I calculate the normal distribution parameters from the earlier estimates.

Normal distribution \(\mathcal{N}(\mu, \sigma^2)\)

I express the parameters of the remaining two models as functions of moments.

Exponential distribution \(\text{Exp}(\lambda, L)\)

Uniform distribution \(\mathcal{U}(a, b)\)

def add_est(data, m1, m2, bin_count, mean, std, title, ax):
    X = np.linspace(data.min()-10, data.max()+10, 1000)
    counts, bins = np.histogram(data, bins=bin_count)
    bins = bins.round()
    
    lamb_hat = 1 / np.sqrt(m2 - m1*m1)
    l_hat = m1 - np.sqrt(m2 - m1*m1)
    
    a_hat = m1 - np.sqrt(3 * (m2 - m1*m1))
    b_hat = m1 + np.sqrt(3 * (m2 - m1*m1))

    Y_norm = st.norm.pdf(X, loc=mean, scale=std)
    
    X_exp = np.linspace(l_hat, data.max()+10, 1000)
    Y_exp = st.expon.pdf(X_exp, loc=l_hat, scale=1/lamb_hat)
    
    sns.histplot(data, bins=bins, stat="density", ax=ax, color="grey", alpha=0.3)
    
    sns.lineplot(x=X, y=Y_norm, ax=ax, label=f"N ( {mean:.0f}, {std*std:.0f} )", linewidth=2.5, color="royalblue")
    
    sns.lineplot(x=X_exp, y=Y_exp, ax=ax, label=f"Exp ( {lamb_hat:.4f}, {l_hat:.0f} )", linewidth=2.5, color="darkorange")
    sns.lineplot(x=[X.min(), l_hat], y=0.0003, ax=ax, linewidth=2.5, color="darkorange")
    
    sns.lineplot(x=[a_hat, b_hat], y=1 / (b_hat - a_hat), ax=ax, label=f"U ( {a_hat:.0f}, {b_hat:.0f} )", linewidth=2.5, color="green")
    sns.lineplot(x=[X.min(), a_hat], y=0.0003, ax=ax, linewidth=2.5, color="green")
    sns.lineplot(x=[b_hat, X.max()], y=0.0003, ax=ax, linewidth=2.5, color="green")
    
    ax.legend(loc="upper right")
    ax.set_title(title)

fig, axes = plt.subplots(1, 2, figsize=(10, 5), layout="constrained", sharey=True)
add_est(S, S_m1, S_m2, 9, S_mean, S_std, "Survived", axes[0])
add_est(P, P_m1, P_m2, 6, P_mean, P_std, "Perished", axes[1])

Conclusion

Both groups look approximately normal, and the normal distribution fits best of the three.

simulation

For each group, I simulate 100 values from the fitted normal distribution and compare the histogram to the original data.

def plot_norm(data, mean, std, bin_count):
    fig, axes = plt.subplots(1, 2, figsize=(10, 5), layout="constrained", sharey=True)

    counts, bins = np.histogram(data, bins=bin_count)
    bins = bins.round()
    w = bins[1] - bins[0]
    
    ax = axes[0]
    sns.histplot(data, bins=bins, stat="density", ax=ax)
    ax.set_title("Original data")

    X = np.linspace(data.min()-10, data.max()+10, 1000)
    while X.max() > bins.max():
        bins = np.append(bins, [bins.max() + w])
    while X.min() < bins.min():
        bins = np.insert(bins, 0, bins.min() - w)
    
    ax = axes[1]
    new_data = st.norm.rvs(mean, std, 100, random_state=42)
    sns.histplot(new_data, bins=bins, stat="density", ax=ax, color="tomato")
    ax.set_title("Simulated data")
    ax.set_xlabel("Humerus")

survived

plot_norm(S, S_mean, S_std, 9)

perished

plot_norm(P, P_mean, P_std, 6)

Conclusion

The histograms of the simulated samples look similar to the observed ones.

confidence interval

Now I calculate a two-sided \(95\%\) confidence interval for the mean in each group.

The two-sided \(100 \cdot (1 - \alpha)\%\) confidence interval for \(\mu\) with unknown variance is calculated as

def interval(data, n, mean, std, alpha):
    t_a_half = st.t.isf(alpha / 2, n - 1)
    delta = t_a_half * std / np.sqrt(n)
    return np.array([mean - delta, mean + delta])

S_ci = interval(S, S_n, S_mean, S_std, 0.05).round(4)
P_ci = interval(P, P_n, P_mean, P_std, 0.05).round(4)

pd.DataFrame({
    "Survived": [f"( {S_ci[0]} , {S_ci[1]} )"],
    "Perished": [f"( {P_ci[0]} , {P_ci[1]} )"],
}, index=["Two-sided 95% confidence interval for \(\mu\)"]).T

mean value test

For each group, I test at the \(5\%\) level whether the mean equals \(\mu_0 = K\), against the two-sided alternative. \(K\) is a parameter of the course assignment and is far below the observed lengths (around 740), so the outcome of this test is clear in advance.

pd.DataFrame([K], index=["K"], columns=[""])

survived

Test

At significance level \(\alpha = 5\%\), the test is

The corresponding confidence interval is

The tested value \(\mu_0 = 16\) is not inside the interval.

Test conclusion

At the \(5\%\) level, we reject the null hypothesis. The true mean humerus length of surviving sparrows is not 16.

perished

Test

At significance level \(\alpha = 5\%\), the test is

The corresponding confidence interval is

The tested value \(\mu_0 = 16\) is not inside the interval.

Test conclusion

At the \(5\%\) level, we reject the null hypothesis. The true mean humerus length of sparrows that died is not 16.

test for equality of means

At the \(5\%\) level, I test whether the two groups have the same mean.

two-sample t-test

pd.DataFrame([np.round(S_std / P_std, 4)], index=["\(s_{\text{survived}} / s_{\text{perished}}\)"], columns=[""])

The variances look similar enough to assume equal variances.

I use a two-sample t-test for normal distributions.

The test statistic \(T\) is compared to the critical value from Student’s \(T\)-distribution.

Where

SP_std = np.sqrt(( (S_n - 1) * S_var + (P_n - 1) * P_var )/( S_n + P_n - 2 ))
T_stat = (S_mean - P_mean) / SP_std * np.sqrt((S_n * P_n)/(S_n + P_n))
T_crit_05_57 = st.t.isf(0.05, 57)
T_crit_025_57 = st.t.isf(0.025, 57)

pd.DataFrame(np.array([T_stat, T_crit_05_57]).round(4),
             index=["\(T\)", "\(t_{0.05, 57}\)"], columns=["Value"])

Test

At significance level \(\alpha = 5\%\), the test is

The test statistic is \(T = 1.777\).

Since

We reject the null hypothesis at the \(5\%\) level.

conclusion

Surviving sparrows had a longer mean humerus than those that died. The one-sided test rejects equal means at 5% in favor of the survivors.

alternative conclusion

pd.DataFrame(np.array([T_stat, T_crit_025_57]).round(4),
             index=["\(T\)", "\(t_{0.025, 57}\)"], columns=["Value"])

With a two-sided alternative instead:

The two-sided \(p\)-value is above \(5\%\), so the two-sided test does not reject \(H_0\). The result therefore depends on the one-sided alternative, which is only justified if the direction was hypothesized before looking at the data.

p_greater = st.mstats.ttest_ind(S, P, equal_var=True, alternative="greater")
p_two_sided = st.mstats.ttest_ind(S, P, equal_var=True, alternative="two-sided")
pd.DataFrame(np.array([p_greater, p_two_sided]).round(4),
             index=["One-sided test", "Two-sided test"], columns=["T statistic", "p-value"])