1. Introduction & Dataset

In David Hamilton’s classic 1987 dataset, the dependent variable \(Y\) depends on the predictor variables collectively as a group, but exhibits almost zero linear correlation when evaluated individually. This statistical anomaly is known as Complete Suppression.

# 1. Input the original Hamilton (1987) dataset
hamilton_data <- data.frame(
  Y  = c(12.37, 12.66, 12.00, 11.93, 11.06, 13.03, 13.13, 11.44, 12.86, 10.84, 11.20, 11.56, 10.83, 12.63, 12.46),
  X1 = c(2.23,  2.57,  3.87,  3.10,  3.39,  2.83,  3.02,  2.14,  3.04,  3.26,  3.39,  2.35,  2.76,  3.90,  3.16),
  X2 = c(9.66,  8.94,  4.40,  6.64,  4.91,  8.52,  8.04,  9.05,  7.71,  5.11,  5.05,  8.51,  6.59,  4.90,  6.96)
)

# Display the 15 observations in a clean table
knitr::kable(hamilton_data, caption = "Hamilton (1987) Original Data Table")
Hamilton (1987) Original Data Table
Y X1 X2
12.37 2.23 9.66
12.66 2.57 8.94
12.00 3.87 4.40
11.93 3.10 6.64
11.06 3.39 4.91
13.03 2.83 8.52
13.13 3.02 8.04
11.44 2.14 9.05
12.86 3.04 7.71
10.84 3.26 5.11
11.20 3.39 5.05
11.56 2.35 8.51
10.83 2.76 6.59
12.63 3.90 4.90
12.46 3.16 6.96

2. Bivariate Correlation Matrix

If we look at the pairwise correlations, \(X_1\) shows virtually zero linear correlation with \(Y\) (\(r = 0.0025\)). Meanwhile, the two predictors \(X_1\) and \(X_2\) are strongly negatively correlated with each other (\(r = -0.8998\)).

print("Bivariate Correlation Matrix:")
## [1] "Bivariate Correlation Matrix:"
print(cor(hamilton_data))
##              Y           X1         X2
## Y  1.000000000  0.002497966  0.4340688
## X1 0.002497966  1.000000000 -0.8997765
## X2 0.434068758 -0.899776481  1.0000000

3. Individual Simple Linear Regressions

3.1 Individual Model for X1 (\(Y \sim X_1\))

Notice that the Adjusted R-squared is negative (-0.0769). Because \(X_1\) fails to explain any meaningful variance in \(Y\), the sample size and variable count penalty causes the calculation to fall below zero.

model1 <- lm(Y ~ X1, data = hamilton_data)
summary(model1)
## 
## Call:
## lm(formula = Y ~ X1, data = hamilton_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.16910 -0.67912 -0.00326  0.64412  1.12993 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.988755   1.266891   9.463  3.4e-07 ***
## X1           0.003747   0.416083   0.009    0.993    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8324 on 13 degrees of freedom
## Multiple R-squared:  6.24e-06,   Adjusted R-squared:  -0.07692 
## F-statistic: 8.112e-05 on 1 and 13 DF,  p-value: 0.993

ANOVA Table for \(Y \sim X_1\)

anova(model1)
## Analysis of Variance Table
## 
## Response: Y
##           Df Sum Sq Mean Sq F value Pr(>F)
## X1         1 0.0001 0.00006   1e-04  0.993
## Residuals 13 9.0085 0.69296

3.2 Individual Model for X2 (\(Y \sim X_2\))

\(X_2\) on its own explains a modest portion of the variance in \(Y\) (\(R^2 = 0.1884\)).

model2 <- lm(Y ~ X2, data = hamilton_data)
summary(model2)
## 
## Call:
## lm(formula = Y ~ X2, data = hamilton_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08999 -0.63345  0.00023  0.61458  1.04033 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.6319     0.8109  13.111 7.18e-09 ***
## X2            0.1955     0.1125   1.737    0.106    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7499 on 13 degrees of freedom
## Multiple R-squared:  0.1884, Adjusted R-squared:  0.126 
## F-statistic: 3.018 on 1 and 13 DF,  p-value: 0.106

4. The Joint Multiple Linear Regression (\(Y \sim X_1 + X_2\))

When we combine \(X_1\) and \(X_2\) in a single model, the statistical magic happens. Both Multiple and Adjusted R-squared jump to a perfect 1.000 (100%). \(X_1\) acts as a Suppressor Variable—it purges the irrelevant variance out of \(X_2\), allowing the true collective predictive power to emerge cleanly.

model <- lm(Y ~ X1 + X2, data = hamilton_data)
print("Multiple Regression Summary:")
## [1] "Multiple Regression Summary:"
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2, data = hamilton_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.013632 -0.009451 -0.002279  0.008630  0.016325 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.515414   0.061142  -73.85   <2e-16 ***
## X1           3.097008   0.012274  252.31   <2e-16 ***
## X2           1.031859   0.003684  280.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01072 on 12 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 3.922e+04 on 2 and 12 DF,  p-value: < 2.2e-16

5. Interactive 3D Scatter Plot

Instructors/Students: Click and drag your mouse over the plot to rotate it in 3D space.

Notice how when you rotate the view to a specific angle, all 15 data points align perfectly into a completely flat plane. This geometric perfection is why our combined model achieves an \(R^2 = 1.0\).