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.\]
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.
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.
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.
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.
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)