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

Data entry

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)

flood_wide <- data.frame(method1, method2, method3, method4)

flood <- pivot_longer(
  flood_wide,
  cols=everything(),
  names_to="method",
  values_to="discharge"
)

flood$method <- factor(flood$method)
flood

## # A tibble: 24 × 2
##    method  discharge
##    <fct>       <dbl>
##  1 method1      0.34
##  2 method2      0.91
##  3 method3      6.31
##  4 method4     17.2 
##  5 method1      0.12
##  6 method2      2.94
##  7 method3      8.37
##  8 method4     11.8 
##  9 method1      1.23
## 10 method2      2.14
## # ℹ 14 more rows

descriptive_table <- flood %>%
  group_by(method) %>%
  summarise(
    n=n(),
    mean=mean(discharge),
    median=median(discharge),
    sd=sd(discharge),
    variance=var(discharge),
    .groups="drop"
  )

knitr::kable(descriptive_table, digits=3,
             caption="Descriptive statistics for the four methods")

Descriptive statistics for the four methods
method	n	mean	median	sd	variance
method1	6	0.710	0.520	0.661	0.437
method2	6	2.627	2.610	1.192	1.421
method3	6	7.930	7.805	1.647	2.713
method4	6	14.718	15.585	2.796	7.815

Question 1 Linear effects model and hypotheses

Let \(Y_{ij}\) be the \(j\)th peak-discharge estimate obtained using estimation method \(i\). The one-way fixed-effects model is

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

where

\(i=1,2,3,4\) identifies the estimation method;
\(j=1,2,\ldots,6\) identifies the observation;
\(\mu\) is the overall mean discharge;
\(\tau_i\) is the effect of method \(i\);
\(\epsilon_{ij}\) is the random error; and
\(\epsilon_{ij}\overset{iid}{\sim}N(0,\sigma^2)\).

The usual constraint is

\[\sum_{i=1}^{4}\tau_i=0.\]

The hypotheses are

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

or equivalently

\[H_0:\tau_1=\tau_2=\tau_3=\tau_4=0.\]

The alternative is

\[H_a:\text{at least one method mean differs}.\]

Question 2 Normality and constant variance

Graphical checks

par(mfrow=c(1,2))
boxplot(discharge~method, data=flood,
        xlab="Estimation method", ylab="Peak discharge",
        main="Peak discharge by method")
qqnorm(flood$discharge, main="Q-Q plot of pooled raw data")
qqline(flood$discharge, col="red", lwd=2)

par(mfrow=c(1,1))

The boxplot shows that the center and spread increase as the method number increases. The sample variances are approximately 0.437, 1.421, 2.713, and 7.815. Method 4 has much greater variation than Method 1. Therefore, constant variance does not appear reasonable for the raw response.

The pooled raw observations are right-skewed. However, ANOVA normality applies to the errors within the treatment groups, so the model residuals should also be examined.

shapiro.test(flood$discharge)

## 
##  Shapiro-Wilk normality test
## 
## data:  flood$discharge
## W = 0.8859, p-value = 0.01094

leveneTest(discharge~method, data=flood, center=median)

## 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

The pooled Shapiro-Wilk test gives approximately \(p=0.0109\), which indicates that the pooled raw observations are not normally distributed. Levene’s test gives approximately \(p=0.0135\), so the equal-variance assumption is rejected at \(\alpha=0.05\).

Conclusion: The main problem in the raw data is nonconstant variance. A variance-stabilizing transformation is appropriate.

Question 3 Parametric model using raw data

raw_model <- aov(discharge~method, data=flood)
summary(raw_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 raw-data ANOVA produces approximately

\[F_0=76.287,\qquad p=4.00\times10^{-11}.\]

Because \(p<0.05\), reject \(H_0\). The four estimation methods do not produce the same mean peak-discharge estimate.

Raw-model residuals

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

par(mfrow=c(1,1))

shapiro.test(residuals(raw_model))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(raw_model)
## W = 0.95693, p-value = 0.3798

leveneTest(discharge~method, data=flood, center=median)

## 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

The residual Q-Q plot is reasonably straight, and the Shapiro-Wilk test for the residuals gives approximately \(p=0.380\). Thus, residual normality is not a serious concern. The residual-versus-fitted plot shows a larger spread at the higher fitted values, and Levene’s test gives approximately \(p=0.0135\).

Conclusion: Although the raw ANOVA rejects equal means, the constant- variance assumption is not satisfied. The raw-model result should therefore be interpreted with caution, and a transformation should be used.

Question 4 Box Cox transformation and transformed ANOVA

Select the transformation

boxcox(raw_model, lambda=seq(-2,2,by=0.05),
       main="Box-Cox profile log-likelihood")

bc <- boxcox(raw_model, lambda=seq(-2,2,by=0.01), plotit=FALSE)
lambda_hat <- bc$x[which.max(bc$y)]
lambda_hat

## [1] 0.54

The Box-Cox estimate is approximately \(\widehat{\lambda}=0.54\). The nearby convenient value \(\lambda=0.50\) corresponds to the square-root transformation:

\[Y^*=\sqrt{Y}.\]

Because 0.50 is close to the estimated value, use the square-root transformation.

Fit the transformed model

flood$sqrt_discharge <- sqrt(flood$discharge)

sqrt_model <- aov(sqrt_discharge~method, data=flood)
summary(sqrt_model)

##             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

For the transformed response, the hypotheses are

\[H_0:\mu^*_1=\mu^*_2=\mu^*_3=\mu^*_4\]

against

\[H_a:\text{at least one transformed mean differs}.\]

The transformed ANOVA gives approximately

\[F_0=81.166,\qquad p=2.27\times10^{-11}.\]

Because \(p<0.05\), reject \(H_0\).

Conclusion: The four methods do not produce equivalent peak-discharge estimates.

Check the transformed-model assumptions

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

par(mfrow=c(1,1))

shapiro.test(residuals(sqrt_model))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sqrt_model)
## W = 0.95878, p-value = 0.4144

leveneTest(sqrt_discharge~method, data=flood, center=median)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.2387 0.8683
##       20

The transformed residuals are approximately normal: the Shapiro-Wilk p-value is about 0.414. The variance is also reasonably constant: the Levene p-value is about 0.868. Both values are greater than 0.05.

Conclusion: The square-root transformation corrects the unequal-variance problem, and the transformed ANOVA model is adequate.

Question 5 Kruskal Wallis test

The Kruskal-Wallis test compares the locations of the four method distributions without requiring normally distributed errors.

The hypotheses are

\[H_0:\text{the four methods have the same population distribution}\]

and

\[H_a:\text{at least one method has a different population distribution}.\]

kruskal.test(discharge~method, data=flood)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  discharge by method
## Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05

The test gives approximately

\[H=21.156,\qquad df=3,\qquad p=0.0000977.\]