Exercises for Section 3.5 Two-sample tests#
This notebook contains the solutions to the exercises from Section 3.5 Two-sample tests in the No Bullshit Guide to Statistics.
Notebook setup#
# Ensure required Python modules are installed
%pip install --quiet numpy scipy seaborn pandas ministats
Note: you may need to restart the kernel to use updated packages.
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Figures setup
plt.clf() # needed otherwise `sns.set_theme` doesn't work
sns.set_theme(
context="paper",
style="whitegrid",
palette="colorblind",
rc={"figure.figsize": (6, 1.6)},
)
%config InlineBackend.figure_format = "retina"
<Figure size 640x480 with 0 Axes>
# Simplified int and float __repr__
np.set_printoptions(legacy='1.25')
# Download datasets/ directory if necessary
from ministats import ensure_datasets
ensure_datasets()
datasets/ directory already exists.
\(\def\stderr#1{\mathbf{se}_{#1}}\) \(\def\stderrhat#1{\hat{\mathbf{se}}_{#1}}\) \(\newcommand{\Mean}{\textbf{Mean}}\) \(\newcommand{\Var}{\textbf{Var}}\) \(\newcommand{\Std}{\textbf{Std}}\) \(\newcommand{\Freq}{\textbf{Freq}}\) \(\newcommand{\RelFreq}{\textbf{RelFreq}}\) \(\newcommand{\DMeans}{\textbf{DMeans}}\) \(\newcommand{\Prop}{\textbf{Prop}}\) \(\newcommand{\DProps}{\textbf{DProps}}\) $\( \newcommand{\CI}[1]{\textbf{CI}_{#1}} \newcommand{\CIL}[1]{\textbf{L}_{#1}} \newcommand{\CIU}[1]{\textbf{U}_{#1}} \newcommand{\ci}[1]{\textbf{ci}_{#1}} \newcommand{\cil}[1]{\textbf{l}_{#1}} \newcommand{\ciu}[1]{\textbf{u}_{#1}} \)\( (this cell contains the macro definitions like \)\stderr{\overline{\mathbf{x}}}\(, \)\stderrhat{}\(, \)\Mean$, …)
Exercises#
E3.12#
In an earlier exercise, we calculated the difference between mean sleep scores for doctors working in rural and urban locations: \(\hat{d} = 7\). In this exercise, we want to describe the uncertainty in this estimate.
a) Obtain an approximation to the sampling distribution of the difference between means in terms of Student’s \(t\)-distribution centred at \(\hat{d}\), with scale \(\stderrhat{\hat{d}}\), and \(\nu_d\) degrees of freedom.
b) Approximate the sampling distribution using the bootstrap.
c) Compare your answers from part a) and b) graphically.
# Load the doctors dataset and get the "rur" and "urb" scores
doctors = pd.read_csv("datasets/doctors.csv")
scoresR = doctors[doctors["loc"]=="rur"]["score"]
scoresU = doctors[doctors["loc"]=="urb"]["score"]
# observed value
dscores = scoresR.mean() - scoresU.mean()
dscores
6.992885375494076
# a) Analytical approximation
# Obtain the sample sizes and stds for the two groups
nR = ...
stdR = ...
nU = ...
stdU = ...
# Compute the standard error of the difference between means
seDhatscores = ...
# Calculate the degrees of freedom parameter
from ministats import calcdf
...
# Build a probability model based on Student's t-distribution
from scipy.stats import t as tdist
...
# Plot the PDF of the sampling distribution `rvDscores`
from ministats import plot_pdf
...
Ellipsis
# b) Bootstrap estimate
# Set the random seed to a particular value to get a reproducible result
np.random.seed(51)
# Compute bootstrap estimates for mean in each group
from ministats import gen_boot_dist
meanR_boot = ...
meanU_boot = ...
# Compute the difference between means of the bootstrap samples
dmeans_boot = ...
# Plot a histogram of the bootstrap distribution
...
Ellipsis
# c) Plot of the PDF of the analytical model
# and the bootstrap distribution on the same axis
...
Ellipsis
E3.19#
Construct a 95% confidence interval for the difference between mean scores
of doctors living in rural ("rur") and urban ("urb") locations
in the doctors dataset
datasets/doctors.csv.
a) Use an analytical approximation based on Student’s \(t\)-distribution.
b) Use the percentile bootstrap method.
# Define the `dmeans` estimator function
def dmeans(xsample, ysample):
dhat = np.mean(xsample) - np.mean(ysample)
return dhat
# Load the doctors dataset and get the "rur" and "urb" scores
doctors = pd.read_csv("datasets/doctors.csv")
scoresR = doctors[doctors["loc"]=="rur"]["score"]
scoresU = doctors[doctors["loc"]=="urb"]["score"]
# Calculate the difference between mean scores
dscores = dmeans(scoresR, scoresU)
dscores
6.992885375494076
a) Analytical approximation#
# Obtain the sample sizes and stds of the two groups
nR = ...
stdR = ...
nU = ...
stdU = ...
# Calculate the standard error of the difference between means
...
# Calculate the degrees of freedom
from ministats import calcdf
dfD = ...
# Find the 0.025 and 0.975 quantiles of the standard t-dist.
from scipy.stats import t as tdist
...
# Construct a 95% CI for the difference between means
...
Ellipsis
b) Bootstrap estimate#
# Set the random seed to a particular value to get a reproducible result
np.random.seed(43)
# Compute bootstrap estimates for mean in each group
...
# Compute the difference between bootstrap means
dmeans_boot = ...
# Construct a 95% CI for the difference between means
...
Ellipsis
E3.20#
Construct an 80% confidence interval for the difference between means
of the debate ("debate") and lecture ("lecture") curriculum variants in the students dataset.
a) Use an analytical approximation based on Student’s \(t\)-distribution.
b) Use the percentile bootstrap method.
# Load the students dataset and get the "debate" and "lecture" scores
students = pd.read_csv("datasets/students.csv")
scoresD = students[students["curriculum"]=="debate"]["score"]
scoresL = students[students["curriculum"]=="lecture"]["score"]
# Calculate the observed difference between mean scores
dhat = dmeans(scoresD, scoresL)
dhat
8.319642857142867
a) Analytical approximation#
# Obtain the sample sizes and stds of the two groups
nD = ...
stdD = ...
nL = ...
stdL = ...
# Calculate the standard error of the difference between means
seDscores = ...
# Calculate the degrees of freedom
from ministats import calcdf
...
# Find the 0.1 and 0.9 quantiles of the standard t-dist.
from scipy.stats import t as tdist
...
# Construct an 80% CI for the difference between means
...
Ellipsis
b) Bootstrap estimate#
# Set the random seed to a particular value to get a reproducible result
np.random.seed(42)
# Compute bootstrap estimates for mean in each group
...
# Compute the difference between bootstrap means
...
# Construct an 80% CI for the difference between means
...
Ellipsis
E3.38#
Use the permutation test to check for a difference between sleep scores
of doctors living in rural (rur) and urban (urb) locations
in the doctors dataset
datasets/doctors.csv.
a) Formulate two competing hypotheses \(H_0\) and \(H_A\).
b) Calculate the observed difference between means \(\hat{d}\).
c) Simulate \(10\,000\) observations
from the distribution of the difference between means \(\widehat{D}_0\)
by calling resample_under_H0.
d) Compute the \(p\)-value of \(\hat{d}\) (or a more extreme value) under \(\widehat{D}_0\).
e) Make a decision if location affects sleep scores using \(\alpha\!=\!0.05\).
# Load the doctors dataset and select the two groups
doctors = pd.read_csv("datasets/doctors.csv")
scoresU = doctors[doctors["loc"]=="urb"]["score"]
scoresR = doctors[doctors["loc"]=="rur"]["score"]
def dmeans(xsample, ysample):
dhat = np.mean(xsample) - np.mean(ysample)
return dhat
# b) Calculate the observed difference between means
dhat = ...
# Set the random seed to a particular value to get a reproducible result
np.random.seed(42)
# c) Simulate observations from the distribution of the difference between means under H0
from ministats.hypothesis_tests import resample_under_H0
P = 10000
pdhats0 = []
for j in range(0, P):
...
# d) Compute the $p$-value of `dhat` (or a more extreme value)
from ministats import tailvalues
tails = ...
pvalue = ...
E3.39#
Apply the permutation test to compare the treat and ctrl groups in the IQ scores dataset
datasets/iqs2.csv.
# Load the dataset and select the two groups
iqs2 = pd.read_csv("datasets/iqs2.csv")
treated = iqs2[iqs2["group"]=="treat"]["iq"]
controls = iqs2[iqs2["group"]=="ctrl"]["iq"]
# Calculate the observed difference between means
dhat_iqs = dmeans(treated, controls)
dhat_iqs
1.5577507598784166
# Set the random seed to a particular value to get a reproducible result
np.random.seed(43)
# Obtain the sampling distribution of the difference between means under H0
from ministats.hypothesis_tests import resample_under_H0
pdhats0 = ...
# Compute the p-value of `dhat` (or a more extreme value)
from ministats import tailvalues
tails = ...
pvalue = ...
E3.40#
Compare the sleep scores for rural and urban in the doctors dataset
(datasets/doctors.csv).
Use Welch’s two-sample \(t\)-test based on the two-sided alternative hypothesis.
a) Formulate two competing hypotheses \(H_0\) and \(H_A\).
b) Calculate the difference between means \(\hat{d} = \mathbf{dmeans}(\mathbf{x}_R,\mathbf{x}_U)\) and the sample standard deviations \(s_{\mathbf{x}_R}\) and \(s_{\mathbf{x}_U}\).
c) Compute the estimated standard error of the sampling distribution \(\stderr{\hat{d}}\) and the degrees of freedom parameter \(\nu_d\).
d) Calculate the \(t\)-statistic.
e) Compute the \(p\)-value of the \(t\)-statistic under \(T_0 \sim \mathcal{T}(\nu_d)\).
f) Check your answer using the ministats function ttest_dmeans.
g) Check your answer using the scipy.stats function ttest_ind.
# Load the doctors dataset and select the two groups
doctors = pd.read_csv("datasets/doctors.csv")
scoresU = doctors[doctors["loc"]=="urb"]["score"]
scoresR = doctors[doctors["loc"]=="rur"]["score"]
# b) Calculate the observed difference between means ...
dhat = ...
# ... the sample sizes ...
nR = ...
nU = ...
# ... and the standard deviations
stdR = ...
stdU = ...
# c) Compute the standard error of the estimator D
seD = ...
# ... and the degrees of freedom
from ministats import calcdf
...
Ellipsis
# d) Compute the value of the t-statistic
obst = ...
# e) Calculate the p-value
from scipy.stats import t as tdist
rvT0 = ...
from ministats import tailprobs
pvalue = ...
# f) Check your answer using the `ministats` function `ttest_dmeans`
from ministats import ttest_dmeans
...
Ellipsis
# g) Check your answer using the `scipy.stats` function `ttest_ind`
from scipy.stats import ttest_ind
...
Ellipsis
E3.41#
Apply Welch’s two-sample \(t\)-test
to compare the treat and ctrl groups in the IQ scores dataset
datasets/iqs2.csv.
# Load the dataset and select the two groups
iqs2 = pd.read_csv("datasets/iqs2.csv")
treated = iqs2[iqs2["group"]=="treat"]["iq"]
controls = iqs2[iqs2["group"]=="ctrl"]["iq"]
# Calculate the observed difference between means
...
Ellipsis
E3.42#
In Example 2T,
we applied Welch’s \(t\)-test to Charlotte’s students dataset
datasets/students.csv.
Repeat the analysis of this dataset using Student’s two-sample \(t\)-test with pooled variance.
Hint: Perform the calculations from scratch then check your answer using scipy.stats.ttest_ind used with the option equal_var=True.
# Load the students dataset and select the two groups
students = pd.read_csv("datasets/students.csv")
scoresD = students[students["curriculum"]=="debate"]["score"]
scoresL = students[students["curriculum"]=="lecture"]["score"]
E3.43#
Your friend Ben, who totally doesn’t have a drinking problem, decided to scientifically compare the price of a pint of beer in the West End and the East End of the city. He sent you the following summary statistics for the prices in the West: \(\overline{\mathbf{x}}_W=10\), \(s_{\mathbf{x}_W}=1\), and \(n_W=7\). The corresponding stats from the East are \(\overline{\mathbf{x}}_E=8\), \(s_{\mathbf{x}_E}=1.5\), and \(n_E=6\). Run Welch’s two-sample \(t\)-test on this data to detemine if there is a difference between the West and East prices.
# Sample statistics from the problem statement
meanW, stdW, nW = 10, 1, 7
meanE, stdE, nE = 8, 1.5, 6
# Compute the difference between means
dhat_beer = ...
# Compute the estimated standard error
seD = ...
# ... and the degrees of freedom
from ministats import calcdf
dfD = ...
# Compute the value of the t-statistic
obst = ...
# Calculate the p-value
from scipy.stats import t as tdist
rvT0 = ...
from ministats import tailprobs
pvalue = ...
from scipy.stats import ttest_ind_from_stats
...
Ellipsis