1 Data Input

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.

The data will be put into a data frame using the code chunk below.

method<-c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4)
obs<-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)
dat<-data.frame(method, obs)
dat$method<-as.factor(dat$method)

2 Part A

Write the linear effects equation and the hypothesis you are testing.

The linear effects equation shows the observation equivalent to the grand mean, the effect of the ith population, and the random error associated with the measurement. This equation can be seen below.

\[ x_{ij}=\mu+\tau_i+\epsilon_{ij} \]

The null hypothesis for this experiment is that the effects for the individual populations is equal to zero and does not affect the grand mean. The equation for the null hypothesis is \(H_o:\tau_i=0\) for all i populations.

The alternative hypothesis for this experiment is that at least one of the effects for the individual populations is not equal to zero and does have an affect on the grand mean. The equation for the alternative hypothesis is \(H_1:\tau_i\neq0\) for at least one of the i populations.

3 Part B

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

The data will be checked for normality and constant variance, which are the two strong assumptions for a valid analysis of variance test. The plots can be seen in the code chunk below.

qqnorm(dat$obs)

hist(dat$obs)

boxplot(obs~method, data = dat)

Based on the normal probability plot, histogram, and box plots, the sampled data does not follow a normal distribution and does not have constant variance. The normal probability plot has some curvature and the histogram indicates a right-skewed distribution. The box plots vary widely in size leaving no doubt that this data set does not have constant variance.

4 Part C

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

Even though the data does not meet the assumptions, we will still run an analysis of variance to see how the test handles violations of the assumptions. The test can be seen in the code chunk below.

anova<-aov(obs~method, data = dat)
summary(anova)
##             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
plot(anova, 1)

plot(anova, 2)

The analysis returned a massive f statistic of 76.29 and p-value of 4e-11. Also, the mean squared treatment is significantly larger than the mean squared error. Since the populations have unequal variances (clearly seen in the Residuals vs Fitted Plot), this could be invalidating the results of the ANOVA test.

5 Part D

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

Since the populations have unequal variances and do not follow a normal distribution, a transformation can be performed to attempt to get the data to meet those assumptions. A Box Cox transformation will be performed and can be seen in the code chunk below.

#install.packages("MASS")
library(MASS)
boxcox(obs~method, data = dat)

lambda = 0.60
x1<-dat$obs^lambda
anova1<-aov(x1~method, data = dat)
summary(anova1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## method       3  63.71  21.236   85.76 1.36e-11 ***
## Residuals   20   4.95   0.248                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(anova1, 1)

plot(anova1, 2)

This analysis yielded an f statistic of 85.76 and p-value of 1.36e-11, which means we would reject the null hypothesis of the means being equal. The transformation brought the data much closer to constant variance and normality but still varies slightly. It is possible that a non-parametric test would better suit this data set.

6 Part E

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

We will now perform a non-parametric test, the Kruskal-Wallis rank sum test, on the data set in the code chunk below.

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

The Kruskal-Wallis test yielded a p-value of 9.771e-5, which means we reject the null hypothesis that the means of the populations are equal at the 95% confidence level.

7 Complete R Code

# Data Input
method<-c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4)
obs<-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)
dat<-data.frame(method, obs)
dat$method<-as.factor(dat$method)

# Part B
qqnorm(dat$obs)
hist(dat$obs)
boxplot(obs~method, data = dat)

# Part C
anova<-aov(obs~method, data = dat)
summary(anova)
plot(anova, 1)
plot(anova, 2)

# Part D
#install.packages("MASS")
library(MASS)
boxcox(obs~method, data = dat)
lambda = 0.60
x1<-dat$obs^lambda
anova1<-aov(x1~method, data = dat)
summary(anova1)
plot(anova1, 1)
plot(anova1, 2)

# Part E
kruskal.test(obs~method, data = dat)