10: Two-Sample Tests

Gather Data

The data are from the US Center for Disease Control and Prevention (CDC) and is stored at OpenIntro. The CDC annually surveys the health status and habits of 350,000 people in the US, and the dataset is a random sample of 20,000 people from the survey conducted in 2000.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

url = "https://raw.githubusercontent.com/DeAnzaDataScience/CIS11/refs/heads/main/datasets_notes/cdc.csv"
health = pd.read_csv(url)
display(health.head())
print("Number of rows:", health.shape[0])

Output?

Number of rows: 20000

The rows from left to right are: general health, any exercise in the past month, health plan coverage, smoke at least 100 cigarettes in respondent's life, height, weight, desired weight, age, and gender.

We will look at two features: exercise habit or exerany (0 for no exercise in past month, 1 for exercise in past month), and weight (a quantitative variable). We want to know whether people who exercise have a lower body weight distribution than people who don't exercise.

Exploratory Data Analysis (EDA)

It might seem reasonable that people who exercise have lower body weight than those who don't exercise. However, a person's weight is influenced by their height. Taller people naturally tend to weigh more than shorter people, no matter how much they exercise.

# find the BMI from the weight and height
BMI = health.weight / (health.height ** 2) * 703

# add to the health DataFrame
health['bmi'] = BMI
display(health.head(10))

Output?

Prepare Data

Next we concentrate on the two features, exercise habit and BMI, by creating a new DataFrame exercise that only has the exerany and bmi columns.

exercise_bmi = health[['exerany', 'bmi']]
display(exercise_bmi.head())
print("Number of rows:", exercise_bmi.shape[0])

Output?

We find the count and proportion of exercise and non-exercise groups.

# Find the exercise and no exercise groups
exercise = exercise_bmi[exercise_bmi["exerany"] == 1]
no_exercise = exercise_bmi[exercise_bmi["exerany"] == 0]

print("No exercise (0): count =", len(exercise), "  proportion =", round(len(exercise)/len(exercise_bmi),2))
print("Exercise (1): count =", len(no_exercise), "  proportion = ", round(len(no_exercise)/len(exercise_bmi),2))

Output?

No exercise (0): count = 14914   proportion = 0.75
Exercise (1): count = 5086   proportion =  0.25

Then we plot the distribution of the BMI of the two groups.

plt.figure(figsize=(6, 4))
plt.hist(no_exercise["bmi"], bins=30, density=True, alpha=0.5, label="No exercise")
plt.hist(exercise["bmi"], bins=30, density=True, alpha=0.5, label="Exercise")
plt.legend()
plt.title("BMI Distribution by Exercise Habit")
plt.xlabel("BMI")
plt.ylabel("Density")
plt.grid()
plt.show()

Output?

It looks like the BMI distributions overlap quite a bit for the exercise and non-exercise groups. We also see that there seems to be a right tail in the distribution, which means that there are some people with BMI values that are much higher than those of most people.

The Test Statistic

To compare the exercise and no exercise group's BMI, we find the difference in the mean BMI of the two groups.

# Calculate the mean BMI for each group
mean_bmi_exercise = exercise["bmi"].mean()
mean_bmi_no_exercise = no_exercise["bmi"].mean()
observed_bmi_difference = mean_bmi_exercise - mean_bmi_no_exercise

print("Mean BMI for exercise group:", round(mean_bmi_exercise, 2))
print("Mean BMI for no exercise group:", round(mean_bmi_no_exercise, 2))
print("Difference (exercise - no exercise):", round(observed_bmi_difference, 2))

Output?

Mean BMI for exercise group: 25.98
Mean BMI for no exercise group: 27.27
Difference (exercise - no exercise): -1.29

There is a small difference: the mean BMI of the exercise group is lower than the mean BMI of the no exercise group by 1.29.

Is the 1.29 difference due to random chance or is there evidence of a difference between the two groups?
We will use hypothesis testing to find the answer.

Hypothesis Testing

We'll use the same steps as in Chapter 9, but this time we're working with two data samples so our model for the null hypothesis will be built differently.

example = pd.DataFrame(exercise_bmi.head(7))
display(example)

Output?

shuffled_exerany = example["exerany"].sample(frac=1, replace=False)
display(shuffled_exerany)

Output?

5    1
2    1
4    0
0    0
6    1
3    1
1    0
Name: exerany, dtype: int64

# add the shuffled column to example,
# and reset the index so it matches the index of the DataFrame
example["shuffled_exerany"] = shuffled_exerany.reset_index(drop=True)
display(example)

Output?

def one_permutation(dataframe) :
  # create a copy of the original dataframe to add a new shuffled_treatment column
  dataframe_permutation = dataframe.copy()

  # create a shuffled_exerany column
  shuffled_exerany = dataframe_permutation["exerany"].sample(frac=1, replace=False)
  dataframe_permutation["shuffled_exerany"] = shuffled_exerany.reset_index(drop=True)

  # separate the data into two shuffled groups
  exercise = dataframe_permutation[dataframe_permutation["shuffled_exerany"] == 1]
  no_exercise = dataframe_permutation[dataframe_permutation["shuffled_exerany"] == 0]

  # calculate the mean BMI of each group
  exercise_mean_bmi = exercise["bmi"].mean()
  no_exercise_mean_bmi = no_exercise["bmi"].mean()

  # return the difference between the mean BMI values
  return exercise_mean_bmi - no_exercise_mean_bmi

Output?

# create a list to store the 1000 differences
differences = []

for i in range(1000):
    one_difference = one_permutation(exercise_bmi)
    differences.append(one_difference)
mean_bmi_difference_array = np.array(differences)

Output?

plt.figure(figsize=(4, 4))
plt.hist(mean_bmi_difference_array, density=True, edgecolor='black')
plt.title("Distribution of Permutation Differences")
plt.xlabel("Difference in Mean BMI")
plt.ylabel("Density")
plt.grid()
plt.show()

Output?

observed_distance = abs(observed_bmi_difference)
simulated_distance = np.abs(mean_bmi_difference_array)
p_value = np.mean(simulated_distance >= observed_distance)
print(p_value)

Output?

0.0

Observational Data

It is important to note that our hypothesis testing above provides evidence that there is an association between BMI value and exercise habit. There seems to be a difference in BMI values among people who exercise and those who don't, but we cannot say that exercise causes the difference in BMI. This is because people were not randomly assigned to the exercise or no exercise group.

If people are randomly assigned to a group, the people who have healthy diets can be assigned to either group, and older and younger people can be assigned to either group. Random assignment tends to balance these other characteristics between the groups, making exercise habit one of the primary difference between the groups. This makes it possible to infer that exercise may cause a difference in mean BMI.

On the other hand, when people are surveyed for their health condition, as in our dataset, then the two groups can be different in many ways. For example, people who exercise may also have healthy diets or are younger. These factors, the confounding variables, may also influence BMI so it's not possible to say that exercise itself causes the difference in BMI.

The dataset for the exercise habit and mean BMI is observational data. They are data gathered by measuring or recording behaviors, events, or phenomena as they naturally occur, without manipulating any chracteristics of the data. Observational data can only show evidence of association between two factors, they cannot be used to determine cause and effect.

Causality and A/B Testing

We now look at an example of experimental data, which is the opposite of observational data. In experimental data people are randomly assigned to two groups, and this makes it possible to infer a cause and effect relationship as discussed above.

The Criteo Uplift Prediction dataset contains data from online advertising experiments conducted by Criteo, a digital advertising company. Users were randomly assigned to either an advertising treatment group or a control group. Criteo could show advertisements to users in the treatment group but withheld its advertising from users in the control group. Then the company recorded whether the user visited the advertiser's website.

This type of randomized experiment is commonly called an A/B test. A/B testing is a type of randomized experiment used to compare two alternatives. Participants are randomly assigned to one of two groups:

• Group A (control): receives the existing or control condition.
• Group B (treatment): receives the new condition being tested.

We then compare an outcome from the two groups to determine whether the treatment made a difference. Because participants are randomly assigned, A/B testing can provide evidence that the treatment caused a difference in the outcome.

A/B testing is widely used in business, technology, and marketing to determine whether a change causes an improvement in a desired outcome.

url = "https://raw.githubusercontent.com/DeAnzaDataScience/CIS11/refs/heads/main/datasets_notes/ad.csv"
ad = pd.read_csv(url)
print("Last 10 rows")
display(ad.tail(10))
print("Number of rows:", ad.shape[0])

Output?

Last 10 rows

Number of rows: 20000

ad.value_counts()

Output?

treatment  visit
1          0        16136
0          0         2916
1          1          828
0          1          120
Name: count, dtype: int64

# select a group's visit column and find the mean
control_visit_proportion = round(ad[ad["treatment"] == 0]["visit"].mean(), 4)
treatment_visit_proportion = round(ad[ad["treatment"] == 1]["visit"].mean(), 4)
observed_visit_difference = round(treatment_visit_proportion - control_visit_proportion, 4)
print("Proportion of control group that visited the website", control_visit_proportion)
print("Proportion of treatment group that visited the website", treatment_visit_proportion)
print("Observed difference:", observed_visit_difference)

Output?

Proportion of control group that visited the website 0.0395
Proportion of treatment group that visited the website 0.0488
Observed difference: 0.0093

def one_permutation(dataframe) :
  # create a copy of the original dataframe to add a new shuffled_treatment column
  dataframe_permutation = dataframe.copy()

  # create a shuffled_treatment column
  shuffled_treatment = dataframe_permutation["treatment"].sample(frac=1, replace=False)
  dataframe_permutation["shuffled_treatment"] = shuffled_treatment.reset_index(drop=True)

  # separate the data into two shuffled groups
  treatment = dataframe_permutation[dataframe_permutation["shuffled_treatment"] == 1]
  control = dataframe_permutation[dataframe_permutation["shuffled_treatment"] == 0]

  # calculate the mean BMI of each group
  treatment_proportion = treatment["visit"].mean()
  control_proportion = control["visit"].mean()

  # return the difference between the mean BMI values
  return treatment_proportion - control_proportion

Output?

# create a list to store the 1000 differences
differences = []

for i in range(1000):
    one_difference = one_permutation(ad)
    differences.append(one_difference)
visit_difference_array = np.array(differences)

Output?

plt.figure(figsize=(6, 4))
plt.hist(visit_difference_array, density=True, edgecolor='black')
plt.title("Distribution of Visit Proportion Differences")
plt.xlabel("Difference in Visit Proportion")
plt.ylabel("Density")
plt.grid()
plt.show()

Output?

observed_distance = abs(observed_visit_difference)
simulated_distance = np.abs(visit_difference_array)
p_value = np.mean(simulated_distance >= observed_distance)
print(p_value)

Output?

0.022

Total Variation Distance

In the previous example we worked with a dataset with two categories of results: visit or no visit. With two categories, the statistic that we want to measure is simple: it's the difference in visit proportion between the two groups. We don't have to separately consider the no-visit proportion because it is determined by the visit proportion.

But suppose there are three categories of results: website visit, no visit, and purchase. Now each group has three proportions, and finding the difference in the proportions of only one category is not sufficient. This is where the Total Variation Distance or TVD comes in. TVD is a way to measure the difference between two distributions while taking into account the proportion of all categories.

To calculate the TVD:

• for each category, find the absolute value of the difference in proportion
• add the absolute values together
• divide the sum by 2

To see how TVD is calculated, we create an example of the three categories of result: visit, no visit, and purchase. We also have their proportions for both treatment and control groups, and the difference between the proportions:

example = pd.DataFrame(
    {
        "control": [0.8, 0.14, 0.06],
        "treatment": [0.7, 0.2, 0.1],
        "difference":[-0.10, 0.06, 0.04]
    },
    index=["no visit", "visit", "purchase"],
)
display(example)

Output?

We note that for the control and treatment groups, each proportion adds up to 1.0.
We also note that the decrease in "no visit" (-0.10) is redistributed to "visit" and "purchase" which together increase by the same amount (0.06 and 0.04). Therefore 0.10 or 10% of the proportion has been redistributed among the categories. We could think of this as 10% of the people changing categories to turn the control distribution into the treatment distribution.

If we simply add the difference of the two groups, we would get: −0.10 + 0.06 + 0.04 = 0.0
The differences cancel out because the proportions in both groups have to add up to 1.0.

Instead we add the absolute value of the difference so the decreases and increases do not cancel each other, but this results in a value that's doubled what it should be: 0.10 + 0.06 + 0.04 = 0.20
This counts the redistribution twice: once when the 0.10 leaves no visit, and again when the same 0.10 appears as 0.06 in visit and 0.04 in purchase.

Therefore we divide by 2 to get the actual difference between the two groups.

# calculate the TVD

absolute_difference = np.abs(example["treatment"] - example["control"])
total = np.sum(absolute_difference)
tvd = total / 2
print("TVD:", round(tvd, 2))

Output?

TVD: 0.1

The TVD gives us one number that summarizes the differences across all the categories.

A TVD of 0 means the distributions are identical. The farther the TVD is from 0, the more different the two distributions are.

Summary

In this chapter we apply our sampling and simulation skills to two-sample hypothesis testing, which is used to compare two groups of data. We also see the importance of random assignment to groups, which allows us to determine whether a treatment may have caused a difference in the outcome.