Introduction

We compare four estimation methods at \(\alpha=0.05\), assuming that they are independent observations.

flood <- data.frame(
  method = factor(rep(1:4, each = 6)),
  dischrg = c(
    0.34, 0.12, 1.23, 0.70, 1.75, 0.12,
    0.91, 2.94, 2.14, 2.36, 2.86, 4.55,
    6.31, 8.37, 9.75, 6.09, 9.82, 7.24,
    17.15, 11.82, 10.97, 17.20, 14.35, 16.82
  )
)
flood
##    method dischrg
## 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

a)

The fixed-effects equation is:

\[Y_{ij} = \mu + \tau_i + \epsilon_{ij}\] where: \[\qquad i=1,2,3,4\quad j=1,\ldots,6.\]

\(\mu\) is the mean

\(\tau_i\) is the method effect

\(\epsilon_{ij}\) is random error

The hypothesis we are testing is:

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

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

b)

Here we check normality within methods and compare their variances.

aggregate(dischrg ~ method, data = flood, FUN = mean)
##   method   dischrg
## 1      1  0.710000
## 2      2  2.626667
## 3      3  7.930000
## 4      4 14.718333
aggregate(dischrg ~ method, data = flood, FUN = var)
##   method  dischrg
## 1      1 0.437040
## 2      2 1.421347
## 3      3 2.712840
## 4      4 7.814937
boxplot(dischrg ~ method, data = flood,
        main = "Discharge by Method", xlab = "Method",
        ylab = "Discharge (cubic feet per second)", col = "blue")

{par(mfrow = c(2, 2))} for (i in levels(flood$method)) { values <- flood$dischrg[flood$method == i] qqnorm(values, main = paste("Method", i), col = "blue") qqline(values, col = "red") } par(mfrow = c(1, 1)) tapply(flood$dischrg, flood$method, function(x) shapiro.test(x)$p.value)

All four of the p-values are nore than 0.05, so normality is approximately normal.

The variances increase from 0.4370 to 7.8149, hence variance is not constant.

c)

model <- aov(dischrg ~ method, data = flood)
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

here F value is 76.2868

\(p=4.004\times10^{-11}\).

Hence we reject \(H_0\),

par(mfrow = c(1, 2))
plot(fitted(model), residuals(model),xlab = "Fitted Values", ylab = "Residuals", main = "Residuals vs Fitted", col = "blue")

qqnorm(residuals(model), main = "Residual Q-Q Plot", col = "blue")

qqline(residuals(model), col = "red")

par(mfrow = c(1, 1))

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

We can see that the residual spread increases with fitted values.

It also appears to be normally distributed, however the variance is not constant.

d)

bxcx <- MASS::boxcox(model, lambda = seq(-2, 2, by = 0.01))

best_lambda <- bxcx$x[which.max(bxcx$y)]
best_lambda
## [1] 0.54

The Box-Cox test gives us \(\lambda\approx0.54\), so we use the nearby square-root transformation (\(\lambda=0.5\)):

\[Z_{ij}=\sqrt{Y_{ij}}.\]

flood$root_discharge <- sqrt(flood$dischrg)
model_root <- aov(root_discharge ~ method, data = flood)
summary(model_root)
##             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

\(H_0\): All four square-root-scale means are equal.

\(H_a\): At least one differs.

par(mfrow = c(1, 2))
plot(fitted(model_root), residuals(model_root),xlab = "Fitted Values", ylab = "Residuals",main = "Transformed Residuals", pch = 19, col = "blue")

qqnorm(residuals(model_root), main = "Transformed Q-Q Plot",
       pch = 19, col = "blue")

qqline(residuals(model_root), col = "red")

par(mfrow = c(1, 1))

shapiro.test(residuals(model_root))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_root)
## W = 0.95878, p-value = 0.4144
aggregate(root_discharge ~ method, data = flood, FUN = var)
##   method root_discharge
## 1      1     0.16358252
## 2      2     0.14878726
## 3      3     0.08584364
## 4      4     0.13885633
flood$root_deviation <- abs(flood$root_discharge -
  ave(flood$root_discharge, flood$method, FUN = median))
summary(aov(root_deviation ~ method, data = flood))
##             Df Sum Sq Mean Sq F value Pr(>F)
## method       3 0.0279 0.00931   0.239  0.868
## Residuals   20 0.7802 0.03901

The transformation improves constant variance (Levene \(p=0.8683\)). Residual normality is reasonable (Shapiro-Wilk \(p=0.4144\)).

F value is 81.1658

\(p=2.266\times10^{-11}\).

Hence we reject \(H_0\).

e)

\(H_0\): All four discharge distributions are equal. \(H_a\): At least one differs.

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

\(H=21.1559\), \(df=3\), \(p=0.00009771\).

Hence we reject \(H_0\)