First Problem

We use pwr.anova.test() from the pwr package with k = 4 fluids, α = 0.05, and power = 0.80. The variance of fluid life is estimated as 3.5, so the standard deviation is σ = √3.5 ≈ 1.871 hrs. The effect size is Cohen’s f = σ_means / σ, where σ_means = √( Σ(μᵢ − μ̄)² / k ).

The smallest mean is 18 hrs and the largest is 20 hrs (a difference of 2 hrs). The remaining two means are placed according to each variability pattern:

Variability Means (hrs) Effect size f n (exact) n per fluid
Minimum 18, 19, 19, 20 0.3780 20.08 21
Intermediate 18, 18.67, 19.33, 20 0.3984 18.18 19
Maximum 18, 18, 20, 20 0.5345 10.57 11

Rounding up, 21 samples of each fluid are needed for minimum variability, 19 samples for intermediate variability, and 11 samples for maximum variability. Minimum variability is the most conservative case, so collecting 21 samples per fluid would satisfy the design criterion no matter how the means are spread between 18 and 20 hrs.

Second Problem

The data are entered in a wide format and then converted to a tidy format using pivot_longer.

# A tibble: 8 × 2
  Fluid  Life
  <fct> <dbl>
1 1      17.6
2 2      16.9
3 3      21.4
4 4      19.3
5 1      18.9
6 2      15.3
7 3      23.6
8 4      21.1
  Fluid    Mean Variance
1     1 18.6500   3.8110
2     2 17.9500   3.4390
3     3 20.9500   3.5310
4     4 18.8167   2.4177

a)

  • H₀: μ₁ = μ₂ = μ₃ = μ₄
  • H₁: at least one μᵢ differs from the others
            Df Sum Sq Mean Sq F value Pr(>F)  
Fluid        3  30.17   10.05   3.047 0.0525 .
Residuals   20  65.99    3.30                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The one-way ANOVA gives F₀ = 3.047 with a p-value of 0.0525. Since the p-value (0.0525) < α = 0.10, we reject H₀ and conclude that the mean life is not the same for all four fluids at the 0.10 level of significance. (This result would not be significant at α = 0.05, so the conclusion depends on the chosen α.)

b)

Model adequacy is checked with the normal probability plot of the residuals, the residuals vs. fitted values plot, and a side-by-side boxplot of the four fluids.


    Shapiro-Wilk normality test

data:  residuals(model)
W = 0.95671, p-value = 0.376
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3   0.137 0.9368
      20

In the normal probability plot, the residuals fall close to the reference line with no strong curvature or outliers, so the normality assumption is reasonable. This is supported by the Shapiro-Wilk test (p = 0.376).

In the residuals vs. fitted plot, the residuals are scattered randomly around zero and the vertical spread is similar at all four fitted values, with no funnel shape. The boxplots also show similar spreads, and the sample variances range only from 2.42 to 3.81. The modified Levene test (p = 0.937) does not reject equal variances, so the constant variance assumption holds.

Assuming the runs were performed in random order (as in a CRD), the independence assumption is also reasonable. Therefore, the model is adequate.

c)

  Tukey multiple comparisons of means
    90% family-wise confidence level

Fit: aov(formula = Life ~ Fluid, data = dat)

$Fluid
          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.2090635

Using Tukey’s test with a familywise error rate of α = 0.10, only the pair Fluid 3 – Fluid 2 has a confidence interval that does not contain zero (difference = 3.00 hrs, 90% CI [0.433, 5.567], adjusted p = 0.044). All other pairwise intervals contain zero, so those differences are not significant. Therefore, fluid 3 has a significantly longer mean life than fluid 2, and no other pairs of fluids differ significantly.

R Code

library(pwr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(car)

# Problem 1: sample size
sigma <- sqrt(3.5)

effect_f <- function(means, sigma) {
  sqrt(sum((means - mean(means))^2) / length(means)) / sigma
}

mu_min <- c(18, 19, 19, 20)
mu_int <- seq(18, 20, length.out = 4)
mu_max <- c(18, 18, 20, 20)

f_min <- effect_f(mu_min, sigma)
f_int <- effect_f(mu_int, sigma)
f_max <- effect_f(mu_max, sigma)

p_min <- pwr.anova.test(k = 4, f = f_min, sig.level = 0.05, power = 0.80)
p_int <- pwr.anova.test(k = 4, f = f_int, sig.level = 0.05, power = 0.80)
p_max <- pwr.anova.test(k = 4, f = f_max, sig.level = 0.05, power = 0.80)

power_table <- data.frame(
  Variability = c("Minimum", "Intermediate", "Maximum"),
  Means = c(paste(round(mu_min, 2), collapse = ", "),
            paste(round(mu_int, 2), collapse = ", "),
            paste(round(mu_max, 2), collapse = ", ")),
  Effect_size_f = round(c(f_min, f_int, f_max), 4),
  n_exact = round(c(p_min$n, p_int$n, p_max$n), 2),
  n_per_fluid = ceiling(c(p_min$n, p_int$n, p_max$n))
)

knitr::kable(power_table, col.names = c("Variability", "Means (hrs)", "Effect size f", "n (exact)", "n per fluid"))

# Problem 2: data in tidy format
wide <- data.frame(
  Fluid_1 = c(17.6, 18.9, 16.3, 17.4, 20.1, 21.6),
  Fluid_2 = c(16.9, 15.3, 18.6, 17.1, 19.5, 20.3),
  Fluid_3 = c(21.4, 23.6, 19.4, 18.5, 20.5, 22.3),
  Fluid_4 = c(19.3, 21.1, 16.9, 17.5, 18.3, 19.8)
)

dat <- wide %>%
  pivot_longer(cols = everything(), names_to = "Fluid", names_prefix = "Fluid_", values_to = "Life")
dat$Fluid <- as.factor(dat$Fluid)

head(dat, 8)

dat %>%
  group_by(Fluid) %>%
  summarise(Mean = round(mean(Life), 4), Variance = round(var(Life), 4)) %>%
  as.data.frame()

# a) One-way ANOVA
model <- aov(Life ~ Fluid, data = dat)
summary(model)

# b) Model adequacy
res <- data.frame(Fitted = fitted(model), Residuals = residuals(model))

ggplot(res, aes(sample = Residuals)) + 
  stat_qq(color = "blue") +                              
  stat_qq_line() +
  labs(x = "Theoretical Quantiles", y = "Residuals", title = "Normal Probability Plot of Residuals")

ggplot(res, aes(x = Fitted, y = Residuals)) +
  geom_point(color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Fitted Values", y = "Residuals", title = "Residuals vs Fitted Values")

boxplot(Life ~ Fluid, data = dat,
        main = "Side by Side Boxplot: Fluid Life at 35kV",
        xlab = "Fluid Type",
        ylab = "Life (hr)",
        col = c("blue", "red", "green", "orange"))

shapiro.test(residuals(model))
leveneTest(Life ~ Fluid, data = dat)

# c) Tukey's test (familywise alpha = 0.10)
tukey <- TukeyHSD(model, conf.level = 0.90)
tukey

par(mar = c(5, 6, 4, 2))
plot(tukey, las = 1)