Assignment 9 - CRD/One-Way ANOVA and Multiple Comparisons in R
G M Farhan Tanvir
2026-09-29
Question 1: Sample Size Determination
We wish to design a new experiment to test for a significant difference between the mean effective life of four insulating fluids at an accelerated load of 35kV. The minimum mean life is 18 hours and the maximum mean life is 20 hours. The desired type I error probability is 0.05 and the desired power is 0.80.
# Number of populations
pop_k <- 4
# Difference between maximum and minimum means
d <- 20 - 18
# Minimum variability case
min_f <- d * sqrt(1/(2*pop_k))
# Intermediate variability case
intermediate_f <- (d/2) * sqrt((pop_k+1)/(3*(pop_k-1)))
# Maximum variability case
if (pop_k %% 2 == 0) {
max_f <- d/2
} else {
max_f <- d * sqrt(pop_k^2 - 1)/(2*pop_k)
}
# Power calculations
min_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = min_f,
sig.level = 0.05,
power = 0.80
)
intermediate_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = intermediate_f,
sig.level = 0.05,
power = 0.80
)
max_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = max_f,
sig.level = 0.05,
power = 0.80
)
min_result##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 6.515965
## f = 0.7071068
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each groupintermediate_result##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 5.979231
## f = 0.745356
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each groupmax_result##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 3.856403
## f = 1
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each groupResults
The required sample sizes are:
- Minimum variability: 7 observations per fluid
- Intermediate variability: 6 observations per fluid
- Maximum variability: 4 observations per fluid
Question 2: One-Way ANOVA
The experimenter collected six observations for each type of fluid.
# Enter the response data
fluid1 <- c(17.6, 18.9, 16.3, 17.4, 20.1, 21.6)
fluid2 <- c(16.9, 15.3, 18.6, 17.1, 19.5, 20.3)
fluid3 <- c(21.4, 23.6, 19.4, 18.5, 20.5, 22.3)
fluid4 <- c(19.3, 21.1, 16.9, 17.5, 18.3, 19.8)
# Create a tidy data frame
fluid.type <- factor(c(
rep(1,6),
rep(2,6),
rep(3,6),
rep(4,6)
))
response <- c(fluid1, fluid2, fluid3, fluid4)
df <- data.frame(response, fluid.type)
df## response fluid.type
## 1 17.6 1
## 2 18.9 1
## 3 16.3 1
## 4 17.4 1
## 5 20.1 1
## 6 21.6 1
## 7 16.9 2
## 8 15.3 2
## 9 18.6 2
## 10 17.1 2
## 11 19.5 2
## 12 20.3 2
## 13 21.4 3
## 14 23.6 3
## 15 19.4 3
## 16 18.5 3
## 17 20.5 3
## 18 22.3 3
## 19 19.3 4
## 20 21.1 4
## 21 16.9 4
## 22 17.5 4
## 23 18.3 4
## 24 19.8 42a. One-Way ANOVA Test
The hypotheses are:
\[ H_0:\mu_1=\mu_2=\mu_3=\mu_4 \]
\[ H_a:\text{At least one population mean differs} \]
The significance level is:
\[ \alpha=0.10 \]
# Fit the one-way ANOVA model
model <- aov(response ~ fluid.type, data = df)
# ANOVA table
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## fluid.type 3 30.16 10.05 3.047 0.0525 .
## Residuals 20 65.99 3.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# One-way ANOVA with additional digits
fluid_compare <- oneway.test(
response ~ fluid.type,
data = df,
var.equal = TRUE
)
fluid_compare##
## One-way analysis of means
##
## data: response and fluid.type
## F = 3.0473, num df = 3, denom df = 20, p-value = 0.05246The test gives an F-statistic of approximately 3.0473 and a p-value of approximately 0.05246.
Since:
\[ 0.05246 < 0.10 \]
we reject the null hypothesis.
Conclusion: There is sufficient evidence at the 0.10 significance level that the mean life of the four fluids is not the same.
2b. Model Adequacy
Normal Probability Plot of Residuals
plot(model, 2)
The residuals appear approximately linear on the normal probability plot, so the normality assumption appears reasonable.
Normal Probability Plots by Fluid
qqnorm(
fluid1,
main = "Normal Probability Plot - Fluid 1",
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles"
)
qqline(fluid1)
qqnorm(
fluid2,
main = "Normal Probability Plot - Fluid 2",
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles"
)
qqline(fluid2)
qqnorm(
fluid3,
main = "Normal Probability Plot - Fluid 3",
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles"
)
qqline(fluid3)
qqnorm(
fluid4,
main = "Normal Probability Plot - Fluid 4",
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles"
)
qqline(fluid4)
The observations in the individual normal probability plots are reasonably close to the reference lines, supporting the normality assumption.
Side-by-Side Boxplot
boxplot(
response ~ fluid.type,
data = df,
main = "Life of Insulating Fluids",
xlab = "Fluid Type",
ylab = "Life (hours)"
)
The boxplots show that the spreads of the four fluid types are reasonably similar. Therefore, the constant variance assumption appears reasonable.
Overall, the ANOVA model appears adequate because the residuals are approximately normal and the treatment variances appear reasonably constant.
2c. Tukey Multiple Comparisons
Because the ANOVA null hypothesis was rejected, Tukey’s test is used to determine which fluid means differ.
The familywise error rate is:
\[ \alpha=0.10 \]
Therefore, the confidence level is:
\[ 1-\alpha=0.90 \]
tukey_result <- TukeyHSD(
model,
conf.level = 0.90
)
tukey_result## Tukey multiple comparisons of means
## 90% family-wise confidence level
##
## Fit: aov(formula = response ~ fluid.type, data = df)
##
## $fluid.type
## diff lwr upr p adj
## 2-1 -0.7000000 -3.2670196 1.8670196 0.9080815
## 3-1 2.3000000 -0.2670196 4.8670196 0.1593262
## 4-1 0.1666667 -2.4003529 2.7336862 0.9985213
## 3-2 3.0000000 0.4329804 5.5670196 0.0440578
## 4-2 0.8666667 -1.7003529 3.4336862 0.8413288
## 4-3 -2.1333333 -4.7003529 0.4336862 0.2090635Tukey Confidence Interval Plot
plot(tukey_result)
A pair of means is significantly different when its confidence interval does not contain zero.
The Tukey test shows that:
\[ \boxed{\text{Fluid 3 and Fluid 2 significantly differ}} \]
The other pairs do not significantly differ at the 0.10 familywise error rate.
Complete R Code
pop_k <- 4
d <- 20 - 18
min_f <- d * sqrt(1/(2*pop_k))
intermediate_f <- (d/2) * sqrt((pop_k+1)/(3*(pop_k-1)))
if (pop_k %% 2 == 0) {
max_f <- d/2
} else {
max_f <- d * sqrt(pop_k^2 - 1)/(2*pop_k)
}
min_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = min_f,
sig.level = 0.05,
power = 0.80
)
intermediate_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = intermediate_f,
sig.level = 0.05,
power = 0.80
)
max_result <- pwr::pwr.anova.test(
k = pop_k,
n = NULL,
f = max_f,
sig.level = 0.05,
power = 0.80
)
min_result
intermediate_result
max_result
fluid1 <- c(17.6, 18.9, 16.3, 17.4, 20.1, 21.6)
fluid2 <- c(16.9, 15.3, 18.6, 17.1, 19.5, 20.3)
fluid3 <- c(21.4, 23.6, 19.4, 18.5, 20.5, 22.3)
fluid4 <- c(19.3, 21.1, 16.9, 17.5, 18.3, 19.8)
fluid.type <- factor(c(
rep(1,6),
rep(2,6),
rep(3,6),
rep(4,6)
))
response <- c(fluid1, fluid2, fluid3, fluid4)
df <- data.frame(response, fluid.type)
model <- aov(response ~ fluid.type, data = df)
summary(model)
fluid_compare <- oneway.test(
response ~ fluid.type,
data = df,
var.equal = TRUE
)
fluid_compare
plot(model, 2)
qqnorm(fluid1)
qqline(fluid1)
qqnorm(fluid2)
qqline(fluid2)
qqnorm(fluid3)
qqline(fluid3)
qqnorm(fluid4)
qqline(fluid4)
boxplot(
response ~ fluid.type,
data = df,
main = "Life of Insulating Fluids",
xlab = "Fluid Type",
ylab = "Life (hours)"
)
tukey_result <- TukeyHSD(
model,
conf.level = 0.90
)
tukey_result
plot(tukey_result)