1 Problem

A civil engineer is interested in determining whether four different methods of estimating flood flow frequency produce equivalent estimates of peak discharge when applied to the same watershed. Each procedure is used six times on the watershed, and the resulting discharge data (in cubic feet per second) are shown below.

2 Data Preparation

E1<-c(0.34,0.12,1.23,0.70,1.75,0.12)
E2<-c(0.91,2.94,2.14,2.36,2.86,4.55)
E3<-c(6.31,8.37,9.75,6.09,9.82,7.24)
E4<-c(17.15,11.82,10.97,17.20,14.35,16.82)
df<-data.frame(E1,E2,E3,E4)
library(tidyr)
df<-pivot_longer(df,c(E1,E2,E3,E4))
df$name<-as.factor(df$name)
df
## # A tibble: 24 × 2
##    name  value
##    <fct> <dbl>
##  1 E1     0.34
##  2 E2     0.91
##  3 E3     6.31
##  4 E4    17.2 
##  5 E1     0.12
##  6 E2     2.94
##  7 E3     8.37
##  8 E4    11.8 
##  9 E1     1.23
## 10 E2     2.14
## # ℹ 14 more rows

3 Solution

3.1 Part a

Write the linear effects equation and the hypothesis you are testing

\[x_{ij}=\mu+\tau_{i}+\varepsilon_{ij}\]

\[H_{0}: \tau_{i}=0\text{ for all } i\]

\[H_{a}: \tau_{i} \neq 0\text{ for some } i\]

3.2 Part b

Does it appear the data is normally distributed? Does it appear that the variance is constant?

aov.model <- aov(value ~ name, data = df)
qqnorm(residuals(aov.model))
qqline(residuals(aov.model))

boxplot(value ~ name, data = df)

The normal probability plot of the residuals falls reasonably close to the line, so there is no clear evidence against normality. The box plot shows the spread increasing with the group median; therefore, the variance does not appear to be constant.

3.3 Part c

(parametric) Fit the model using aov() in R on the raw data. Examine the residuals.

summary(aov.model)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## name         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
plot(aov.model, which = 1)

This test violates the one of the key assumptions of ANOVA:

  1. The variance is constant.

As shown in the residuals plot, the variance is not constant across all four estimation methods. Therefore, the results from this ANOVA test on the raw data cannot be trusted. Either some type of transformation needs to be performed on the data or use a non-parametric approach.

3.4 Part d

(parametric) Select an appropriate transformation using Box Cox, transform the data and test hypothesis in R (alpha=0.05)

library(MASS)
boxcox(value~name, data = df)

Based on the results from the boxcox analysis, the lambda used in this transformation will be:

\[\lambda=0.5\]

A lambda value of 0.5 lies within the 95% confidence interval for the boxcox test.

lambda<-0.5
df$value_l<-df$value^lambda
aov.model2<-aov(value_l~name,data = df)
summary(aov.model2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## name         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

Using an alpha value of 0.05 and the P value from the test = 2.27e-11, there is overwhelming evidence that the four estimation methods differ; therefore, the null hypothesis can be rejected.

plot(aov.model2, which = c(1:2))

To ensure the ANOVA model on the transformed data is adequate, the residuals will be checked for constant variance and normality. The variance appears to be fairly constant shown by the small deviations of the red line away from zero in the residuals vs fitted chart. There appears to be some minor s shaping to the normal probability plot, but overall the data appears to be very close to normally distributed.

3.5 Part e

(nonparametric) Perform a Kruskal-Wallis test in R (alpha=0.05)

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

Based on the results from the Kruskal-Wallis test, the p-value = 9.771e-05. This p-value is less than alpha; therefore, there is strong evidence that we can reject the null hypothesis.

4 Complete R Code

# Data Cleansing
E1<-c(0.34,0.12,1.23,0.70,1.75,0.12)
E2<-c(0.91,2.94,2.14,2.36,2.86,4.55)
E3<-c(6.31,8.37,9.75,6.09,9.82,7.24)
E4<-c(17.15,11.82,10.97,17.20,14.35,16.82)
df<-data.frame(E1,E2,E3,E4)
library(tidyr)
df<-pivot_longer(df,c(E1,E2,E3,E4))
df$name<-as.factor(df$name)
df

# Solution

## Part a
**Write the linear effects equation and the hypothesis you are testing**
$$x_{ij}=\mu+\tau_{i}+\varepsilon_{ij}$$
$$H_{0}: \tau_{i}=0\text{ for all } i$$
$$H_{a}: \tau_{i} \neq 0\text{ for some } i$$

## Part b
**Does it appear the data is normally distributed? Does it appear that the variance is constant?**
aov.model <- aov(value ~ name, data = df)
qqnorm(residuals(aov.model))
qqline(residuals(aov.model))
boxplot(value ~ name, data = df)

## Part c
**(parametric) Fit the model using aov() in R on the raw data.  Examine the residuals.**
summary(aov.model)
plot(aov.model, which = 1)

## Part d
**(parametric) Select an appropriate transformation using Box Cox, transform the data and test hypothesis in R (alpha=0.05)**
library(MASS)
boxcox(value~name, data = df)
lambda<-0.5
df$value_l<-df$value^lambda
aov.model2<-aov(value_l~name,data = df)
summary(aov.model2)
plot(aov.model2, which = c(1:2))

## Part e
kruskal.test(value~name, data = df)