Assignment10
Priscilla Garcia
Last compiled on October 05, 2026 at 4:42 PM - CDT
1 Exercise 1
A civil engineer is interested in determining whether four different methods of estimating flood flow frequency produce equivalent estimates of peak discharge when applied to the same watershed. Each procedure is used six times on the watershed, and the resulting discharge data (in cubic feet per second) are shown below.
- Write the linear effects equation and the hypothesis you are testing.
Linear Effects Equation
\[Y_{ij} = \mu + \tau_i + \epsilon_{ij}\] Where:\(Y_{ij}\) = The \(j\)-th peak discharge observation (\(j = 1, 2, \dots, 6\)) under the \(i\)-th flood estimation method (\(i = 1, 2, 3, 4\)). \(\mu\) = The overall mean peak discharge. \(\tau_i\) = The effect of the \(i\)-th estimation method. \(\epsilon_{ij}\) = Random experimental error component, assumed to be \(\text{N}(0, \sigma^2)\).
Hypothesis:
Null Hypothesis (\(H_0\)): \(\tau_1 = \tau_2 = \tau_3 = \tau_4 = 0\) (All methods produce equivalent mean peak discharge estimates).
Alternative Hypothesis (\(H_1\)): At least one \(\tau_i \neq 0\) (At least one estimation method has a different mean peak discharge).
library(ggplot2)
library(MASS)
library(car)
# Data
method1 <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12)
method2 <- c(0.91, 2.94, 2.14, 2.36, 2.86, 4.55)
method3 <- c(6.31, 8.37, 9.75, 6.09, 9.82, 7.24)
method4 <- c(17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
discharge <- c(method1, method2, method3, method4)
method <- factor(rep(c("Method 1", "Method 2", "Method 3", "Method 4"), each = 6))
df <- data.frame(discharge, method)
print(df)## discharge method
## 1 0.34 Method 1
## 2 0.12 Method 1
## 3 1.23 Method 1
## 4 0.70 Method 1
## 5 1.75 Method 1
## 6 0.12 Method 1
## 7 0.91 Method 2
## 8 2.94 Method 2
## 9 2.14 Method 2
## 10 2.36 Method 2
## 11 2.86 Method 2
## 12 4.55 Method 2
## 13 6.31 Method 3
## 14 8.37 Method 3
## 15 9.75 Method 3
## 16 6.09 Method 3
## 17 9.82 Method 3
## 18 7.24 Method 3
## 19 17.15 Method 4
## 20 11.82 Method 4
## 21 10.97 Method 4
## 22 17.20 Method 4
## 23 14.35 Method 4
## 24 16.82 Method 4- Does it appear the data is normally distributed? Does it appear that the variance is constant?
# Fit Raw Model
model_raw <- aov(discharge ~ method, data = df)
# Residual Plots
par(mfrow = c(1, 2))
plot(model_raw, which = 1) # Residuals vs Fitted
plot(model_raw, which = 2) # Q-Q plot
par(mfrow = c(1, 1))
summary(model_raw)## 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 ' ' 1The plots indicate a normal distribution, aligned with the diagonal reference line. However, the Residuals vs Fitted plot demonstrated that the variance is not constant, especially for the Method 4, where the distribution is wider.
- (parametric) Fit the model using aov() in R on the raw data. Examine the residuals.
The Residuals vs. Fitted plot exhibits an increasing spread as fitted values increase.The ANOVA test on the raw data yields as F statistic of 76.29 and a p-value far below the significance level of alpha = 0.05. Therefore, we can reject the null hypothesis. Additionally, the evaluation of the residuals reveals a violation of the constant variance assumption. Because of that, it is necessary apply a variance-stabilizing transformation, for example a Box-Cox.
- (parametric) Select an appropriate transformation using Box Cox, transform the data and test hypothesis in R (alpha=0.05).
# Lambda
bc <- boxcox(discharge ~ method, data = df, lambda = seq(-1, 2, by = 0.1))
lambda <- bc$x[which.max(bc$y)]
cat("Optimal Lambda:", lambda, "\n")## Optimal Lambda: 0.5454545# Box-Cox transformation
if (lambda == 0) {
df$trans_discharge <- log(df$discharge)
} else {
df$trans_discharge <- (df$discharge^lambda - 1) / lambda
}
# Fit model on transformed data
model_trans <- aov(trans_discharge ~ method, data = df)
summary(model_trans)## Df Sum Sq Mean Sq F value Pr(>F)
## method 3 149.66 49.89 83.69 1.71e-11 ***
## Residuals 20 11.92 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Applying the transformation with Box-Cox, we can identify the optimal
power transformation (lambda) that stabilizes the variance across
groups. The lambda found was 0.55
After the transformation, we can reject H0 and conclude the methods
differ significantly, p-value < 0.05.
- (nonparametric) Perform a Kruskal-Wallace test in R (alpha=0.05).
kruskal.test(discharge ~ method, data = df)##
## Kruskal-Wallis rank sum test
##
## data: discharge by method
## Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05Because p < 0.05, we can reject H0. The non-parametric test indicates that the four methods yield significantly different peak discharge estimates.
2 Complete R Code
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message = FALSE)
library(ggplot2)
library(MASS)
library(car)
# Data
method1 <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12)
method2 <- c(0.91, 2.94, 2.14, 2.36, 2.86, 4.55)
method3 <- c(6.31, 8.37, 9.75, 6.09, 9.82, 7.24)
method4 <- c(17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
discharge <- c(method1, method2, method3, method4)
method <- factor(rep(c("Method 1", "Method 2", "Method 3", "Method 4"), each = 6))
df <- data.frame(discharge, method)
print(df)
# Fit Raw Model
model_raw <- aov(discharge ~ method, data = df)
# Residual Plots
par(mfrow = c(1, 2))
plot(model_raw, which = 1) # Residuals vs Fitted
plot(model_raw, which = 2) # Q-Q plot
par(mfrow = c(1, 1))
summary(model_raw)
# Lambda
bc <- boxcox(discharge ~ method, data = df, lambda = seq(-1, 2, by = 0.1))
lambda <- bc$x[which.max(bc$y)]
cat("Optimal Lambda:", lambda, "\n")
# Box-Cox transformation
if (lambda == 0) {
df$trans_discharge <- log(df$discharge)
} else {
df$trans_discharge <- (df$discharge^lambda - 1) / lambda
}
# Fit model on transformed data
model_trans <- aov(trans_discharge ~ method, data = df)
summary(model_trans)
kruskal.test(discharge ~ method, data = df)