Introduction

This document contains the analysis of an experiment in which a civil engineer compares four different methods of estimating flood flow frequency. Each method is used six times on the same watershed, and the peak discharge (in cubic feet per second) is recorded. The goal is to determine whether the four methods produce equivalent estimates of peak discharge. First, the data is checked for normality and constant variance. Then a one-way ANOVA is fitted on the raw data, a variance stabilizing transformation is selected with Box-Cox and the hypothesis is tested on the transformed data, and finally the nonparametric Kruskal-Wallis test is performed. All of the data manipulation, analysis and plotting is done within R.

Data

m1<-c(0.34,0.12,1.23,0.70,1.75,0.12)
m2<-c(0.91,2.94,2.14,2.36,2.86,4.55)
m3<-c(6.31,8.37,9.75,6.09,9.82,7.24)
m4<-c(17.15,11.82,10.97,17.20,14.35,16.82)
obs<-c(m1,m2,m3,m4)
x<-as.factor(c(rep(1,6),rep(2,6),rep(3,6),rep(4,6)))
cat("sample mean of the peak discharge of each method \n",c(mean(m1),mean(m2),mean(m3),mean(m4)),"\n")
## sample mean of the peak discharge of each method 
##  0.71 2.626667 7.93 14.71833
cat("sample standard deviation of the peak discharge of each method \n",c(sd(m1),sd(m2),sd(m3),sd(m4)),"\n")
## sample standard deviation of the peak discharge of each method 
##  0.66109 1.192202 1.64707 2.795521

a. Linear Effects Equation and Hypothesis

The linear effects model is:

\[y_{ij} = \mu + \tau_i + \epsilon_{ij}, \qquad i = 1,2,3,4 \qquad j = 1,2,...,6\]

where \(y_{ij}\) is the \(j\)th observation of peak discharge with the \(i\)th estimation method, \(\mu\) is the overall mean, \(\tau_i\) is the effect of the \(i\)th method and \(\epsilon_{ij}\) is the random error, assumed \(NID(0,\sigma^2)\).

The hypothesis being tested is:

\[H_0: \tau_1 = \tau_2 = \tau_3 = \tau_4 = 0\] \[H_1: \tau_i \neq 0 \text{ for at least one } i\]

b. Normality and Constant Variance

boxplot(obs~x,xlab="Estimation Method",ylab="Peak Discharge (cfs)",main="Boxplot of Observations")

The boxplot shows that the spread of the observations is not the same for all of the methods. The boxes get larger as the mean of the method increases, going from a standard deviation of about 0.66 cfs in method 1 to about 2.80 cfs in method 4. This means that the variance increases with the mean, so the data does not appear to have constant variance. The boxes are reasonably symmetric and there are no extreme outliers, so the data does not show a strong departure from normality. These assumptions are checked again with the residuals of the model in part c.

c. ANOVA on the Raw Data

model<-aov(obs~x)
summary(model)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## x            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 F statistic is about 76.3 with a p-value much smaller than 0.05, so the raw data suggests that the methods do not produce the same mean peak discharge. However, this conclusion depends on the assumptions of the model, so the residuals are examined below.

meanx<-c(rep(mean(m1),6),rep(mean(m2),6),rep(mean(m3),6),rep(mean(m4),6))
res<-obs-meanx
qqnorm(res)

plot(meanx,res,xlab="method average",ylab="residual",main="constant variance wrt method means")

The points in the normal probability plot fall approximately along a straight line, so the normality assumption is reasonable. The plot of the residuals against the method means, however, shows a funnel shape: the residuals of methods with small means are close to zero, while the residuals of methods with large means are spread much wider. The constant variance assumption is violated, and since ANOVA is not robust to this violation, a variance stabilizing transformation is needed.

d. Box-Cox Transformation

library(MASS)
boxcox(obs~x)

The maximum of the log-likelihood is close to \(\lambda = 0.5\), and the 95% confidence interval on \(\lambda\) goes from about 0.36 to 0.76. Since the value of 1 is not in the confidence interval, a transformation is needed. Because \(\lambda\) is not zero, the data is raised to the power \(\lambda = 0.5\), which is the square root transformation.

lambda=0.5
obs2<-obs^(lambda)
boxcox(obs2~x)

After the transformation, the confidence interval on \(\lambda\) contains the value of 1, so no further transformation is needed.

boxplot(obs2~x,xlab="Estimation Method",ylab="Square Root of Peak Discharge",main="Boxplot of Transformed Observations")

model2<-aov(obs2~x)
summary(model2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            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
meanx2<-c(rep(mean(m1^lambda),6),rep(mean(m2^lambda),6),rep(mean(m3^lambda),6),rep(mean(m4^lambda),6))
res2<-obs2-meanx2
qqnorm(res2)

plot(meanx2,res2,xlab="method average",ylab="residual",main="constant variance wrt method means")

With the transformed data the boxes have a similar size, and the residuals have about the same spread for all of the methods, so the constant variance assumption is now satisfied. The normal probability plot is still approximately a straight line. The F statistic is about 81.2 with a p-value much smaller than alpha = 0.05, so we reject the null hypothesis. At least one of the four methods produces a different mean estimate of peak discharge, so the methods are not equivalent.

e. Kruskal-Wallis Test

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

The Kruskal-Wallis test does not assume that the data is normally distributed, and works with the ranks of the observations. The test statistic is about 21.2 with 3 degrees of freedom and a p-value of about 0.0001, which is smaller than alpha = 0.05. So we reject the null hypothesis, and the nonparametric test agrees with the parametric test on the transformed data: the four estimation methods do not produce equivalent estimates of peak discharge.

Complete Code

m1<-c(0.34,0.12,1.23,0.70,1.75,0.12)
m2<-c(0.91,2.94,2.14,2.36,2.86,4.55)
m3<-c(6.31,8.37,9.75,6.09,9.82,7.24)
m4<-c(17.15,11.82,10.97,17.20,14.35,16.82)
obs<-c(m1,m2,m3,m4)
x<-as.factor(c(rep(1,6),rep(2,6),rep(3,6),rep(4,6)))
cat("sample mean of the peak discharge of each method \n",c(mean(m1),mean(m2),mean(m3),mean(m4)),"\n")
cat("sample standard deviation of the peak discharge of each method \n",c(sd(m1),sd(m2),sd(m3),sd(m4)),"\n")
boxplot(obs~x,xlab="Estimation Method",ylab="Peak Discharge (cfs)",main="Boxplot of Observations")
model<-aov(obs~x)
summary(model)
meanx<-c(rep(mean(m1),6),rep(mean(m2),6),rep(mean(m3),6),rep(mean(m4),6))
res<-obs-meanx
qqnorm(res)
plot(meanx,res,xlab="method average",ylab="residual",main="constant variance wrt method means")
library(MASS)
boxcox(obs~x)
lambda=0.5
obs2<-obs^(lambda)
boxcox(obs2~x)
boxplot(obs2~x,xlab="Estimation Method",ylab="Square Root of Peak Discharge",main="Boxplot of Transformed Observations")
model2<-aov(obs2~x)
summary(model2)
meanx2<-c(rep(mean(m1^lambda),6),rep(mean(m2^lambda),6),rep(mean(m3^lambda),6),rep(mean(m4^lambda),6))
res2<-obs2-meanx2
qqnorm(res2)
plot(meanx2,res2,xlab="method average",ylab="residual",main="constant variance wrt method means")
kruskal.test(obs~x)