SD01 — Statistical Thinking
Course: Statistics and Probability for Data Science (DS14270).
Source deck: Samer El Kababji, SD01 — Statistical Thinking.
This lesson uses a small fictional example to translate the ideas in the deck into R. It uses group averages rather than a regression model.
Suppose we want to know whether a new medication reduces systolic blood pressure after eight weeks.
The table below is constructed for teaching. Each row describes a group of patients rather than an individual patient.
blood_pressure <- data.frame(
medication = c("Standard", "Standard", "New", "New"),
age_group = c("Younger", "Older", "Younger", "Older"),
patients = c(20, 80, 90, 10),
mean_bp = c(130, 140, 125, 135)
)
blood_pressure medication age_group patients mean_bp
1 Standard Younger 20 130
2 Standard Older 80 140
3 New Younger 90 125
4 New Older 10 135Age group is a covariate because it is an additional variable that may help us understand differences in blood pressure.
The overall mean blood pressure for each medication combines its younger and older patients. Each age group contributes according to its patient count:
\[ \bar y=\frac{\sum_j n_j\bar y_j}{\sum_j n_j}. \]
Here \(j\) runs over the age groups within one medication, \(n_j\) is the patient count, and \(\bar y_j\) is that age group’s mean blood pressure. Apply the equation separately for New and Standard.
Worked example — comparing medications. For New, combine 90 younger patients with a mean of 125 mmHg and 10 older patients with a mean of 135 mmHg:
\[ \bar y_{\text{New}} =\frac{90\times125+10\times135}{90+10} =\frac{12600}{100} =126\text{ mmHg}. \]
For Standard, combine 20 younger patients with a mean of 130 mmHg and 80 older patients with a mean of 140 mmHg:
\[ \bar y_{\text{Standard}} =\frac{20\times130+80\times140}{20+80} =\frac{13800}{100} =138\text{ mmHg}. \]
These calculations produce the table below: weighted_bp is the numerator, patients is the denominator, and mean_bp is their quotient. Younger and older patients are combined within each medication, giving one row for New and one for Standard.
overall_means <- aggregate(
cbind(weighted_bp = blood_pressure$mean_bp * blood_pressure$patients,
patients = blood_pressure$patients),
by = list(medication = blood_pressure$medication),
FUN = sum
)
overall_means$mean_bp <- overall_means$weighted_bp / overall_means$patients
overall_means medication weighted_bp patients mean_bp
1 New 12600 100 126
2 Standard 13800 100 138overall_difference <- with(
overall_means,
mean_bp[medication == "New"] - mean_bp[medication == "Standard"]
)
overall_difference[1] -12The new-medication group has an average blood pressure 12 mmHg lower than the standard-medication group. This is an association in the observed data. It does not yet tell us how much of the difference was caused by the medication.
The treatment groups have very different age compositions.
age_counts <- xtabs(patients ~ medication + age_group, data = blood_pressure)
age_counts age_group
medication Older Younger
New 10 90
Standard 80 20prop.table(age_counts, margin = 1) age_group
medication Older Younger
New 0.1 0.9
Standard 0.8 0.2Most patients receiving the standard medication are older, while most patients receiving the new medication are younger. Older patients also have higher blood pressure in both medication groups. Age is therefore connected to both:
In this example, age is a confounder. It mixes the medication comparison with an age comparison.
We can compare medications separately within each age group:
\[ \text{difference within an age group} =\bar y_{\text{New}}-\bar y_{\text{Standard}}. \]
Worked example — compare within each age group. For younger patients, \(125-130=-5\) mmHg; for older patients, \(135-140=-5\) mmHg. A negative difference means the new-medication group has lower average blood pressure.
younger <- subset(blood_pressure, age_group == "Younger")
older <- subset(blood_pressure, age_group == "Older")
within_age <- data.frame(
age_group = c("Younger", "Older"),
difference = c(
younger$mean_bp[younger$medication == "New"] -
younger$mean_bp[younger$medication == "Standard"],
older$mean_bp[older$medication == "New"] -
older$mean_bp[older$medication == "Standard"]
)
)
within_age age_group difference
1 Younger -5
2 Older -5Within both age groups, the difference is only −5 mmHg. The larger overall difference of −12 mmHg partly reflected the unequal age composition.
bp_matrix <- xtabs(mean_bp ~ medication + age_group, data = blood_pressure)
barplot(bp_matrix, beside = TRUE, ylim = c(0, 150),
col = c("steelblue", "tomato"),
ylab = "Mean systolic blood pressure (mmHg)",
legend.text = rownames(bp_matrix),
args.legend = list(x = "topleft"))
This simple separation by age is useful for explaining confounding. In a real study, many covariates may matter at the same time, and some age groups may contain very few patients. More advanced statistical methods are then needed.
Two variables are associated when their values tend to differ together. An association can be measured from observational data, but it does not by itself show that changing one variable will change the other.
A causal effect asks a different question: what would happen to the same kind of patients under the new medication compared with the standard medication? We cannot give both medications to the same person under identical conditions, so study design matters.
In a randomized experiment, chance assigns the treatment. This helps make the treatment groups similar in age and other characteristics before treatment. An observational study can reveal useful associations, but differences between the groups may provide other explanations for the result.
The lesson from this example is:
Try it: change the number of older and younger patients so both medication groups have the same age composition. Recalculate the overall difference.