Assignment 10 - VST and Kruskal-Wallace
G M Farhan Tanvir
Assignment 10
VST and Kruskal-Wallace Using R
Data
fluid1 <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12)
fluid2 <- c(0.91, 2.94, 2.14, 2.36, 2.86, 4.55)
fluid3 <- c(6.31, 8.37, 9.75, 6.09, 9.82, 7.24)
fluid4 <- c(17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
response <- c(fluid1, fluid2, fluid3, fluid4)
method <- factor(c(
rep(1, 6),
rep(2, 6),
rep(3, 6),
rep(4, 6)
))
dat <- data.frame(response, method)
head(dat)## response method
## 1 0.34 1
## 2 0.12 1
## 3 1.23 1
## 4 0.70 1
## 5 1.75 1
## 6 0.12 11. Linear Effects Equation and Hypothesis
The linear effects model is
\[ Y_{ij} = \mu + \tau_i + \epsilon_{ij} \]
where \(\mu\) is the overall mean, \(\tau_i\) is the effect of estimation method, and \(\epsilon_{ij}\) is the random error.
The hypotheses are
\[ H_0:\mu_1=\mu_2=\mu_3=\mu_4 \]
\[ H_a:\text{at least one population mean differs} \]
2. Normality and Constant Variance
Normal Probability Plot
par(mfrow=c(2,2))
qqnorm(fluid1, main="Method 1")
qqline(fluid1)
qqnorm(fluid2, main="Method 2")
qqline(fluid2)
qqnorm(fluid3, main="Method 3")
qqline(fluid3)
qqnorm(fluid4, main="Method 4")
qqline(fluid4)
par(mfrow=c(1,1))The normal probability plots show some departures from normality, particularly for the individual groups. The data therefore do not appear to be strongly normally distributed.
Side-by-Side Boxplots
boxplot(
response ~ method,
data = dat,
main = "Flood Discharge by Estimation Method",
xlab = "Estimation Method",
ylab = "Peak Discharge"
)
The boxplots show that the variability differs substantially between methods. The higher-mean methods also have noticeably larger spreads. Thus, the constant variance assumption appears questionable for the raw data.
3. Fit the Raw-Data ANOVA Model
model_raw <- aov(response ~ method, data = dat)
summary(model_raw)## 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 ' ' 1Residual Diagnostics
par(mfrow=c(1,2))
plot(model_raw, which=1,
main="Residuals vs Fitted")
plot(model_raw, which=2,
main="Normal Q-Q Plot")
par(mfrow=c(1,1))The residual plots show evidence of unequal variability, particularly because the spread of the residuals increases with the fitted values. Therefore, the raw-data ANOVA model does not satisfy the constant-variance assumption very well.
4. Box-Cox Transformation
library(MASS)
bc <- boxcox(
lm(response ~ method, data = dat),
lambda = seq(-2, 2, by=0.01),
plotit = TRUE
)
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 0.54The Box-Cox analysis gives an estimated transformation parameter of approximately
\[ \lambda \approx 0.54 \]
Therefore, an appropriate transformation is
\[ Y^*=\frac{Y^{0.54}-1}{0.54} \]
lambda <- 0.54
dat$transformed <- (dat$response^lambda - 1) / lambda
head(dat)## response method transformed
## 1 0.34 1 -0.8176511
## 2 0.12 1 -1.2625143
## 3 1.23 1 0.2190285
## 4 0.70 1 -0.3244294
## 5 1.75 1 0.6533734
## 6 0.12 1 -1.2625143ANOVA Using the Transformed Data
model_trans <- aov(transformed ~ method, data = dat)
summary(model_trans)## Df Sum Sq Mean Sq F value Pr(>F)
## method 3 147.22 49.07 83.43 1.76e-11 ***
## Residuals 20 11.76 0.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The transformed-data ANOVA gives a very small p-value, so the null hypothesis is rejected at the 0.05 significance level.
Residual Diagnostics After Transformation
par(mfrow=c(1,2))
plot(model_trans, which=1,
main="Transformed Residuals vs Fitted")
plot(model_trans, which=2,
main="Transformed Normal Q-Q Plot")
par(mfrow=c(1,1))After the transformation, the residual spread is much more consistent across fitted values and the normal Q-Q plot is improved. Therefore, the transformed model better satisfies the ANOVA assumptions.
The hypothesis test is
\[ H_0:\mu_1=\mu_2=\mu_3=\mu_4 \]
versus
\[ H_a:\text{at least one mean differs} \]
summary(model_trans)## Df Sum Sq Mean Sq F value Pr(>F)
## method 3 147.22 49.07 83.43 1.76e-11 ***
## Residuals 20 11.76 0.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Because the p-value is less than 0.05, we reject \(H_0\).
Conclusion: There is statistically significant evidence that the mean flood discharge differs among the four estimation methods.
5. Kruskal-Wallis Test
The Kruskal-Wallis test is a nonparametric alternative to one-way ANOVA.
kruskal.test(response ~ method, data = dat)##
## Kruskal-Wallis rank sum test
##
## data: response by method
## Kruskal-Wallis chi-squared = 21.156, df = 3, p-value = 9.771e-05The Kruskal-Wallis test gives a p-value much less than 0.05.
Therefore,
\[ \boxed{\text{Reject }H_0} \]
Conclusion: There is sufficient evidence at the 0.05 significance level that the distributions of peak discharge are not the same for all four estimation methods.
Summary
The raw data show evidence of unequal variance and some departures from normality. A Box-Cox transformation with a parameter of approximately 0.54 substantially improves the ANOVA assumptions. Both the transformed parametric ANOVA and the nonparametric Kruskal-Wallis test indicate a significant difference among the four estimation methods.
Complete Code
fluid1 <- c(0.34, 0.12, 1.23, 0.70, 1.75, 0.12)
fluid2 <- c(0.91, 2.94, 2.14, 2.36, 2.86, 4.55)
fluid3 <- c(6.31, 8.37, 9.75, 6.09, 9.82, 7.24)
fluid4 <- c(17.15, 11.82, 10.97, 17.20, 14.35, 16.82)
response <- c(fluid1, fluid2, fluid3, fluid4)
method <- factor(c(
rep(1, 6),
rep(2, 6),
rep(3, 6),
rep(4, 6)
))
dat <- data.frame(response, method)
par(mfrow=c(2,2))
qqnorm(fluid1, main="Method 1")
qqline(fluid1)
qqnorm(fluid2, main="Method 2")
qqline(fluid2)
qqnorm(fluid3, main="Method 3")
qqline(fluid3)
qqnorm(fluid4, main="Method 4")
qqline(fluid4)
par(mfrow=c(1,1))
boxplot(
response ~ method,
data = dat,
main = "Flood Discharge by Estimation Method",
xlab = "Estimation Method",
ylab = "Peak Discharge"
)
model_raw <- aov(response ~ method, data = dat)
summary(model_raw)
par(mfrow=c(1,2))
plot(model_raw, which=1,
main="Residuals vs Fitted")
plot(model_raw, which=2,
main="Normal Q-Q Plot")
par(mfrow=c(1,1))
library(MASS)
bc <- boxcox(
lm(response ~ method, data = dat),
lambda = seq(-2, 2, by=0.01),
plotit = TRUE
)
lambda <- bc$x[which.max(bc$y)]
lambda
lambda <- 0.54
dat$transformed <- (dat$response^lambda - 1) / lambda
model_trans <- aov(transformed ~ method, data = dat)
summary(model_trans)
par(mfrow=c(1,2))
plot(model_trans, which=1,
main="Transformed Residuals vs Fitted")
plot(model_trans, which=2,
main="Transformed Normal Q-Q Plot")
par(mfrow=c(1,1))
kruskal.test(response ~ method, data = dat)