Problem

Four methods of estimating flood flow frequency are each applied six times to the same watershed, and the resulting peak discharge (cfs) is recorded. The estimation method is the single factor with 4 levels and n = 6 observations per level (N = 24). The data are entered in a wide format and then converted to a tidy format using pivot_longer.

# A tibble: 8 × 2
  Method Discharge
  <fct>      <dbl>
1 1           0.34
2 2           0.91
3 3           6.31
4 4          17.2 
5 1           0.12
6 2           2.94
7 3           8.37
8 4          11.8

a)

The linear effects model is

\[y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \qquad i = 1, 2, 3, 4, \quad j = 1, 2, \ldots, 6\]

where \(y_{ij}\) is the \(j\)-th peak discharge estimate obtained with method \(i\), \(\mu\) is the overall mean, \(\tau_i\) is the effect of the \(i\)-th estimation method, and \(\varepsilon_{ij} \sim NID(0, \sigma^2)\) is the random error.

The hypotheses being tested are:

  • \(H_0: \tau_1 = \tau_2 = \tau_3 = \tau_4 = 0\) (all four methods produce equivalent mean estimates of peak discharge)
  • \(H_1: \tau_i \neq 0\) for at least one \(i\)

b)

The raw data are examined with a side-by-side boxplot and the mean, variance, and standard deviation of each method. A Shapiro-Wilk test is also run within each method as a rough check of normality.

  Method    Mean Variance     SD Shapiro_p
1      1  0.7100   0.4370 0.6611    0.2976
2      2  2.6267   1.4213 1.1922    0.8033
3      3  7.9300   2.7128 1.6471    0.2901
4      4 14.7183   7.8149 2.7955    0.1238

Normality: Within each method the six observations are spread fairly symmetrically and the boxplots show no outliers. The Shapiro-Wilk p-values for the four methods (0.298, 0.803, 0.290, 0.124) are all greater than 0.05, so it appears the data are approximately normally distributed. With only six observations per method this check has limited power, so normality is examined more carefully with the residuals in part c).

Constant variance: The variance does not appear to be constant. The boxes get wider as the mean discharge increases, and the sample variance grows from 0.437 for method 1 to 7.815 for method 4, which is about 18 times larger. Since the spread increases with the mean, a variance-stabilizing transformation is likely needed.

c)

            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 ' ' 1

On the raw data, the one-way ANOVA gives F₀ = 76.29 with a p-value of about 4.0 × 10⁻¹¹. Before relying on this result, the residuals are examined.


    Shapiro-Wilk normality test

data:  residuals(model)
W = 0.95693, p-value = 0.3798
Levene's Test for Homogeneity of Variance (center = median)
      Df F value  Pr(>F)  
group  3  4.5762 0.01348 *
      20                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In the normal probability plot, the residuals fall reasonably close to the reference line, and the Shapiro-Wilk test (p = 0.380) does not reject normality, so the normality assumption is acceptable.

The residuals vs. fitted plot, however, shows a clear funnel (megaphone) shape: the residuals are tightly grouped around zero at the smallest fitted value (0.71) and spread much wider at the largest fitted value (14.72). The modified Levene test gives p = 0.0135 < 0.05, so we reject the hypothesis of equal variances. The constant variance assumption is violated, so the ANOVA on the raw data is not fully reliable even though its F statistic is very large. This motivates the transformation in part d).

d)

The Box-Cox method is used to find the power transformation \(y^* = y^\lambda\) that maximizes the log-likelihood.

Estimated lambda: 0.54 
95% CI for lambda: [ 0.33 , 0.76 ]

The Box-Cox procedure gives λ̂ = 0.54 with a 95% confidence interval of [0.33, 0.76]. The interval contains λ = 0.5 but not λ = 1 (no transformation) or λ = 0 (log transformation). Since λ = 0.5 is inside the interval and is easier to interpret, the square root transformation \(y^* = \sqrt{y}\) is selected.

The hypotheses are the same as in part a), now stated for the transformed response, and are tested at α = 0.05.

            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 ' ' 1

On the transformed data, the one-way ANOVA gives F₀ = 81.17 with a p-value of about 2.3 × 10⁻¹¹. Since the p-value < α = 0.05, we reject H₀ and conclude that the four estimation methods do not produce equivalent estimates of peak discharge.

The residuals of the transformed model are checked to confirm the transformation worked.

  Method Variance_sqrt
1      1        0.1636
2      2        0.1488
3      3        0.0858
4      4        0.1389

    Shapiro-Wilk normality test

data:  residuals(model_t)
W = 0.95878, p-value = 0.4144
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3  0.2387 0.8683
      20

After the transformation, the funnel shape is gone: the residuals are scattered with a similar spread at every fitted value. The variances of \(\sqrt{y}\) now range only from 0.086 to 0.164, and the modified Levene test (p = 0.868) does not reject equal variances. Normality still holds (Shapiro-Wilk p = 0.414). The transformed model is adequate, so the conclusion of the test above is trustworthy.

e)

The Kruskal-Wallis test does not require the normality assumption. It replaces each observation by its rank among all N = 24 observations (tied values get the average rank) and tests

  • \(H_0\): the four methods have identical distributions of peak discharge
  • \(H_1\): at least one method tends to give larger (or smaller) values than the others
  Method Rank_Sum Mean_Rank
1      1       23    3.8333
2      2       55    9.1667
3      3       93   15.5000
4      4      129   21.5000

    Kruskal-Wallis rank sum test

data:  Discharge by Method
Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05
[1] 7.814728

The rank sums are 23, 55, 93, and 129 for methods 1 to 4, and the ranks barely overlap: every observation from method 3 ranks above all of method 1 and 2, and every observation from method 4 ranks above all of method 3. The test statistic is H = 21.16 (with the tie correction), with a p-value of about 9.8 × 10⁻⁵. Since H = 21.16 > χ²₀.₀₅,₃ = 7.81 (equivalently, p-value < α = 0.05), we reject H₀ and conclude that the four estimation methods do not produce equivalent peak discharge estimates.

This agrees with the parametric analysis on the square-root transformed data in part d). All approaches lead to the same conclusion: the choice of estimation method has a significant effect on the estimated peak discharge.

R Code

library(MASS)
library(dplyr)
library(tidyr)
library(ggplot2)
library(car)

# Data in tidy format
wide <- data.frame(
  Method_1 = c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12),
  Method_2 = c(0.91, 2.94, 2.14, 2.36, 2.86, 4.55),
  Method_3 = c(6.31, 8.37, 9.75, 6.09, 9.82, 7.24),
  Method_4 = c(17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
)

dat <- wide %>%
  pivot_longer(cols = everything(), names_to = "Method", names_prefix = "Method_", values_to = "Discharge")
dat$Method <- as.factor(dat$Method)

head(dat, 8)

# b) Normality and constant variance of the raw data
boxplot(Discharge ~ Method, data = dat,
        main = "Side by Side Boxplot: Peak Discharge by Method",
        xlab = "Estimation Method",
        ylab = "Peak Discharge (cfs)",
        col = c("blue", "red", "green", "orange"))

dat %>%
  group_by(Method) %>%
  summarise(Mean = round(mean(Discharge), 4),
            Variance = round(var(Discharge), 4),
            SD = round(sd(Discharge), 4),
            Shapiro_p = round(shapiro.test(Discharge)$p.value, 4)) %>%
  as.data.frame()

# c) One-way ANOVA on the raw data and residual analysis
model <- aov(Discharge ~ Method, data = dat)
summary(model)

res <- data.frame(Fitted = fitted(model), Residuals = residuals(model))

ggplot(res, aes(sample = Residuals)) + 
  stat_qq(color = "blue") +                              
  stat_qq_line() +
  labs(x = "Theoretical Quantiles", y = "Residuals", title = "Normal Probability Plot of Residuals (Raw Data)")

ggplot(res, aes(x = Fitted, y = Residuals)) +
  geom_point(color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Fitted Values", y = "Residuals", title = "Residuals vs Fitted Values (Raw Data)")

shapiro.test(residuals(model))
leveneTest(Discharge ~ Method, data = dat)

# d) Box-Cox transformation
bc <- boxcox(Discharge ~ Method, data = dat, lambda = seq(-1, 2, 0.01))

lambda_hat <- bc$x[which.max(bc$y)]
lambda_ci <- range(bc$x[bc$y > max(bc$y) - qchisq(0.95, 1) / 2])

cat("Estimated lambda:", lambda_hat, "\n")
cat("95% CI for lambda: [", lambda_ci[1], ",", lambda_ci[2], "]\n")

# Square root transformation (lambda = 0.5) and ANOVA (alpha = 0.05)
dat$SqrtDischarge <- sqrt(dat$Discharge)

model_t <- aov(SqrtDischarge ~ Method, data = dat)
summary(model_t)

res_t <- data.frame(Fitted = fitted(model_t), Residuals = residuals(model_t))

ggplot(res_t, aes(sample = Residuals)) + 
  stat_qq(color = "blue") +                              
  stat_qq_line() +
  labs(x = "Theoretical Quantiles", y = "Residuals", title = "Normal Probability Plot of Residuals (Square Root)")

ggplot(res_t, aes(x = Fitted, y = Residuals)) +
  geom_point(color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Fitted Values", y = "Residuals", title = "Residuals vs Fitted Values (Square Root)")

dat %>%
  group_by(Method) %>%
  summarise(Variance_sqrt = round(var(SqrtDischarge), 4)) %>%
  as.data.frame()

shapiro.test(residuals(model_t))
leveneTest(SqrtDischarge ~ Method, data = dat)

# e) Kruskal-Wallis test (alpha = 0.05)
dat %>%
  mutate(Rank = rank(Discharge)) %>%
  group_by(Method) %>%
  summarise(Rank_Sum = sum(Rank), Mean_Rank = round(mean(Rank), 4)) %>%
  as.data.frame()

kruskal.test(Discharge ~ Method, data = dat)

qchisq(0.95, df = 3)