Assignment 10
IE 5342 | Variance Stabilizing Transformation and Kruskal-Wallis
Jeffrey Mateo Hanna Hernandez
2026-10-02
Introduction
A civil engineer wants to know if four methods of estimating flood flow frequency give the same peak discharge on the same watershed. Each method was used six times, so there are 24 observations in total.
Data
The data is entered in tidy format: one column with the method and one column with the discharge.
discharge <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12,
0.91, 2.94, 2.14, 2.36, 2.86, 4.55,
6.31, 8.37, 9.75, 6.09, 9.82, 7.24,
17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
method <- factor(rep(1:4, each = 6))
dat <- data.frame(method, discharge)
head(dat)## method discharge
## 1 1 0.34
## 2 1 0.12
## 3 1 1.23
## 4 1 0.70
## 5 1 1.75
## 6 1 0.12str(dat)## 'data.frame': 24 obs. of 2 variables:
## $ method : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 2 2 2 2 ...
## $ discharge: num 0.34 0.12 1.23 0.7 1.75 0.12 0.91 2.94 2.14 2.36 ...tapply(dat$discharge, dat$method, mean)## 1 2 3 4
## 0.710000 2.626667 7.930000 14.718333tapply(dat$discharge, dat$method, sd)## 1 2 3 4
## 0.661090 1.192202 1.647070 2.795521a. Linear effects model and hypothesis
The linear effects model for a completely randomized design is:
y_ij = mu + tau_i + e_ij
where
y_ijis observation j of method imuis the overall meantau_iis the effect of method i (how much that method moves away from the overall mean)e_ijis the error, assumed normal with mean 0 and constant variance- i goes from 1 to 4 (the four methods) and j goes from 1 to 6 (the six observations of each method)
The hypotheses are:
- H0: tau1 = tau2 = tau3 = tau4 = 0 (no method has an effect, so all four have the same mean discharge)
- Ha: at least one tau_i is not 0 (at least one method gives a different mean discharge)
b. Are the data normal? Is the variance constant?
boxplot(discharge ~ method, data = dat,
col = "lightblue",
main = "Peak Discharge by Estimation Method",
xlab = "Estimation Method",
ylab = "Peak Discharge (cubic feet per second)")
tapply(dat$discharge, dat$method, var)## 1 2 3 4
## 0.437040 1.421347 2.712840 7.814937c. Fit the model with aov() and look at the residuals
model <- aov(discharge ~ method, data = dat)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## method 3 708.7 236.2 76.29 4e-11 ***
## Residuals 20 61.9 3.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qqnorm(residuals(model), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(model))
plot(fitted(model), residuals(model),
main = "Residuals vs Fitted Values",
xlab = "Fitted Value (mean discharge of the method)",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
Comments
The normal Q-Q plot of the residuals is close to the line, so normality is not the main problem here.
The residuals vs fitted plot is the one that shows the problem. On the left the points are very close to the zero line, and as we move to the right they get farther and farther away. The residuals open up like a funnel, which means the variance grows with the mean. This breaks the constant variance condition that the model needs.
So the model on the raw data is not adequate and the data needs a transformation.
d. Box Cox transformation and test
The Box Cox method looks for the power lambda that works
best for the data. The boxcox() function is in the MASS
package.
library(MASS)
boxcox(model, lambda = seq(-2, 2, 0.1))
bc <- boxcox(model, lambda = seq(-2, 2, 0.01), plotit = FALSE)
bc$x[which.max(bc$y)]## [1] 0.54The best lambda is 0.54, which is very close to 0.5. A lambda of 0.5 is the square root, so we use the square root transformation. We pick 0.5 instead of 0.54 because it is easier to explain and 0.5 is inside the confidence interval of the plot.
dat$sqrt_discharge <- sqrt(dat$discharge)
model2 <- aov(sqrt_discharge ~ method, data = dat)
summary(model2)## Df Sum Sq Mean Sq F value Pr(>F)
## method 3 32.69 10.898 81.17 2.27e-11 ***
## Residuals 20 2.69 0.134
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Checking the residuals again:
qqnorm(residuals(model2), main = "Normal Q-Q Plot of the Residuals (square root)")
qqline(residuals(model2))
plot(fitted(model2), residuals(model2),
main = "Residuals vs Fitted Values (square root)",
xlab = "Fitted Value",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
tapply(dat$sqrt_discharge, dat$method, sd)## 1 2 3 4
## 0.4044534 0.3857295 0.2929908 0.3726343Comments
The funnel is gone. In the new residuals vs fitted plot the points are spread about the same amount from left to right. The standard deviations are now between 0.29 and 0.40, instead of between 0.66 and 2.80 on the raw data, so the variance is much more constant. The Q-Q plot still looks fine.
Now the test. Using alpha = 0.05, the F value is 81.17 and the p-value is 2.27e-11, which is almost zero. Since the p-value is much smaller than 0.05, we reject H0.
There is strong evidence that the four methods do not give the same mean peak discharge.
e. Kruskal-Wallis test
The Kruskal-Wallis test is the nonparametric version of the one-way ANOVA. It works with the ranks of the data instead of the values, so it does not need the data to be normal or the variance to be constant. We run it on the raw data.
kruskal.test(discharge ~ method, data = dat)##
## Kruskal-Wallis rank sum test
##
## data: discharge by method
## Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05Comments
The test statistic is 21.156 with 3 degrees of freedom and the p-value is 9.771e-05. Using alpha = 0.05, the p-value is much smaller, so we reject H0 again.
Both tests reach the same conclusion: the four estimation methods do not produce equivalent estimates of peak discharge. The Kruskal-Wallis test is a good check here because it does not need the transformation, and it agrees with the result we got from the transformed ANOVA.
Complete code
# ---- data ----
discharge <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12,
0.91, 2.94, 2.14, 2.36, 2.86, 4.55,
6.31, 8.37, 9.75, 6.09, 9.82, 7.24,
17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
method <- factor(rep(1:4, each = 6))
dat <- data.frame(method, discharge)
head(dat)
str(dat)
tapply(dat$discharge, dat$method, mean)
tapply(dat$discharge, dat$method, sd)
# ---- b. normality and constant variance ----
boxplot(discharge ~ method, data = dat,
col = "lightblue",
main = "Peak Discharge by Estimation Method",
xlab = "Estimation Method",
ylab = "Peak Discharge (cubic feet per second)")
tapply(dat$discharge, dat$method, var)
# ---- c. model on the raw data ----
model <- aov(discharge ~ method, data = dat)
summary(model)
qqnorm(residuals(model), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(model))
plot(fitted(model), residuals(model),
main = "Residuals vs Fitted Values",
xlab = "Fitted Value (mean discharge of the method)",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
# ---- d. Box Cox and transformed model ----
library(MASS)
boxcox(model, lambda = seq(-2, 2, 0.1))
bc <- boxcox(model, lambda = seq(-2, 2, 0.01), plotit = FALSE)
bc$x[which.max(bc$y)]
dat$sqrt_discharge <- sqrt(dat$discharge)
model2 <- aov(sqrt_discharge ~ method, data = dat)
summary(model2)
qqnorm(residuals(model2), main = "Normal Q-Q Plot of the Residuals (square root)")
qqline(residuals(model2))
plot(fitted(model2), residuals(model2),
main = "Residuals vs Fitted Values (square root)",
xlab = "Fitted Value",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
tapply(dat$sqrt_discharge, dat$method, sd)
# ---- e. Kruskal-Wallis ----
kruskal.test(discharge ~ method, data = dat)
Comments
The variance is clearly not constant. Method 1 has a variance of 0.44 and method 4 has a variance of 7.81, almost 18 times bigger. In the boxplot the boxes get taller from left to right.
The problem is that the spread grows together with the mean. Method 1 has the smallest mean (0.71) and the smallest box, method 4 has the biggest mean (14.72) and the biggest box. This is the typical case where a transformation helps.
For normality, each group only has 6 observations, which is too few to judge alone, so the normality is checked with the residuals of the model in part c.