statistical analysis of data on natural selection

i’m looking at whether sparrow arm bone length has anything to do with surviving a winter storm. the dataset has humerus measurements for birds that made it and birds that didn’t. i’ll run through point estimates, distributions, fitting a model, simulation, confidence intervals, and a couple of 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 bad winter storm, some sparrows survived and some didn’t. i’m checking whether the ones that died had shorter arm bones.

data description

the data covers 59 adult male sparrows. 24 died and 35 survived.

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 roughly normal, and the normal fit seems to match best.

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 simulated histograms look pretty close to the real 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 the task parameter \(K\), against the two-sided alternative.

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 \(0.16 \text{in} \approx 4 \text{mm}\).

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 \(0.16 \text{in} \approx 4 \text{mm}\).

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"])

návrh

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

the data suggests surviving sparrows had a longer mean humerus than the ones that died. the one-sided test rejects equality in favor of survivors being larger.

alternative conclusion

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

if we look at it two-sided instead

the two-sided \(p\)-value is above \(5\%\), so \(H_A\) would be statistically insignificant.

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"])