a) Linear Effects Equation and Hypotheses

For this completely randomized design, the linear effects model is

\[Y_{ij} = \mu + \tau_i + \epsilon_{ij}\]

where:

The hypotheses: \[H_0:\mu_1=\mu_2=\mu_3=\mu_4\] versus \[H_a:\text{At least one estimation-method mean is different.}\]

Where,

\[\alpha=0.05.\]

b) Examination of Normality and Constant Variance

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)

flow <- c(method1, method2, method3, method4)

method <- factor(rep(1:4, each = 6))

data <- data.frame(method, flow)

data
##    method  flow
## 1       1  0.34
## 2       1  0.12
## 3       1  1.23
## 4       1  0.70
## 5       1  1.75
## 6       1  0.12
## 7       2  0.91
## 8       2  2.94
## 9       2  2.14
## 10      2  2.36
## 11      2  2.86
## 12      2  4.55
## 13      3  6.31
## 14      3  8.37
## 15      3  9.75
## 16      3  6.09
## 17      3  9.82
## 18      3  7.24
## 19      4 17.15
## 20      4 11.82
## 21      4 10.97
## 22      4 17.20
## 23      4 14.35
## 24      4 16.82

Before interpreting the ANOVA, the assumptions of normality and constant variance will be examined. First, a one-way ANOVA model is fitted to the original data.

qqnorm(data$flow,main = "Normal Q-Q Plot",xlab = "Theoretical Quantiles",ylab = "Sample Quantiles")

qqline(data$flow)

boxplot(flow ~ method,data = data,xlab = "Estimation Method",ylab = "Peak Discharge",main = "Boxplot of Peak Discharge by Estimation Method")

model <- aov(flow ~ method, data = data)

The diagnostic plots are:

par(mfrow = c(2, 2))
plot(model)

par(mfrow = c(1, 1))

Based on the diagnostic plots, the normality assumption is reasonable. However, the residual spread changes across the fitted values, indicating that the constant-variance assumption is questionable. So, a variance-stabilizing transformation should be considered.

c) Parametric ANOVA on the Raw Data

A one-way ANOVA is performed on the original peak-discharge data.

model <- aov(flow ~ method, data = data)
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 ' ' 1

The hypotheses are

\[H_0:\mu_1=\mu_2=\mu_3=\mu_4\]

and

\[H_a:\text{At least one mean is different.}\]

At the \(\alpha=0.05\) significance level, the p-value from the ANOVA is less than 0.05. So, the null hypothesis is rejected. Therefore, there is sufficient statistical evidence to conclude that the mean peak-discharge estimates are not the same for all four estimation methods.

Diagnostic plots for the ANOVA model:

par(mfrow = c(2, 2))
plot(model)

par(mfrow = c(1, 1))

The residual plots show that the spread of the residuals increases as the fitted values increase, indicating that the constant variance assumption is not satisfied. The Q-Q plot also shows some departures from normality, particularly in the tails. Therefore, the ANOVA assumptions are questionable, and a Box-Cox transformation should be considered.

d) Box-Cox Transformation and ANOVA

A Box-Cox analysis is used to identify an appropriate variance-stabilizing transformation.

library(MASS)

model.lm <- lm(flow ~ method, data = data)

boxcox(model.lm)

The Box-Cox plot indicates that \(\lambda\) is approximately \(0.5\), a square-root transformation is appropriate:

data$sqrt_flow <- sqrt(data$flow)

The one-way ANOVA is then repeated using the transformed response.

model.sqrt <- aov(sqrt_flow ~ method, data = data)
summary(model.sqrt)
##             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

At \(\alpha=0.05\), the p-value is less than 0.05, the null hypothesis is rejected. This indicates that significant differences among the estimation methods remain after accounting for the nonconstant variance through transformation.

par(mfrow = c(2, 2))
plot(model.sqrt)

par(mfrow = c(1, 1))

After applying the square-root transformation, the residuals show a more uniform spread across the fitted values, indicating that the constant variance assumption has improved. The Q-Q plot is also reasonably close to a straight line, with only minor deviations in the tails. So, the transformed data appear to satisfy the ANOVA assumptions more adequately than the raw data.

e) Kruskal-Wallis Test

The hypotheses are

\[H_0:\text{The distributions of peak-discharge estimates are the same for all four methods}\]

versus

\[H_a:\text{At least one estimation method differs}\]

The Kruskal-Wallis test is :

kruskal.test(flow ~ method, data = data)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  flow by method
## Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05

At

\[\alpha=0.05,\]

Since the p-value is less than 0.05, the null hypothesis is rejected.

So, there is sufficient statistical evidence to conclude that the four estimation methods do not all produce equivalent peak-discharge estimates.

Complete R Code

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)

flow <- c(method1, method2, method3, method4)
method <- factor(rep(1:4, each = 6))

data <- data.frame(method, flow)

data

qqnorm(data$flow, main = "NPP for Peak Discharge")
qqline(data$flow)

boxplot(flow ~ method,data = data,xlab = "Estimation Method",ylab = "Peak Discharge",main = "Peak Discharge by Estimation Method")


model <- aov(flow ~ method, data = data)

summary(model)

par(mfrow = c(2, 2))
plot(model)

par(mfrow = c(1, 1))

library(MASS)

model.lm <- lm(flow ~ method, data = data)

boxcox(model.lm)

data$sqrt_flow <- sqrt(data$flow)

model.sqrt <- aov(sqrt_flow ~ method, data = data)

summary(model.sqrt)

par(mfrow = c(2, 2))
plot(model.sqrt)

par(mfrow = c(1, 1))

kruskal.test(flow ~ method, data = data)