208251 — Regression Analysis and Non-Parametric Statistics: Lab Materials
Donlapark Ponnoprat
Department of Statistics
Faculty of Science, Chiang Mai University
Department of Statistics
Faculty of Science, Chiang Mai University
September 03, 2026
Every dataset used in these labs is already built into R (or generated by a short line of code).
How to use these lab materials
- Read the explanation text before running each code chunk.
- Pay attention to interpretation notes — these help you understand what the output means statistically.
- In each Practice section you will see worked examples; the Assignment section is for you to complete independently.
- Two datasets recur across several labs so you can build familiarity with them:
trees(Case Study A) andmtcars(Case Study B), both built into base R. A few others (attitude, and one simulated dataset) appear once, in the assignments.
Lab 1: Correlation Analysis
Objectives: Students are able to use R to:
- Perform descriptive statistics for two quantitative variables
- Construct and interpret a scatterplot
- Compute and interpret the Pearson correlation coefficient
- Perform a hypothesis test for the population correlation
- Find a confidence interval for the population correlation
Key concepts: The Pearson correlation coefficient \(r\) measures the strength and direction of the linear relationship between two quantitative variables \(X\) and \(Y\). It ranges from −1 to +1:
\[r = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum(x_i-\bar{x})^2 \; \sum(y_i-\bar{y})^2}}\]
| Value of \(|r|\) | Interpretation |
|---|---|
| 0.00 – 0.19 | Very weak |
| 0.20 – 0.39 | Weak |
| 0.40 – 0.59 | Moderate |
| 0.60 – 0.79 | Strong |
| 0.80 – 1.00 | Very strong |
The sign of \(r\) indicates direction: positive means both variables increase together; negative means one increases as the other decreases.
Important: Correlation measures linear association only. A correlation of 0 does not necessarily mean there is no relationship — only that there is no linear relationship.
Practice I: Cherry Tree Volume (Case Study A —
trees)
R comes with a built-in dataset called trees: 31 felled
black cherry trees, each with measured girth (trunk
diameter), height, and volume of
usable timber. No import is needed — trees is available the
moment R starts.
Variables:
Girth(trunk diameter, inches) andVolume(timber volume, cubic feet). Research question: Is there a significant linear correlation between trunk girth and timber volume?
Step 1 — Load and inspect the data
# trees is a built-in dataset - no download or setwd() required
data(trees)
head(trees) # first 6 rows## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...## Girth Height Volume
## Min. : 8.30 Min. :63 Min. :10.20
## 1st Qu.:11.05 1st Qu.:72 1st Qu.:19.40
## Median :12.90 Median :76 Median :24.20
## Mean :13.25 Mean :76 Mean :30.17
## 3rd Qu.:15.25 3rd Qu.:80 3rd Qu.:37.30
## Max. :20.60 Max. :87 Max. :77.00Note:
str()tells you how many observations (rows) and variables (columns) are in your data, and what type each variable is (numeric, character, factor, etc.).
Step 2 — Descriptive statistics
x <- trees$Girth # independent variable (predictor)
y <- trees$Volume # dependent variable (response)
# Summary statistics
summary(x)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.30 11.05 12.90 13.25 15.25 20.60## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.20 19.40 24.20 30.17 37.30 77.00## [1] 3.138139## [1] 16.43785Step 3 — Create a scatterplot
A scatterplot is always the first step — it reveals the form, direction, and strength of the relationship before any numerical summary.
plot(x, y,
xlab = "Girth (in.)",
ylab = "Volume (cu. ft.)",
main = "Scatterplot: Volume vs. Girth",
pch = 16, # filled circle points
col = "steelblue")
What to look for in the scatterplot:
- Form: Is the pattern linear (straight) or curved?
- Direction: Positive (both increase together) or negative (one increases as the other decreases)?
- Strength: Are the points tightly clustered around a line (strong) or widely scattered (weak)?
- Outliers: Are there any unusual observations far from the general pattern?
Step 4 — Compute the correlation coefficient
## [1] 0.9671194Interpretation: Report the value of \(r\) and describe its sign and magnitude using the reference table above. Running the code above gives \(r \approx 0.967\): “The sample correlation coefficient is \(r = 0.967\), indicating a very strong positive linear relationship between girth and volume.”
Step 5 — Hypothesis test for the population correlation
We test whether there is a statistically significant linear correlation in the population:
\[H_0: \rho = 0 \quad \text{(no linear correlation)} \qquad \text{vs.} \qquad H_a: \rho \neq 0\]
The test statistic is:
\[t = \frac{r\sqrt{n-2}}{\sqrt{1-r^2}} \sim t_{n-2} \text{ under } H_0\]
##
## Pearson's product-moment correlation
##
## data: x and y
## t = 20.478, df = 29, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9322519 0.9841887
## sample estimates:
## cor
## 0.9671194Interpretation of
cor.test()output:
Output field What it means corSample correlation coefficient \(r\) tTest statistic value dfDegrees of freedom (\(n - 2\)) p-valueReject \(H_0\) if p-value \(< \alpha = 0.05\) 95 percent confidence interval95% CI for the population correlation \(\rho\) 4-step conclusion template: 1. \(H_0: \rho = 0\) vs. \(H_a: \rho \neq 0\) 2. \(t \approx 20.5\), \(df = 29\) 3. \(p\text{-value} < 0.001\) 4. Since p-value \(< 0.05\), we reject \(H_0\). There is sufficient evidence of a significant linear correlation between girth and volume.
Practice II: Car Weight and Fuel Economy (Case Study B —
mtcars)
R’s built-in mtcars dataset records fuel consumption and
10 design/performance variables for 32 cars (1973–74 models). It will be
used as a running case study across several labs, so
it’s worth getting to know well now.
Variables:
wt(weight, 1000 lbs),mpg(miles per gallon),am(transmission: 0 = automatic, 1 = manual),cyl(number of cylinders: 4, 6, or 8),gear(number of forward gears: 3, 4, or 5). Research question: Is there a significant linear correlation between a car’s weight and its fuel economy?
Step 1 — Load the data
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...Step 2 — Explore the data
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...# Frequency table for transmission type (am: 0 = automatic, 1 = manual)
table1 <- table(mtcars$am)
table1##
## 0 1
## 19 13##
## 0 1
## 0.59375 0.40625# Two-way frequency table: number of cylinders by transmission type
table2 <- table(mtcars$cyl, mtcars$am)
table2##
## 0 1
## 4 3 8
## 6 4 3
## 8 12 2##
## 0 1
## 4 0.09375 0.25000
## 6 0.12500 0.09375
## 8 0.37500 0.06250## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.513 2.581 3.325 3.217 3.610 5.424## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.40 15.43 19.20 20.09 22.80 33.90# Boxplots
par(mfrow = c(1, 2))
boxplot(mtcars$mpg,
main = "Distribution of MPG", ylab = "Miles per Gallon")
boxplot(mtcars$mpg ~ mtcars$am,
xlab = "Transmission (0 = auto, 1 = manual)", ylab = "MPG",
main = "MPG by Transmission",
col = c("lightpink", "lightblue"))
par(mfrow = c(1, 1))
boxplot(mtcars$mpg ~ mtcars$gear,
xlab = "Number of Gears", ylab = "MPG",
main = "MPG by Number of Gears")
Interpretation: Boxplots let you compare distributions across groups. Look for differences in medians, spread (IQR), and the presence of outliers.
Step 3 — Scatterplot
x2 <- mtcars$wt
y2 <- mtcars$mpg
plot(x2, y2,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
main = "Scatterplot: MPG vs. Weight",
pch = 16, col = "steelblue")
Step 4 — Correlation coefficient and hypothesis test
## [1] -0.8676594# Hypothesis test: H0: rho = 0 vs Ha: rho ≠ 0
cor.test(x2, y2, alternative = "two.sided", conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: x2 and y2
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9338264 -0.7440872
## sample estimates:
## cor
## -0.8676594Running the code gives \(r \approx -0.868\) (\(p < 0.001\)) — a very strong negative linear relationship: heavier cars tend to have lower fuel economy.
Assignment: Employee Survey Ratings (Case Study C —
attitude)
R’s built-in attitude dataset comes from a survey of
clerical employees at a large financial organisation: 30 employees rated
their department (0–100 scale) on several dimensions, including overall
performance rating and perceived opportunity for
learning.
Variables:
learning(perceived learning opportunity score, 0–100) andrating(overall performance rating, 0–100).
Complete all questions using R. Show your R code and output, and write your interpretation in full sentences.
## 'data.frame': 30 obs. of 7 variables:
## $ rating : num 43 63 71 61 81 43 58 71 72 67 ...
## $ complaints: num 51 64 70 63 78 55 67 75 82 61 ...
## $ privileges: num 30 51 68 45 56 49 42 50 72 45 ...
## $ learning : num 39 54 69 47 66 44 56 55 67 47 ...
## $ raises : num 61 63 76 54 71 54 66 70 71 62 ...
## $ critical : num 92 73 86 84 83 49 68 66 83 80 ...
## $ advance : num 45 47 48 35 47 34 35 41 31 41 ...Perform descriptive statistics for both
learningandrating(mean, SD, min, max). Summarise the results.Make a scatterplot of
learning(\(x\)) vs.rating(\(y\)). Describe the form, direction, strength, and any outliers.Find the sample correlation coefficient \(r\) between
learningandrating. Interpret its value.Test whether the population correlation is different from zero at \(\alpha = 0.05\). Use the 4-step procedure: (i) state hypotheses, (ii) give the test statistic and degrees of freedom, (iii) give the p-value, (iv) state your conclusion in context.
Find the 95% confidence interval for the population correlation \(\rho\) and interpret it.
Lab 2: Simple Linear Regression
Objectives: Students are able to use R to:
- Fit a simple linear regression model using least squares
- Interpret the regression coefficients
- Perform inference on regression parameters (t-tests and confidence intervals)
- Assess model fit using \(r^2\) and \(MS_{res}\)
- Check regression assumptions using residual diagnostics
- Construct confidence and prediction intervals
Key concepts: The simple linear regression model expresses the relationship between a response variable \(Y\) and a single predictor \(X\) as:
\[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \qquad \varepsilon_i \overset{iid}{\sim} N(0, \sigma^2)\]
where:
- \(\beta_0\) = intercept (value of \(E(Y)\) when \(X = 0\))
- \(\beta_1\) = slope (change in \(E(Y)\) for a 1-unit increase in \(X\))
- \(\varepsilon_i\) = random error term (independent, normally distributed with constant variance)
The least-squares estimates \(\hat{\beta}_0\) and \(\hat{\beta}_1\) minimise \(\sum_{i=1}^n (y_i - \hat{y}_i)^2\).
Practice I: Cherry Tree Volume (continued)
We continue with the same trees data used in Lab 1.
Variables:
Girth(in.) as \(X\),Volume(cu. ft.) as \(Y\). Research question: Can we predict timber volume from trunk girth using a linear model?
Step 1 — Load data
Step 2 — Fit the regression model
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.065 -3.107 0.152 3.495 9.587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.9435 3.3651 -10.98 7.62e-12 ***
## x 5.0659 0.2474 20.48 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9331
## F-statistic: 419.4 on 1 and 29 DF, p-value: < 2.2e-16## 2.5 % 97.5 %
## (Intercept) -43.825953 -30.060965
## x 4.559914 5.571799Key output from
summary(fit):
Output field What it means Estimate(Intercept)\(\hat{\beta}_0\): predicted \(Y\) when \(X = 0\) Estimate(x)\(\hat{\beta}_1\): change in \(\hat{Y}\) for each 1-unit increase in \(X\) Std. ErrorStandard error of the coefficient estimate t valueTest statistic for \(H_0: \beta_j = 0\) Pr(>|t|)p-value; reject \(H_0\) if \(< 0.05\) Multiple R-squared\(r^2\): proportion of variance in \(Y\) explained by \(X\) F-statisticOverall model F-test (equivalent to t-test in simple regression) Residual standard error\(\hat{\sigma} = \sqrt{MS_{res}}\): estimated SD of errors How to write the fitted equation: Find the
Estimatecolumn in the output and substitute: \[\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x\] Here, \(\hat{\beta}_0 \approx -36.94\) and \(\hat{\beta}_1 \approx 5.07\), so \(\hat{y} = -36.94 + 5.07x\).
Step 3 — Plot the fitted regression line
plot(x, y,
xlab = "Girth (in.)",
ylab = "Volume (cu. ft.)",
main = "Regression Line: Volume vs. Girth",
pch = 16, col = "steelblue")
abline(fit, col = "red", lwd = 2) # overlay regression line
Step 4 — Residual analysis
Residuals are defined as \(e_i = y_i - \hat{y}_i\). Before trusting regression inference, we must verify three assumptions:
- Linearity — no systematic pattern in residuals vs. fitted values
- Normality — residuals are approximately normally distributed
- Constant variance (homoscedasticity) — spread of residuals is roughly equal for all fitted values
res <- resid(fit)
par(mfrow = c(1, 3))
# 1. Residuals vs. Fitted — check linearity and constant variance
plot(fitted(fit), res,
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals vs. Fitted",
pch = 16, col = "steelblue")
abline(h = 0, lty = 2, col = "red")
# 2. Histogram — check approximate normality
hist(res, main = "Histogram of Residuals",
xlab = "Residuals", col = "lightblue", border = "white")
# 3. Normal Q-Q plot — check normality more carefully
qqnorm(res, main = "Normal Q-Q Plot", pch = 16)
qqline(res, col = "red")
What to look for:
- Residuals vs. Fitted: points should scatter randomly around the zero line. A curve suggests non-linearity; a fan shape suggests non-constant variance.
- Histogram: should be roughly bell-shaped and symmetric around zero.
- Q-Q plot: points should fall close to the red diagonal line. Systematic deviations indicate non-normality.
Note: With girth and volume, a slight curve is visible in the residuals vs. fitted plot — a cylinder is only an approximation for a tree trunk’s shape, and real trees taper. This is a nice discussion point: sometimes real data hint that a straight line isn’t the perfect model.
Step 5 — Prediction for a new value of x
Two types of intervals are available:
- Confidence interval (CI): estimates the mean response \(E(Y \mid X = x^*)\) — how high is the average volume for all trees with girth \(= x^*\)?
- Prediction interval (PI): estimates a single future observation at \(X = x^*\) — what volume might an individual tree with girth \(= x^*\) have?
The PI is always wider than the CI because it includes both estimation uncertainty and individual variability.
## 1
## 28.91267# 95% Confidence interval for E(Y | x = 13)
predict(fit, newdata = new, interval = "confidence", level = 0.95)## fit lwr upr
## 1 28.91267 27.34573 30.47962# 95% Prediction interval for a new individual at x = 13
predict(fit, newdata = new, interval = "prediction", level = 0.95)## fit lwr upr
## 1 28.91267 20.07634 37.74901Output columns:
fit= point estimate,lwr= lower bound,upr= upper bound.
Step 6 — Plot fitted line with CI and PI bands
new_range <- data.frame(x = sort(x)) # sorted x values for smooth curves
pred.CI <- predict(fit, newdata = new_range, interval = "confidence")
pred.PI <- predict(fit, newdata = new_range, interval = "prediction")
plot(x, y,
xlab = "Girth (in.)", ylab = "Volume (cu. ft.)",
main = "Fitted Line with 95% CI and PI Bands",
pch = 16, col = "steelblue")
# Fitted regression line
lines(new_range$x, pred.CI[, "fit"], col = "black", lwd = 2)
# 95% Confidence interval (red dashed)
lines(new_range$x, pred.CI[, "lwr"], col = "red", lwd = 1.5, lty = 2)
lines(new_range$x, pred.CI[, "upr"], col = "red", lwd = 1.5, lty = 2)
# 95% Prediction interval (blue dotted)
lines(new_range$x, pred.PI[, "lwr"], col = "blue", lwd = 1.5, lty = 3)
lines(new_range$x, pred.PI[, "upr"], col = "blue", lwd = 1.5, lty = 3)
legend("topleft",
legend = c("Fitted line", "95% CI", "95% PI"),
col = c("black", "red", "blue"),
lty = c(1, 2, 3), lwd = 2, bty = "n")
Why is the PI always wider than the CI? The CI captures only the uncertainty in estimating the mean of \(Y\), while the PI must additionally account for the natural scatter of individual observations around that mean. The gap between CI and PI bands is largest near the extremes of \(x\) and narrowest near \(\bar{x}\).
Practice II: Car Weight and Fuel Economy (continued)
We continue with the same mtcars data used in Lab 1
Practice II.
Variables:
wt(1000 lbs) as \(X\),mpgas \(Y\). Research question: Can we predict fuel economy from weight using a linear model?
Step 1 — Load data
Step 2 — Fit the model
##
## Call:
## lm(formula = y2 ~ x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## x2 -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10## 2.5 % 97.5 %
## (Intercept) 33.450500 41.119753
## x2 -6.486308 -4.202635Step 3 — Plot fitted line
plot(x2, y2,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
main = "Regression Line: MPG vs. Weight",
pch = 16, col = "steelblue")
abline(fit2, col = "red", lwd = 2)
Step 4 — Residual analysis
res2 <- resid(fit2)
par(mfrow = c(1, 3))
plot(fitted(fit2), res2,
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals vs. Fitted", pch = 16, col = "steelblue")
abline(h = 0, lty = 2, col = "red")
hist(res2, main = "Histogram of Residuals",
xlab = "Residuals", col = "lightblue", border = "white")
qqnorm(res2, main = "Normal Q-Q Plot", pch = 16)
qqline(res2, col = "red")
Step 5 — CI and PI bands
new2 <- data.frame(x2 = sort(x2))
pred.CI2 <- predict(fit2, newdata = new2, interval = "confidence")
pred.PI2 <- predict(fit2, newdata = new2, interval = "prediction")
plot(x2, y2,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
main = "Fitted Line with 95% CI and PI Bands",
pch = 16, col = "steelblue",
ylim = range(c(pred.PI2[, "lwr"], pred.PI2[, "upr"])))
lines(new2$x2, pred.CI2[, "fit"], col = "black", lwd = 2)
lines(new2$x2, pred.CI2[, "lwr"], col = "red", lty = 2, lwd = 1.5)
lines(new2$x2, pred.CI2[, "upr"], col = "red", lty = 2, lwd = 1.5)
lines(new2$x2, pred.PI2[, "lwr"], col = "blue", lty = 3, lwd = 1.5)
lines(new2$x2, pred.PI2[, "upr"], col = "blue", lty = 3, lwd = 1.5)
legend("topleft",
legend = c("Fitted line", "95% CI", "95% PI"),
col = c("black", "red", "blue"),
lty = c(1, 2, 3), lwd = 2, bty = "n")
Assignment: Employee Survey Ratings (continued)
We continue with the attitude data from Lab 1. Now we
extend the analysis from correlation to regression.
Obtain the least-squares regression line for predicting
rating(\(Y\)) fromlearning(\(X\)). Write the fitted equation.Interpret \(\hat{\beta}_0\) and \(\hat{\beta}_1\) in the context of this problem.
Find \(MS_{res}\) (= residual standard error squared) and state its units.
Test whether
learningis a significant linear predictor ofratingat \(\alpha = 0.05\). Use the 4-step procedure.Find the 95% confidence intervals for \(\beta_0\) and \(\beta_1\).
Find \(r^2\) and interpret its value in context.
Check the regression assumptions using residual plots (residuals vs. fitted, histogram, Q-Q plot). Are the assumptions satisfied?
Find the 95% confidence interval and 95% prediction interval for the rating of a department with a learning score of 60. Explain the difference between the two intervals.
Lab 3: Multiple Linear Regression
Objectives: Students are able to use R to:
- Perform descriptive statistics for multiple linear regression
- Select independent variables into the model
- Perform multiple linear regression analysis and make inference on regression parameters
- Interpret the results
Key concepts: The multiple linear regression model with \(p\) predictors is:
\[Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + \varepsilon_i\]
We estimate the coefficients by least squares and evaluate the model with:
- Overall F-test: tests \(H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\)
- Individual t-tests: tests \(H_0: \beta_j = 0\) for each predictor
- \(R^2\) (multiple): proportion of variance in \(Y\) explained by all predictors together
- Adjusted \(R^2\): penalises for the number of predictors; better for comparing models of different sizes
Practice: Predicting Fuel Economy from Several Car Characteristics
We continue with mtcars. Rather than predicting
mpg from weight alone, we now bring in several quantitative
predictors at once:
| Variable | Description | Type |
|---|---|---|
disp |
Engine displacement (cu. in.) | Quantitative |
hp |
Gross horsepower | Quantitative |
drat |
Rear axle ratio | Quantitative |
wt |
Weight (1000 lbs) | Quantitative |
qsec |
Quarter-mile time (seconds) | Quantitative |
mpg |
Miles per gallon | Response |
Step 1 — Load data
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000Step 2 — Descriptive statistics and correlation matrix
Since all variables here are quantitative, we examine pairwise scatterplots and correlations.
## mpg disp hp drat wt qsec
## mpg 1.0000000 -0.8475514 -0.7761684 0.68117191 -0.8676594 0.41868403
## disp -0.8475514 1.0000000 0.7909486 -0.71021393 0.8879799 -0.43369788
## hp -0.7761684 0.7909486 1.0000000 -0.44875912 0.6587479 -0.70822339
## drat 0.6811719 -0.7102139 -0.4487591 1.00000000 -0.7124406 0.09120476
## wt -0.8676594 0.8879799 0.6587479 -0.71244065 1.0000000 -0.17471588
## qsec 0.4186840 -0.4336979 -0.7082234 0.09120476 -0.1747159 1.00000000# Pairwise scatterplot matrix
pairs(mtcars[, c("mpg", "disp", "hp", "drat", "wt", "qsec")],
main = "Scatterplot Matrix for mtcars",
pch = 16, col = "steelblue")
Interpretation: Look at the first row (or column) to see which predictors appear most linearly related to
mpg. Also check for high correlations between predictors (e.g.disp,hp, andwtare all strongly related to each other) — this is called multicollinearity and can cause instability in coefficient estimates.
Step 3 — Full model (all five predictors)
##
## Call:
## lm(formula = mpg ~ disp + hp + drat + wt + qsec, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5404 -1.6701 -0.4264 1.1320 5.4996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.53357 10.96423 1.508 0.14362
## disp 0.00872 0.01119 0.779 0.44281
## hp -0.02060 0.01528 -1.348 0.18936
## drat 2.01578 1.30946 1.539 0.13579
## wt -4.38546 1.24343 -3.527 0.00158 **
## qsec 0.64015 0.45934 1.394 0.17523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.558 on 26 degrees of freedom
## Multiple R-squared: 0.8489, Adjusted R-squared: 0.8199
## F-statistic: 29.22 on 5 and 26 DF, p-value: 6.892e-10## Analysis of Variance Table
##
## Response: mpg
## Df Sum Sq Mean Sq F value Pr(>F)
## disp 1 808.89 808.89 123.6185 2.23e-11 ***
## hp 1 33.67 33.67 5.1449 0.031854 *
## drat 1 30.15 30.15 4.6073 0.041340 *
## wt 1 70.51 70.51 10.7754 0.002933 **
## qsec 1 12.71 12.71 1.9422 0.175233
## Residuals 26 170.13 6.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Key outputs:
anova(full)gives the sequential (Type I) ANOVA table, showing the reduction in \(SS_{res}\) as each predictor is added.- The overall F-test p-value tests whether at least one \(\beta_j \neq 0\). Here it is highly significant (\(R^2 \approx 0.85\)), even though several individual predictors are not.
- Some individual predictors may be non-significant even if the overall model is — this is common with correlated predictors (see the multicollinearity note above).
Step 4 — Variable selection
With many predictors we use automatic selection to find the best subset. All three methods use AIC (Akaike Information Criterion) — lower AIC = better model.
null <- lm(mpg ~ 1, data = mtcars) # intercept-only (no predictors)
full <- lm(mpg ~ disp + hp + drat + wt + qsec, data = mtcars)
# Forward selection: starts empty, adds one predictor at a time
fw.fit <- step(null,
scope = list(lower = null, upper = full),
direction = "forward",
trace = 1)## Start: AIC=115.94
## mpg ~ 1
##
## Df Sum of Sq RSS AIC
## + wt 1 847.73 278.32 73.217
## + disp 1 808.89 317.16 77.397
## + hp 1 678.37 447.67 88.427
## + drat 1 522.48 603.57 97.988
## + qsec 1 197.39 928.66 111.776
## <none> 1126.05 115.943
##
## Step: AIC=73.22
## mpg ~ wt
##
## Df Sum of Sq RSS AIC
## + hp 1 83.274 195.05 63.840
## + qsec 1 82.858 195.46 63.908
## + disp 1 31.639 246.68 71.356
## <none> 278.32 73.217
## + drat 1 9.081 269.24 74.156
##
## Step: AIC=63.84
## mpg ~ wt + hp
##
## Df Sum of Sq RSS AIC
## <none> 195.05 63.840
## + drat 1 11.3659 183.68 63.919
## + qsec 1 8.9885 186.06 64.331
## + disp 1 0.0571 194.99 65.831# Backward elimination: starts full, removes one predictor at a time
be.fit <- step(full, direction = "backward", trace = 1)## Start: AIC=65.47
## mpg ~ disp + hp + drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## - disp 1 3.974 174.10 64.205
## <none> 170.13 65.466
## - hp 1 11.886 182.01 65.627
## - qsec 1 12.708 182.84 65.772
## - drat 1 15.506 185.63 66.258
## - wt 1 81.394 251.52 75.978
##
## Step: AIC=64.21
## mpg ~ hp + drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## - hp 1 9.418 183.52 63.891
## - qsec 1 9.578 183.68 63.919
## <none> 174.10 64.205
## - drat 1 11.956 186.06 64.331
## - wt 1 113.882 287.99 78.310
##
## Step: AIC=63.89
## mpg ~ drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## <none> 183.52 63.891
## - drat 1 11.942 195.46 63.908
## - qsec 1 85.720 269.24 74.156
## - wt 1 275.686 459.21 91.241# Stepwise regression: both directions (generally preferred)
sw.fit <- step(full, direction = "both", trace = 1)## Start: AIC=65.47
## mpg ~ disp + hp + drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## - disp 1 3.974 174.10 64.205
## <none> 170.13 65.466
## - hp 1 11.886 182.01 65.627
## - qsec 1 12.708 182.84 65.772
## - drat 1 15.506 185.63 66.258
## - wt 1 81.394 251.52 75.978
##
## Step: AIC=64.21
## mpg ~ hp + drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## - hp 1 9.418 183.52 63.891
## - qsec 1 9.578 183.68 63.919
## <none> 174.10 64.205
## - drat 1 11.956 186.06 64.331
## + disp 1 3.974 170.13 65.466
## - wt 1 113.882 287.99 78.310
##
## Step: AIC=63.89
## mpg ~ drat + wt + qsec
##
## Df Sum of Sq RSS AIC
## <none> 183.52 63.891
## - drat 1 11.942 195.46 63.908
## + hp 1 9.418 174.10 64.205
## + disp 1 1.506 182.02 65.627
## - qsec 1 85.720 269.24 74.156
## - wt 1 275.686 459.21 91.241Comparing the three methods: They do not always agree. In this dataset, forward selection tends to settle on
{wt, hp}, while backward elimination and stepwise (both) agree on{drat, wt, qsec}. Stepwise regression is generally preferred as it reconsiders variables at each step.
Step 5 — Fitted model from stepwise regression
##
## Call:
## lm(formula = mpg ~ drat + wt + qsec, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1152 -1.8273 -0.2696 1.0502 5.5010
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.3945 8.0689 1.412 0.16892
## drat 1.6561 1.2269 1.350 0.18789
## wt -4.3978 0.6781 -6.485 5.01e-07 ***
## qsec 0.9462 0.2616 3.616 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.56 on 28 degrees of freedom
## Multiple R-squared: 0.837, Adjusted R-squared: 0.8196
## F-statistic: 47.93 on 3 and 28 DF, p-value: 3.723e-11## Analysis of Variance Table
##
## Response: mpg
## Df Sum Sq Mean Sq F value Pr(>F)
## drat 1 522.48 522.48 79.715 1.105e-09 ***
## wt 1 334.33 334.33 51.008 9.009e-08 ***
## qsec 1 85.72 85.72 13.078 0.001163 **
## Residuals 28 183.52 6.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 2.5 % 97.5 %
## (Intercept) -5.1339149 27.922817
## drat -0.8571278 4.169418
## wt -5.7868161 -3.008776
## qsec 0.4102527 1.482155How to write the fitted equation: If stepwise selected
drat,wt, andqsec: \[\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \, \text{drat} + \hat{\beta}_2 \, \text{wt} + \hat{\beta}_3 \, \text{qsec}\] Substitute theEstimatevalues fromsummary()output (approximately \(\hat{y} = 11.39 + 1.66\,\text{drat} - 4.40\,\text{wt} + 0.95\,\text{qsec}\)).
Assignment: Predicting Fuel Economy (continued)
Using the results from the Practice section above, answer the following:
Which predictors were selected by each method (forward, backward, stepwise)? Do they agree?
Write the fitted regression equation from the stepwise model.
Test the overall fit at \(\alpha = 0.05\) using the 4-step hypothesis testing procedure.
Interpret each significant regression coefficient in the context of the problem.
Find the 95% CI for each coefficient in the selected model.
Find \(R^2\) and adjusted \(R^2\). What do they tell you?
Lab 4: Multiple Linear Regression with Dummy Variables
Objectives: Students are able to use R to:
- Recognise when and why dummy variables are needed
- Transform a qualitative variable into dummy variables
- Fit a multiple regression model with dummy variables
- Interpret all regression coefficients, including dummy variable coefficients
- Derive and compare the fitted equations for each group
Key concepts — Dummy variables: When a predictor is qualitative (categorical), we cannot enter it directly into a regression model as a number. Instead, we create \(k - 1\) dummy (indicator) variables for a variable with \(k\) groups. The omitted group is called the reference category.
Each dummy variable’s coefficient \(\hat{\beta}_j\) represents the estimated difference in the mean response between that group and the reference group, holding all other predictors constant.
Why \(k - 1\) dummies, not \(k\)? Using \(k\) dummies would create perfect multicollinearity (the dummy columns would always sum to 1, which equals the intercept column). Dropping one group — the reference — avoids this problem. The intercept then represents the mean response for the reference group.
General dummy coding scheme for a variable with \(k = 3\) groups (A = reference, B, C):
| Group | \(D_B\) | \(D_C\) |
|---|---|---|
| A (reference) | 0 | 0 |
| B | 1 | 0 |
| C | 0 | 1 |
The model becomes: \[\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 X_1 + \hat{\beta}_2 D_B + \hat{\beta}_3 D_C\]
This gives separate parallel regression lines for each group:
| Group | Fitted equation | Intercept |
|---|---|---|
| A (reference) | \(\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 X_1\) | \(\hat{\beta}_0\) |
| B | \(\hat{Y} = (\hat{\beta}_0 + \hat{\beta}_2) + \hat{\beta}_1 X_1\) | \(\hat{\beta}_0 + \hat{\beta}_2\) |
| C | \(\hat{Y} = (\hat{\beta}_0 + \hat{\beta}_3) + \hat{\beta}_1 X_1\) | \(\hat{\beta}_0 + \hat{\beta}_3\) |
All three lines share the same slope \(\hat{\beta}_1\) but have different intercepts.
Practice I: Fuel Economy by Weight and Number of Cylinders
Research question: Can fuel economy be predicted from weight and number of cylinders?
Variables:
| Variable | Name in data | Type |
|---|---|---|
| Fuel economy (mpg) | mpg |
Response (\(Y\)) |
| Weight (1000 lbs) | wt |
Quantitative predictor (\(X_1\)) |
| Number of cylinders | cyl |
Qualitative predictor (4 / 6 / 8) |
Step 1 — Load data
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...Step 2 — Descriptive statistics
Before fitting, explore the relationship between each predictor and the response separately.
y <- mtcars$mpg
x1 <- mtcars$wt # quantitative predictor
x2 <- factor(mtcars$cyl) # qualitative predictor (4, 6, or 8 cylinders)
# Summaries
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.40 15.43 19.20 20.09 22.80 33.90## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.513 2.581 3.325 3.217 3.610 5.424## x2
## 4 6 8
## 11 7 14## [1] -0.8676594# Scatterplot coloured by cylinder group
plot(x1, y,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
pch = 16, col = as.integer(x2),
main = "MPG vs. Weight, by Cylinders")
legend("topright",
legend = levels(x2),
col = 1:length(levels(x2)),
pch = 16, bty = "n")
# Boxplot: mpg by cylinder group
boxplot(y ~ x2,
xlab = "Number of Cylinders", ylab = "Miles per Gallon",
main = "MPG by Cylinders",
col = c("lightblue", "lightgreen", "lightyellow"))
What to look for:
- Do the groups (4, 6, 8 cylinders) appear to have different average fuel economy in the boxplot? If yes, cylinder count is likely a useful predictor.
- In the scatterplot, do the three groups of points appear to sit at different vertical levels? If so, they may have different intercepts.
- Is the relationship between weight and mpg roughly linear within each group?
# Alternative using ggplot2 (install once: install.packages("ggplot2"))
library(ggplot2)
# Scatterplot coloured by cylinder group
ggplot(mtcars, aes(x = wt, y = mpg, colour = factor(cyl))) +
geom_point(size = 2) +
labs(x = "Weight (1000 lbs)", y = "Miles per Gallon",
colour = "Cylinders",
title = "MPG vs. Weight, by Cylinders") +
theme_minimal()
Step 3 — Create dummy variables
With \(k = 3\) groups (4, 6, 8 cylinders), we need \(k - 1 = 2\) dummy variables. We choose 4 cylinders as the reference category.
x2.dummy <- ifelse(mtcars$cyl == 6, 1, 0) # D_Six
x3.dummy <- ifelse(mtcars$cyl == 8, 1, 0) # D_Eight
# Verify the coding is correct: print side-by-side
data.frame(Cylinders = mtcars$cyl,
D_Six = x2.dummy,
D_Eight = x3.dummy)## Cylinders D_Six D_Eight
## 1 6 1 0
## 2 6 1 0
## 3 4 0 0
## 4 6 1 0
## 5 8 0 1
## 6 6 1 0
## 7 8 0 1
## 8 4 0 0
## 9 4 0 0
## 10 6 1 0
## 11 6 1 0
## 12 8 0 1
## 13 8 0 1
## 14 8 0 1
## 15 8 0 1
## 16 8 0 1
## 17 8 0 1
## 18 4 0 0
## 19 4 0 0
## 20 4 0 0
## 21 4 0 0
## 22 8 0 1
## 23 8 0 1
## 24 8 0 1
## 25 8 0 1
## 26 4 0 0
## 27 4 0 0
## 28 4 0 0
## 29 8 0 1
## 30 6 1 0
## 31 8 0 1
## 32 4 0 0Dummy coding summary:
Cylinders \(D_{\text{Six}}\) \(D_{\text{Eight}}\) 4 (reference) 0 0 6 1 0 8 0 1 You can verify: when both dummies are 0, the observation must be a 4-cylinder car. ✓
Step 4 — Fit the full model
We fit the model with \(X_1\)
(wt), \(D_{\text{Six}}\),
and \(D_{\text{Eight}}\) as
predictors.
##
## Call:
## lm(formula = y ~ x1 + x2.dummy + x3.dummy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5890 -1.2357 -0.5159 1.3845 5.7915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9908 1.8878 18.006 < 2e-16 ***
## x1 -3.2056 0.7539 -4.252 0.000213 ***
## x2.dummy -4.2556 1.3861 -3.070 0.004718 **
## x3.dummy -6.0709 1.6523 -3.674 0.000999 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.557 on 28 degrees of freedom
## Multiple R-squared: 0.8374, Adjusted R-squared: 0.82
## F-statistic: 48.08 on 3 and 28 DF, p-value: 3.594e-11## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 847.73 847.73 129.6650 5.079e-12 ***
## x2.dummy 1 7.00 7.00 1.0713 0.3095128
## x3.dummy 1 88.26 88.26 13.4999 0.0009992 ***
## Residuals 28 183.06 6.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 2.5 % 97.5 %
## (Intercept) 30.123824 37.857764
## x1 -4.749898 -1.661328
## x2.dummy -7.094824 -1.416341
## x3.dummy -9.455418 -2.686301Reading the
summary()output:
(Intercept): \(\hat{\beta}_0 \approx 34.0\) — estimated mean mpg for the 4-cylinder group whenwt = 0(a mathematical extrapolation, not a realistic car).x1: \(\hat{\beta}_1 \approx -3.21\) — estimated change in mpg for each additional 1000 lbs, within any cylinder group (slope is shared).x2.dummy: \(\hat{\beta}_2 \approx -4.26\) — estimated difference in mean mpg between 6-cylinder and 4-cylinder cars, holding weight constant.x3.dummy: \(\hat{\beta}_3 \approx -6.07\) — estimated difference in mean mpg between 8-cylinder and 4-cylinder cars, holding weight constant.All three coefficients are statistically significant (\(p < 0.01\)). Both dummy coefficients are negative, meaning that — at the same weight — cars with more cylinders get worse fuel economy.
Step 5 — Write the fitted equation and derive group-specific equations
From the Estimate column of
summary(fit.full), substitute values into:
\[\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 D_{\text{Six}} + \hat{\beta}_3 D_{\text{Eight}}\]
Then derive the three separate equations by substituting the dummy values:
## (Intercept) x1 x2.dummy x3.dummy
## 33.990794 -3.205613 -4.255582 -6.070860# The three group equations are derived as follows:
# 4-cylinder (D_Six=0, D_Eight=0):
# yhat = b0 + b1*x1
# 6-cylinder (D_Six=1, D_Eight=0):
# yhat = (b0 + b2) + b1*x1
# 8-cylinder (D_Six=0, D_Eight=1):
# yhat = (b0 + b3) + b1*x1Interpretation example (from the actual fitted model): \(\hat{\beta}_0 \approx 34.0\), \(\hat{\beta}_1 \approx -3.21\), \(\hat{\beta}_2 \approx -4.26\), \(\hat{\beta}_3 \approx -6.07\).
- 4-cylinder: \(\hat{y} = 34.0 - 3.21\,x_1\)
- 6-cylinder: \(\hat{y} = 29.7 - 3.21\,x_1\) (intercept decreases by 4.26)
- 8-cylinder: \(\hat{y} = 27.9 - 3.21\,x_1\) (intercept decreases by 6.07)
Interpretation: “At the same weight, 8-cylinder cars are estimated to get about 6.07 fewer miles per gallon on average than 4-cylinder cars.”
Step 6 — Plot the three fitted lines
# Compute fitted values
yhat <- fitted(fit.full)
# Get coefficient estimates
b <- coef(fit.full)
b0 <- b["(Intercept)"]
b1 <- b["x1"]
b2 <- b["x2.dummy"]
b3 <- b["x3.dummy"]
x1_range <- seq(min(x1), max(x1), length.out = 100)
# Base scatterplot coloured by cylinder group
plot(x1, y,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
pch = 16, col = as.integer(x2),
main = "Fitted Regression Lines by Cylinder Group")
# Add a fitted line for each group
lines(x1_range, b0 + b1 * x1_range, col = 1, lwd = 2) # 4-cylinder
lines(x1_range, b0 + b2 + b1 * x1_range, col = 2, lwd = 2) # 6-cylinder
lines(x1_range, b0 + b3 + b1 * x1_range, col = 3, lwd = 2) # 8-cylinder
legend("topright",
legend = c("4-cylinder", "6-cylinder", "8-cylinder"),
col = 1:3, lty = 1, lwd = 2, pch = 16, bty = "n")
What this plot shows: Three parallel lines — same slope, different intercepts. The vertical gap between any two lines equals the corresponding dummy coefficient. If the lines were not parallel, that would suggest an interaction between cylinders and weight (a more advanced model not covered here).
Step 7 — Plot observed vs. predicted
plot(x1, y,
xlab = "Weight (1000 lbs)", ylab = "Miles per Gallon",
pch = 1, main = "Observed vs. Predicted MPG")
points(x1, yhat, pch = 20, col = "red")
legend("topright",
legend = c("Observed y", "Predicted ŷ"),
pch = c(1, 20), col = c("black", "red"), bty = "n")
Assignment: Student Performance — Part 1 (Dummy Coding and Interpretation)
For this assignment we use a simulated dataset of 48
students. It is generated directly in R with set.seed(), so
it is fully reproducible and requires no download.
Variables: Y (post-test score),
Pretest (pre-test score), EQ (emotional
quotient score), Method (teaching method: A, B, C, D).
set.seed(208251)
n <- 48
Method <- factor(rep(c("A", "B", "C", "D"), each = n / 4))
Pretest <- round(rnorm(n, mean = 60, sd = 8), 1)
EQ <- round(rnorm(n, mean = 70, sd = 7), 1)
# True (unknown to students) method effects, added as noise for realism
method_effect <- c(A = 0, B = 4, C = 9, D = 6)[as.character(Method)]
Y <- round(15 + 0.55 * Pretest + 0.30 * EQ + method_effect +
rnorm(n, mean = 0, sd = 5), 1)
student_data <- data.frame(Y, Pretest, EQ, Method)
head(student_data)## Y Pretest EQ Method
## 1 82.9 71.1 74.2 A
## 2 74.9 65.2 78.2 A
## 3 69.3 63.3 72.2 A
## 4 73.9 53.8 77.3 A
## 5 71.3 68.4 60.5 A
## 6 71.0 66.0 70.9 A## 'data.frame': 48 obs. of 4 variables:
## $ Y : num 82.9 74.9 69.3 73.9 71.3 71 50.3 64.7 61.1 65.6 ...
## $ Pretest: num 71.1 65.2 63.3 53.8 68.4 66 47.1 57.9 51.4 52.9 ...
## $ EQ : num 74.2 78.2 72.2 77.3 60.5 70.9 58.9 65.2 70 67.7 ...
## $ Method : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 1 1 1 1 ...Explore the data:
- Create scatterplots of \(Y\)
vs.
Pretestand \(Y\) vs.EQ. Describe each relationship (form, direction, strength). - Create a boxplot of \(Y\) by
Method. Do the teaching methods appear to differ in mean score?
- Create scatterplots of \(Y\)
vs.
Methodis a qualitative variable with \(k = 4\) groups. Create \(k - 1 = 3\) dummy variables. Clearly state:- Which group is the reference category.
- The complete dummy coding table.
- The R code used to create the dummies and verify them.
Fit the full model: \(Y = \beta_0 + \beta_1 \cdot \text{Pretest} + \beta_2 \cdot \text{EQ} + \beta_3 D_B + \beta_4 D_C + \beta_5 D_D + \varepsilon\). Show the R output and write the fitted equation.
Interpret each coefficient in the context of the problem. Pay particular attention to the dummy coefficients — what do they tell you about the effect of teaching method?
Write the separate fitted equation for students taught by each of Method A, B, C, and D. Comment on how the methods compare.
Lab 5: Variable Selection and Model Diagnostics
Objectives: Students are able to use R to:
- Apply automatic variable selection methods (forward, backward, stepwise)
- Understand and compare AIC-based selection criteria
- Perform a comprehensive model diagnostic analysis
- Interpret diagnostic plots and formal tests for regression assumptions
- Identify influential observations and outliers
- Make predictions using the selected model
Key concepts — Variable selection: When we have many candidate predictors, automatic selection methods help find a parsimonious model (few predictors, good fit). All three methods below use AIC (Akaike Information Criterion): lower AIC = better trade-off between fit and complexity.
| Method | Starting model | Strategy |
|---|---|---|
| Forward selection | Intercept only | Add the predictor that reduces AIC the most, repeat |
| Backward elimination | Full model | Remove the predictor that reduces AIC the most, repeat |
| Stepwise (both) | Full model | Add or remove at each step — generally preferred |
Important limitation: Automated selection does not guarantee the “best” model in any absolute sense. Always check whether the selected model makes practical sense and satisfies the regression assumptions.
Key concepts — Model diagnostics: Before using a fitted model for inference or prediction, we must verify four assumptions of the regression model \(Y_i = \beta_0 + \beta_1 X_{i1} + \cdots + \varepsilon_i\):
| Assumption | What it means | How to check |
|---|---|---|
| Linearity | \(E(\varepsilon_i) = 0\); no pattern in residuals | Residuals vs. Fitted plot |
| Normality | \(\varepsilon_i \sim N(0, \sigma^2)\) | Histogram, Q-Q plot, Shapiro-Wilk test |
| Constant variance | \(\text{Var}(\varepsilon_i) = \sigma^2\) (homoscedasticity) | Scale-Location plot, Breusch-Pagan test |
| Independence | Residuals are uncorrelated | Residuals vs. order plot, Durbin-Watson test |
Additionally, we check for influential observations — data points that have an unusually large effect on the fitted model.
Practice I: Fuel Economy by Weight and Cylinders (continued)
We continue with the mtcars dummy-variable model from
Lab 4. The goal is now to select the best subset of predictors and
validate the model assumptions.
Step 1 — Re-load data and recreate variables
data(mtcars)
y <- mtcars$mpg
x1 <- mtcars$wt
x2.dummy <- ifelse(mtcars$cyl == 6, 1, 0)
x3.dummy <- ifelse(mtcars$cyl == 8, 1, 0)Step 2 — Variable selection
null <- lm(y ~ 1) # intercept-only model
full <- lm(y ~ x1 + x2.dummy + x3.dummy) # full model from Lab 4
# Forward selection
fw.fit <- step(null,
scope = list(lower = null, upper = full),
direction = "forward",
trace = 1)## Start: AIC=115.94
## y ~ 1
##
## Df Sum of Sq RSS AIC
## + x1 1 847.73 278.32 73.217
## + x3.dummy 1 619.89 506.16 92.355
## <none> 1126.05 115.943
## + x2.dummy 1 1.08 1124.96 117.913
##
## Step: AIC=73.22
## y ~ x1
##
## Df Sum of Sq RSS AIC
## + x3.dummy 1 33.635 244.69 71.096
## <none> 278.32 73.217
## + x2.dummy 1 7.004 271.32 74.402
##
## Step: AIC=71.1
## y ~ x1 + x3.dummy
##
## Df Sum of Sq RSS AIC
## + x2.dummy 1 61.628 183.06 63.810
## <none> 244.69 71.096
##
## Step: AIC=63.81
## y ~ x1 + x3.dummy + x2.dummy## Start: AIC=63.81
## y ~ x1 + x2.dummy + x3.dummy
##
## Df Sum of Sq RSS AIC
## <none> 183.06 63.810
## - x2.dummy 1 61.628 244.69 71.096
## - x3.dummy 1 88.259 271.32 74.402
## - x1 1 118.204 301.26 77.752# Stepwise regression (both directions) — generally preferred
sw.fit <- step(full,
direction = "both",
trace = 1)## Start: AIC=63.81
## y ~ x1 + x2.dummy + x3.dummy
##
## Df Sum of Sq RSS AIC
## <none> 183.06 63.810
## - x2.dummy 1 61.628 244.69 71.096
## - x3.dummy 1 88.259 271.32 74.402
## - x1 1 118.204 301.26 77.752Reading the
step()output: Each row shows one candidate model. The columns are:
Df: degrees of freedom used by the variable being added/removedAIC: AIC of the candidate model after the changeThe algorithm selects the change that gives the lowest AIC. It stops when no change improves AIC.
Comparing the three methods: In this dataset, all three methods agree that all three predictors (
x1,x2.dummy,x3.dummy) should be retained — the full model already has the lowest AIC. This won’t always happen, but it’s a useful reminder that selection sometimes confirms your original model rather than simplifying it.
Step 3 — Examine the selected model
##
## Call:
## lm(formula = y ~ x1 + x2.dummy + x3.dummy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5890 -1.2357 -0.5159 1.3845 5.7915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9908 1.8878 18.006 < 2e-16 ***
## x1 -3.2056 0.7539 -4.252 0.000213 ***
## x2.dummy -4.2556 1.3861 -3.070 0.004718 **
## x3.dummy -6.0709 1.6523 -3.674 0.000999 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.557 on 28 degrees of freedom
## Multiple R-squared: 0.8374, Adjusted R-squared: 0.82
## F-statistic: 48.08 on 3 and 28 DF, p-value: 3.594e-11## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 847.73 847.73 129.6650 5.079e-12 ***
## x2.dummy 1 7.00 7.00 1.0713 0.3095128
## x3.dummy 1 88.26 88.26 13.4999 0.0009992 ***
## Residuals 28 183.06 6.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 2.5 % 97.5 %
## (Intercept) 30.123824 37.857764
## x1 -4.749898 -1.661328
## x2.dummy -7.094824 -1.416341
## x3.dummy -9.455418 -2.686301Overall F-test (4-step procedure):
- \(H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\) (no predictors are useful) vs. \(H_a\): at least one \(\beta_j \neq 0\)
- \(F = MS_{reg} / MS_{res}\), with \(df_1 = p\) and \(df_2 = n - p - 1\)
- p-value from the
F-statisticline insummary()- Conclusion in context
Step 4 — Comprehensive model diagnostics
R’s built-in plot() for lm objects produces
four diagnostic plots at once. Understanding each one is essential.
Interpreting the four diagnostic plots:
Plot 1 — Residuals vs. Fitted - Checks: linearity and constant variance - Good: points scatter randomly around the horizontal zero line, with roughly equal spread across all fitted values - Warning signs: a curve (non-linearity) or a funnel/fan shape (non-constant variance)
Plot 2 — Normal Q-Q - Checks: normality of residuals - Good: points follow the diagonal reference line closely - Warning signs: S-shaped curve (heavy tails), points bending away at the ends
Plot 3 — Scale-Location (√|standardised residuals| vs. Fitted) - Checks: constant variance (more sensitive than Plot 1) - Good: a roughly horizontal red line with points evenly spread - Warning signs: a clear upward or downward trend in the red line
Plot 4 — Residuals vs. Leverage - Checks: influential observations - Points with high leverage AND large residuals are most dangerous - Points outside the Cook’s distance dashed lines (if visible) are highly influential
Step 5 — Formal tests for each assumption
res <- resid(sw.fit)
# --- Normality ---
# Visual: histogram and Q-Q plot side by side
par(mfrow = c(1, 2))
hist(res, main = "Histogram of Residuals",
xlab = "Residuals", col = "lightblue", border = "white",
freq = FALSE) # density scale
curve(dnorm(x, mean = mean(res), sd = sd(res)),
add = TRUE, col = "red", lwd = 2) # overlay normal curve
qqnorm(res, main = "Normal Q-Q Plot", pch = 16, col = "steelblue")
qqline(res, col = "red", lwd = 2)
par(mfrow = c(1, 1))
# Shapiro-Wilk test (best for n < 50; use for up to ~2000 obs)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.95907, p-value = 0.259# H0: residuals are normally distributed
# Satisfied if p-value > 0.05 (here p ≈ 0.26 — satisfied)
# --- Constant variance (homoscedasticity) ---
# install.packages("lmtest") # run once if not yet installed
library(lmtest)
bptest(sw.fit)##
## studentized Breusch-Pagan test
##
## data: sw.fit
## BP = 6.6551, df = 3, p-value = 0.08374# H0: variance is constant (homoscedastic)
# Satisfied if p-value > 0.05 (here p ≈ 0.08 — satisfied, though borderline)
# --- Independence (no autocorrelation) ---
dwtest(sw.fit)##
## Durbin-Watson test
##
## data: sw.fit
## DW = 1.8064, p-value = 0.2188
## alternative hypothesis: true autocorrelation is greater than 0# H0: residuals are not autocorrelated
# Satisfied if p-value > 0.05 (DW statistic close to 2 is ideal; here DW ≈ 1.81)Summary of formal tests:
Test \(H_0\) Assumption satisfied if Shapiro-Wilk Residuals are normally distributed p-value > 0.05 Breusch-Pagan Variance is constant (homoscedastic) p-value > 0.05 Durbin-Watson No autocorrelation in residuals p-value > 0.05 If an assumption is violated, consider:
Violation Possible remedy Non-linearity Add a polynomial term (e.g. \(X^2\)) or transform \(X\) Non-constant variance Apply log or square-root transformation to \(Y\) Non-normality Transform \(Y\); check for outliers Autocorrelation Include a time trend; use time-series methods
Step 6 — Checking for influential observations
An observation can be unusual in two ways:
- Outlier: unusual response value (large residual)
- High-leverage point: unusual predictor value (far from the centre of the \(X\) space)
- Influential point: changes the fitted model substantially if removed (high Cook’s distance)
# Standardised residuals — values beyond ±2 are potential outliers
std.res <- rstandard(sw.fit)
which(abs(std.res) > 2) # row indices of potential outliers (e.g. Fiat 128, Toyota Corolla)## 18 20
## 18 20# Leverage values — rule of thumb: high if > 2*(p+1)/n
n <- nrow(mtcars)
p <- length(coef(sw.fit)) - 1 # number of predictors
hat.vals <- hatvalues(sw.fit)
which(hat.vals > 2 * (p + 1) / n)## named integer(0)# Cook's distance — rule of thumb: influential if > 4/n
cooks <- cooks.distance(sw.fit)
which(cooks > 4 / n) # e.g. Chrysler Imperial, Fiat 128, Toyota Corolla## 17 18 20
## 17 18 20# Combined influence plot
plot(sw.fit, which = 4,
main = "Cook's Distance",
sub = paste("Threshold = 4/n =", round(4/n, 4)))
abline(h = 4/n, col = "red", lty = 2)
What to do with influential observations: Do NOT automatically delete them. First: 1. Check whether the data point was recorded correctly (data entry error?). 2. Fit the model with and without the point and compare — if results change substantially, report both. 3. If justified, remove the point and note it in your report.
Step 7 — Prediction using the selected model
# Predict mpg for a new car:
# wt = 3.0 (thousand lbs), cyl = 6 (D_Six=1, D_Eight=0)
new_obs <- data.frame(x1 = 3.0,
x2.dummy = 1,
x3.dummy = 0)
# Point estimate
predict(sw.fit, newdata = new_obs)## 1
## 20.11837# 95% confidence interval for the mean mpg at these values
predict(sw.fit, newdata = new_obs, interval = "confidence", level = 0.95)## fit lwr upr
## 1 20.11837 18.1305 22.10625# 95% prediction interval for a single new car
predict(sw.fit, newdata = new_obs, interval = "prediction", level = 0.95)## fit lwr upr
## 1 20.11837 14.51622 25.72052Assignment: Student Performance — Part 2 (Variable Selection and Diagnostics)
This assignment continues directly from Lab 4 Assignment. Re-generate
the same simulated dataset (the set.seed() guarantees you
get identical data) and use the dummy variables you created there.
set.seed(208251)
n <- 48
Method <- factor(rep(c("A", "B", "C", "D"), each = n / 4))
Pretest <- round(rnorm(n, mean = 60, sd = 8), 1)
EQ <- round(rnorm(n, mean = 70, sd = 7), 1)
method_effect <- c(A = 0, B = 4, C = 9, D = 6)[as.character(Method)]
Y <- round(15 + 0.55 * Pretest + 0.30 * EQ + method_effect +
rnorm(n, mean = 0, sd = 5), 1)
student_data <- data.frame(Y, Pretest, EQ, Method)Apply forward selection, backward elimination, and stepwise regression to select the best predictors of \(Y\) from
{Pretest, EQ, D_B, D_C, D_D}. Compare the three methods — do they agree?Write the fitted equation from the stepwise-selected model. Interpret each coefficient in context.
Test the overall fit of the selected model at \(\alpha = 0.05\) using the 4-step procedure. Report \(R^2\) and adjusted \(R^2\) and explain what they tell you.
Produce R’s four built-in diagnostic plots (
plot(model)). For each plot, describe what you see and state what assumption it is checking.Carry out the three formal assumption tests:
- Shapiro-Wilk test for normality
- Breusch-Pagan test for constant variance
- Durbin-Watson test for independence
For each test: state \(H_0\), report the test statistic and p-value, and conclude whether the assumption is satisfied.
Check for influential observations using Cook’s distance. Report any observations that exceed the threshold \(4/n\) and describe what you would do with them.
Lab 6: Non-Parametric Statistics I
Objectives: Students are able to:
- Perform descriptive statistics
- Apply appropriate non-parametric tests to answer research questions
When to use non-parametric tests: Use non-parametric methods when:
- The sample size is small and normality cannot be assumed
- Data are ordinal (ranked) rather than continuous
- Assumptions of the corresponding parametric test are clearly violated
| Parametric test | Non-parametric equivalent | Use case |
|---|---|---|
| One-sample t-test | Wilcoxon signed-rank test | One sample, test about median |
| Independent samples t-test | Mann-Whitney U (Wilcoxon rank-sum) | Two independent samples |
| Paired t-test | Wilcoxon signed-rank test (paired) | Two related samples |
| One-way ANOVA | Kruskal-Wallis test | Three or more independent samples |
| Repeated-measures ANOVA | Friedman test | Three or more related samples |
| Pearson correlation | Spearman rank correlation | Association, non-normal data |
Practice I: Fuel Economy and Car Characteristics
Variables:
mpg,qsec,hp,wt,am(transmission),cyl(cylinders)
Load and prepare data
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...# Convert qualitative variables to factors
mtcars$am <- factor(mtcars$am, labels = c("automatic", "manual"))
mtcars$cyl <- factor(mtcars$cyl)Question A: Is the median mpg equal to 20?
One sample, test about a median → Wilcoxon signed-rank test
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.40 15.43 19.20 20.09 22.80 33.90## [1] 6.026948par(mfrow = c(1, 2))
boxplot(mtcars$mpg, main = "Boxplot of MPG", ylab = "MPG")
hist(mtcars$mpg, main = "Histogram of MPG",
xlab = "MPG", col = "lightblue", border = "white")
par(mfrow = c(1, 1))
# H0: median mpg = 20 Ha: median mpg ≠ 20
wilcox.test(x = mtcars$mpg, mu = 20, alternative = "two.sided")##
## Wilcoxon signed rank exact test
##
## data: mtcars$mpg
## V = 249, p-value = 0.7855
## alternative hypothesis: true location is not equal to 20Interpretation: The Wilcoxon signed-rank test tests \(H_0\): population median =
mu. - p-value < 0.05: reject \(H_0\); the median mpg is significantly different from 20. - p-value ≥ 0.05: insufficient evidence to conclude the median differs from 20.
Question B: Is the median quarter-mile time
(qsec) less than 18 seconds?
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.50 16.89 17.71 17.85 18.90 22.90par(mfrow = c(1, 2))
boxplot(mtcars$qsec, ylab = "Quarter-mile time (sec)", main = "Boxplot")
hist(mtcars$qsec, xlab = "Quarter-mile time (sec)",
main = "Histogram", col = "lightblue", border = "white")
par(mfrow = c(1, 1))
# H0: median = 18 Ha: median < 18
wilcox.test(x = mtcars$qsec, mu = 18, alternative = "less")##
## Wilcoxon signed rank exact test
##
## data: mtcars$qsec
## V = 225.5, p-value = 0.2427
## alternative hypothesis: true location is less than 18Practice II: Comparing Three Timing Methods (Case Study E —
RoundingTimes)
This example uses a classic teaching dataset for the Friedman test,
distributed with R’s own documentation (?friedman.test). It
records how long it took 22 different operators to
round the corner of a piece using three measurement
methods: “Round Out”, “Narrow Angle”, and “Wide Angle”. Because
each operator used all three methods, this is a paired (repeated
measures) design.
Step 1 — Enter the data directly
RoundingTimes <-
matrix(c(5.40, 5.50, 5.55,
5.85, 5.70, 5.75,
5.20, 5.60, 5.50,
5.55, 5.50, 5.40,
5.90, 5.85, 5.70,
5.45, 5.55, 5.60,
5.40, 5.40, 5.35,
5.45, 5.50, 5.35,
5.25, 5.15, 5.00,
5.85, 5.80, 5.70,
5.25, 5.20, 5.10,
5.65, 5.55, 5.45,
5.60, 5.35, 5.45,
5.05, 5.00, 4.95,
5.50, 5.50, 5.40,
5.45, 5.55, 5.50,
5.55, 5.55, 5.35,
5.45, 5.50, 5.55,
5.50, 5.45, 5.25,
5.65, 5.60, 5.40,
5.70, 5.65, 5.55,
6.30, 6.30, 6.25),
nrow = 22,
byrow = TRUE,
dimnames = list(1 : 22,
c("RoundOut", "NarrowAngle", "WideAngle")))
rt <- as.data.frame(RoundingTimes)
head(rt)## RoundOut NarrowAngle WideAngle
## 1 5.40 5.50 5.55
## 2 5.85 5.70 5.75
## 3 5.20 5.60 5.50
## 4 5.55 5.50 5.40
## 5 5.90 5.85 5.70
## 6 5.45 5.55 5.60Step 2 — Explore data
## RoundOut NarrowAngle WideAngle
## Min. :5.050 Min. :5.000 Min. :4.950
## 1st Qu.:5.412 1st Qu.:5.463 1st Qu.:5.350
## Median :5.500 Median :5.525 Median :5.450
## Mean :5.543 Mean :5.534 Mean :5.459
## 3rd Qu.:5.650 3rd Qu.:5.600 3rd Qu.:5.550
## Max. :6.300 Max. :6.300 Max. :6.250boxplot(rt,
xlab = "Measurement Method", ylab = "Time (sec)",
main = "Rounding Time by Method",
col = c("lightblue", "lightgreen", "lightyellow"))
Question A: Is there a difference in timing between “Round Out” and “Narrow Angle”?
# H0: median(RoundOut) = median(NarrowAngle) Ha: ≠
wilcox.test(x = rt$RoundOut,
y = rt$NarrowAngle,
paired = TRUE,
alternative = "two.sided")##
## Wilcoxon signed rank exact test
##
## data: rt$RoundOut and rt$NarrowAngle
## V = 154.5, p-value = 0.288
## alternative hypothesis: true location shift is not equal to 0Question B: Is the timing under “Narrow Angle” greater than under “Wide Angle”?
# H0: median(NarrowAngle) = median(WideAngle) Ha: NarrowAngle > WideAngle
wilcox.test(x = rt$NarrowAngle,
y = rt$WideAngle,
paired = TRUE,
alternative = "greater")##
## Wilcoxon signed rank exact test
##
## data: rt$NarrowAngle and rt$WideAngle
## V = 222, p-value = 0.0005486
## alternative hypothesis: true location shift is greater than 0Why
paired = TRUE? Each operator produced a time under all three methods, so observations are matched by operator. Usingpaired = TRUEaccounts for individual differences between operators and increases the power of the test.
Assignment: Sleep and Precipitation Data (Case Study F — sleep, precip)
R’s built-in precip records average annual precipitation (inches) for 70 US cities (a one-sample design).
## Named num [1:70] 67 54.7 7 48.5 14 17.2 20.7 13 43.4 40.2 ...
## - attr(*, "names")= chr [1:70] "Mobile" "Juneau" "Phoenix" "Little Rock" ...One-sample test. Meteorologists often cite 35 inches as a typical annual precipitation figure for US cities. Is the median annual precipitation across these 70 cities significantly different from 35 inches? Use the Wilcoxon signed-rank test. State your hypotheses, test statistic, p-value, and conclusion (4-step procedure).
One-sample test, one-sided. Is the median annual precipitation less than 40 inches? Repeat the 4-step procedure with the appropriate one-sided alternative.
R’s built-in sleep dataset records the increase in hours of sleep for 10 patients, each tested under two soporific drugs (a paired design, matched by patient).
## 'data.frame': 20 obs. of 3 variables:
## $ extra: num 0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0 2 ...
## $ group: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ ID : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...Paired test. Each patient’s sleep was measured under both drugs, so the two extra measurements per patient are matched. Is there a significant difference in the median extra sleep gained between Drug 1 and Drug 2? Use wilcox.test() with paired = TRUE. State your hypotheses, test statistic, p-value, and conclusion.
Paired test, one-sided. Is Drug 2 associated with significantly more extra sleep than Drug 1? Repeat the 4-step procedure with the appropriate one-sided alternative.
Lab 7: Non-Parametric Statistics II
Objectives: Students are able to use R to:
- Apply the Mann-Whitney U test for two independent groups
- Apply the Kruskal-Wallis test for three or more independent groups
- Apply the Friedman test for three or more related groups
- Interpret the results in context
Choosing the right test:
| Test | Design | Parametric equivalent |
|---|---|---|
| Mann-Whitney U | 2 independent groups | Independent samples t-test |
| Kruskal-Wallis | 3+ independent groups | One-way ANOVA |
| Friedman | 3+ related groups (blocks) | Repeated-measures ANOVA |
Note: A significant Kruskal-Wallis or Friedman result tells us that at least one group differs from the others — it does not tell us which specific pairs differ.
Practice I: ToothGrowth Data
📊 Built-in R dataset:
ToothGrowth— no download requiredVariables:
len(tooth length),supp(supplement: OJ or VC),dose(vitamin C dose in mg/day: 0.5, 1.0, 2.0)
Import and prepare data
data("ToothGrowth")
head(ToothGrowth)
str(ToothGrowth)
# Convert dose to a factor for group comparisons
ToothGrowth$dose <- as.factor(ToothGrowth$dose)Question A: Is there a difference in tooth length between the two supplement types?
Two independent groups → Mann-Whitney U test
tapply(ToothGrowth$len, ToothGrowth$supp, summary)
tapply(ToothGrowth$len, ToothGrowth$supp, sd)
boxplot(len ~ supp, data = ToothGrowth,
xlab = "Supplement Type", ylab = "Tooth Length",
main = "Tooth Length by Supplement Type",
col = c("lightblue", "lightcoral"),
names = c("Orange Juice", "Ascorbic Acid"))
# H0: median length is equal for OJ and VC Ha: ≠
wilcox.test(len ~ supp, data = ToothGrowth,
paired = FALSE,
alternative = "two.sided")Interpretation: The Mann-Whitney U test compares the medians of two independent groups. - p-value < 0.05: reject \(H_0\); supplement type significantly affects tooth length. - p-value ≥ 0.05: insufficient evidence that supplement type affects tooth length.
Question B: Does tooth length differ across the three dose levels?
Three independent groups → Kruskal-Wallis test
tapply(ToothGrowth$len, ToothGrowth$dose, summary)
boxplot(len ~ dose, data = ToothGrowth,
xlab = "Dose (mg/day)", ylab = "Tooth Length",
main = "Tooth Length by Dose Level",
col = c("lightblue", "lightgreen", "lightyellow"))
# H0: median length equal across all dose levels Ha: at least one differs
kruskal.test(len ~ dose, data = ToothGrowth)Interpretation: The Kruskal-Wallis test compares the medians of three or more independent groups. - p-value < 0.05: reject \(H_0\); at least one dose level has a different median tooth length. - p-value ≥ 0.05: insufficient evidence that dose level affects tooth length.
Practice II: PlantGrowth and iris Data
📊 Built-in R datasets:
PlantGrowthandiris— no download requiredPlantGrowth variables:
weight(dried plant weight),group(ctrl, trt1, trt2)iris variables:
Sepal.Length,Sepal.Width,Petal.Length,Petal.Width,Species(setosa, versicolor, virginica)
Step 1 — Import and prepare data
Question A: Does the treatment affect plant weight?
tapply(PlantGrowth$weight, PlantGrowth$group, summary)
boxplot(weight ~ group, data = PlantGrowth,
xlab = "Treatment Group", ylab = "Dried Weight (g)",
main = "Plant Weight by Treatment",
col = c("lightgray", "lightblue", "lightgreen"),
names = c("Control", "Treatment 1", "Treatment 2"))
# H0: median weight equal across all groups Ha: at least one differs
kruskal.test(weight ~ group, data = PlantGrowth)Question B: Does sepal length differ across the three iris species?
tapply(iris$Sepal.Length, iris$Species, summary)
boxplot(Sepal.Length ~ Species, data = iris,
xlab = "Species", ylab = "Sepal Length (cm)",
main = "Sepal Length by Iris Species",
col = c("lightblue", "lightgreen", "lightyellow"))
# H0: median sepal length equal across all species Ha: at least one differs
kruskal.test(Sepal.Length ~ Species, data = iris)Practice III: RoundingTimes Data
📊 Built-in R dataset:
RoundingTimesfrom theMASSpackage — no download requiredDesign: Each of 15 subjects performed the rounding task under all five methods. This is a repeated measures design with 3+ related groups → use the Friedman test.
Step 1 — Import and reshape data
library(MASS)
data("RoundingTimes")
head(RoundingTimes)
# Reshape from wide to long format for friedman.test()
rt <- stack(RoundingTimes)
names(rt) <- c("Time", "Method")
rt$Subject <- factor(rep(1:nrow(RoundingTimes), times = ncol(RoundingTimes)))
head(rt)Step 2 — Explore data
tapply(rt$Time, rt$Method, summary)
boxplot(Time ~ Method, data = rt,
xlab = "Rounding Method", ylab = "Response Time (s)",
main = "Response Time by Rounding Method",
col = c("lightblue","lightgreen","lightyellow","lightcoral","lavender"))Question A: Do the five rounding methods produce different response times?
Three or more related groups → Friedman test
# H0: response time distributions are equal across all methods
# Ha: at least one method differs
friedman.test(Time ~ Method | Subject, data = rt)Interpretation: The Friedman test compares three or more related groups. - p-value < 0.05: reject \(H_0\); at least one rounding method produces a different distribution of response times. - p-value ≥ 0.05: insufficient evidence that rounding method affects response time.
Why Friedman and not Kruskal-Wallis? Each subject completed all five methods, so observations are matched by subject. The Friedman test accounts for this block structure; Kruskal-Wallis would treat all observations as independent and lose power.
Assignment: Lab 7
Use built-in R datasets. For each test, follow the 4-step procedure: (1) state hypotheses, (2) state the test and test statistic, (3) give the p-value, (4) state your conclusion in context.
- Using
ToothGrowth, restrict the data to a dose of 1.0 mg/day only, then test whether tooth length differs between the two supplement types at \(\alpha = 0.05\).
Repeat Question 1 for a dose of 2.0 mg/day. Compare with Question 1 — does dose appear to modify the supplement effect?
Using
iris, test whether petal length (Petal.Length) differs across the three species at \(\alpha = 0.05\).Using
iris, test whether sepal width (Sepal.Width) differs across the three species at \(\alpha = 0.05\).Using
RoundingTimes, test whether rounding method affects response time at \(\alpha = 0.05\). State clearly why the Friedman test is appropriate for this design.
Lab 8: Non-Parametric Statistics III
Objectives: Students are able to use R to:
- Compute and interpret the Spearman rank correlation coefficient
- Test for a significant monotone association between two variables
- Construct a contingency table and check expected frequencies
- Apply the chi-squared test of independence
- Interpret the results in context
Key concepts:
The Spearman rank correlation \(r_s\) measures the strength and direction of a monotone (not necessarily linear) association:
\[r_s = 1 - \frac{6\sum_{i=1}^{n} D_i^2}{n(n^2-1)}, \qquad D_i = R(X_i) - R(Y_i)\]
where \(R(X_i)\) and \(R(Y_i)\) are the ranks of observation \(i\). The range is \(-1 \leq r_s \leq 1\), interpreted with the same strength table as Pearson’s \(r\).
The chi-squared test of independence tests whether two categorical variables are associated:
\[\chi^2 = \sum_{i,j} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}, \qquad E_{ij} = \frac{(\text{Row } i \text{ total}) \times (\text{Col } j \text{ total})}{N}\]
Validity condition: All expected cell counts must be \(\geq 5\). If any are smaller, use Fisher’s exact test (
fisher.test()) instead.
Practice I: trees Data
📊 Built-in R dataset:
trees— no download requiredVariables:
Girth(diameter in inches at 4.5 ft),Height(feet),Volume(cubic feet of timber) for 31 black cherry trees
Import and prepare data
data("trees")
head(trees)
str(trees)
summary(trees)
# Scatterplot matrix
pairs(trees, main = "Scatterplot Matrix: trees Data",
pch = 19, col = "steelblue")Question A: Is there a positive monotone association between girth and volume?
Two quantitative variables, non-normal → Spearman rank correlation
plot(trees$Girth, trees$Volume,
xlab = "Girth (inches)", ylab = "Volume (cubic ft)",
main = "Tree Girth vs. Timber Volume",
pch = 19, col = "steelblue")
# H0: ρ_s ≤ 0 Ha: ρ_s > 0
cor.test(trees$Girth, trees$Volume,
method = "spearman",
alternative = "greater")Interpretation of
cor.test()output:
Output field What it means rhoSample Spearman correlation \(r_s\) STest statistic based on rank differences p-valueReject \(H_0\) if p-value < 0.05
- p-value < 0.05: reject \(H_0\); there is a significant positive monotone association.
- p-value ≥ 0.05: insufficient evidence of a positive monotone association.
Question B: Is there a monotone association between height and volume?
plot(trees$Height, trees$Volume,
xlab = "Height (feet)", ylab = "Volume (cubic ft)",
main = "Tree Height vs. Timber Volume",
pch = 19, col = "steelblue")
# H0: ρ_s = 0 Ha: ρ_s ≠ 0
cor.test(trees$Height, trees$Volume,
method = "spearman",
alternative = "two.sided")Question C: How does Spearman compare with Pearson for girth and volume?
# Pearson correlation (assumes linearity and normality)
cor.test(trees$Girth, trees$Volume, method = "pearson")
# Spearman correlation (monotone, robust to outliers)
cor.test(trees$Girth, trees$Volume, method = "spearman")Spearman vs. Pearson: Spearman measures monotone association (consistently increasing or decreasing), not just linear. It is more robust to outliers and does not require normality. If the relationship is monotone but curved, Spearman will give a higher value than Pearson.
Practice II: HairEyeColor Data
📊 Built-in R dataset:
HairEyeColor— no download requiredVariables:
Hair(Black, Brown, Red, Blond),Eye(Brown, Blue, Hazel, Green),Sex(Male, Female) for 592 statistics students
Step 1 — Import and prepare data
data("HairEyeColor")
str(HairEyeColor)
# Collapse across Sex to get a two-way Hair × Eye table
tab <- margin.table(HairEyeColor, margin = c(1, 2))
tab
prop.table(tab, margin = 1) # row percentagesStep 2 — Explore data
barplot(prop.table(tab, margin = 1),
beside = TRUE,
legend = rownames(tab),
xlab = "Eye Colour", ylab = "Proportion",
main = "Hair Colour Distribution by Eye Colour",
col = c("tan4", "goldenrod", "firebrick", "black"),
args.legend = list(title = "Hair Colour", bty = "n"))Question A: Is there a relationship between hair colour and eye colour?
Two categorical variables → Chi-squared test of independence
# Always check expected counts BEFORE running the test
chisq.test(tab)$expected
# H0: hair colour and eye colour are independent Ha: they are associated
chisq.test(tab)Check the chi-squared conditions: All expected cell counts should be ≥ 5. If any expected count < 5, use Fisher’s exact test:
fisher.test(tab).
Question B: Which combinations drive the association?
# Standardised residuals: |residual| > 2 indicates a cell contributing strongly
chisq.test(tab)$stdresInterpretation: A large positive residual means more observations than expected under independence (those categories occur together more often than chance). A large negative residual means fewer than expected.
Practice III: survey Data
📊 Built-in R dataset:
surveyfrom theMASSpackage — no download requiredVariables:
Sex(Female, Male),Smoke(Never, Occas, Regul, Heavy),Exer(Freq, None, Some),W.Hnd(writing hand) for 237 statistics students
Step 1 — Import and prepare data
library(MASS)
data("survey")
str(survey)
# Remove missing values for the variables of interest
sv <- survey[!is.na(survey$Smoke) & !is.na(survey$Exer) & !is.na(survey$Sex), ]Question A: Is smoking frequency related to exercise frequency?
tab1 <- table(sv$Smoke, sv$Exer)
tab1
prop.table(tab1, margin = 1)
# Check expected counts
chisq.test(tab1)$expected
# H0: smoking and exercise frequency are independent Ha: associated
chisq.test(tab1)
# If any expected count < 5, use Fisher's exact test:
fisher.test(tab1)Question B: Is smoking frequency related to sex?
tab2 <- table(sv$Sex, sv$Smoke)
tab2
prop.table(tab2, margin = 1)
chisq.test(tab2)$expected
# H0: sex and smoking frequency are independent Ha: associated
chisq.test(tab2)Assignment: Lab 8
Use built-in R datasets. For each test, follow the 4-step procedure: (1) state hypotheses, (2) state the test and test statistic, (3) give the p-value, (4) state your conclusion in context.
Using
trees, test whether there is a significant positive monotone association between girth and height at \(\alpha = 0.05\). Interpret \(r_s\) (direction and strength).Using
mtcars, test whether there is a significant negative monotone association between weight (wt) and fuel efficiency (mpg) at \(\alpha = 0.05\).Using
mtcars, test whether there is a significant monotone association (two-sided) between horsepower (hp) and fuel efficiency (mpg) at \(\alpha = 0.05\). Also compute Pearson’s \(r\) for the same pair and compare the two values.Using
surveyfromMASS, test whether exercise frequency (Exer) is associated with sex (Sex) at \(\alpha = 0.05\). Show the contingency table with row percentages and check all expected cell counts before running the test.Using
survey, test whether writing hand (W.Hnd) is associated with sex (Sex) at \(\alpha = 0.05\). Check the validity conditions and choose the appropriate test.
End of Lab Materials