1 Problem 1: Silver powder production

1.1 Introduction

This experiment investigates the effects of ammonium concentration, stir rate, and temperature on silver powder density.

The three factors are fixed effects, each with two levels. Every treatment combination has two observations, giving a total of 16 observations.

The significance level is \(\alpha=0.05\).

1.2 Data

powder <- data.frame(
  Ammonium = factor(
    c(2, 2, 30, 30,
      2, 2, 30, 30,
      2, 2, 30, 30,
      2, 2, 30, 30),
    levels = c(2, 30)
  ),

  StirRate = factor(
    c(100, 100, 100, 100,
      150, 150, 150, 150,
      100, 100, 100, 100,
      150, 150, 150, 150),
    levels = c(100, 150)
  ),

  Temperature = factor(
    c(8, 8, 8, 8,
      8, 8, 8, 8,
      40, 40, 40, 40,
      40, 40, 40, 40),
    levels = c(8, 40)
  ),

  Density = c(
    14.68, 15.18,
    15.12, 17.48,
    7.54, 6.66,
    12.46, 12.62,
    10.95, 17.68,
    12.65, 15.96,
    8.03, 8.84,
    14.96, 14.96
  )
)

knitr::kable(
  powder,
  caption = "Silver powder production data"
)
Silver powder production data
Ammonium StirRate Temperature Density
2 100 8 14.68
2 100 8 15.18
30 100 8 15.12
30 100 8 17.48
2 150 8 7.54
2 150 8 6.66
30 150 8 12.46
30 150 8 12.62
2 100 40 10.95
2 100 40 17.68
30 100 40 12.65
30 100 40 15.96
2 150 40 8.03
2 150 40 8.84
30 150 40 14.96
30 150 40 14.96 
stopifnot(nrow(powder) == 16)

with(
  powder,
  table(Ammonium, StirRate, Temperature)
)
## , , Temperature = 8
## 
##         StirRate
## Ammonium 100 150
##       2    2   2
##       30   2   2
## 
## , , Temperature = 40
## 
##         StirRate
## Ammonium 100 150
##       2    2   2
##       30   2   2

1.3 Part (a): Full factorial model

Let \(A\) denote ammonium concentration, \(B\) denote stir rate, and \(C\) denote temperature.

The full factorial model is

\[ \begin{aligned} Y_{ijkl}={}&\mu+\alpha_i+\beta_j+\gamma_k\\ &+(\alpha\beta)_{ij} +(\alpha\gamma)_{ik} +(\beta\gamma)_{jk}\\ &+(\alpha\beta\gamma)_{ijk} +\varepsilon_{ijkl}. \end{aligned} \]

Here:

  • \(\mu\) is the overall mean.
  • \(\alpha_i\), \(\beta_j\), and \(\gamma_k\) are the main effects.
  • The two-factor terms represent interactions between pairs of factors.
  • The three-factor term represents the interaction among all three factors.
  • \(\varepsilon_{ijkl}\) is the random error for observation \(l\) in treatment combination \(i,j,k\).

The factor effects are fixed. The errors are assumed independent, normally distributed, with mean zero and common variance \(\sigma^2\).

1.4 Part (b): ANOVA and significant factors

Fit the full factorial model, including all main effects, two-factor interactions, and the three-factor interaction.

q1_model <- aov(
  Density ~ Ammonium * StirRate * Temperature,
  data = powder
)

summary(q1_model)
##                               Df Sum Sq Mean Sq F value  Pr(>F)   
## Ammonium                       1  44.39   44.39  11.180 0.01018 * 
## StirRate                       1  70.69   70.69  17.804 0.00292 **
## Temperature                    1   0.33    0.33   0.083 0.78117   
## Ammonium:StirRate              1  28.12   28.12   7.082 0.02875 * 
## Ammonium:Temperature           1   0.02    0.02   0.005 0.94281   
## StirRate:Temperature           1  10.13   10.13   2.551 0.14889   
## Ammonium:StirRate:Temperature  1   1.52    1.52   0.383 0.55341   
## Residuals                      8  31.76    3.97                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display the ANOVA results in a table.

q1_anova <- anova(q1_model)

q1_results <- data.frame(
  Effect = rownames(q1_anova),
  Df = q1_anova$Df,
  Sum_Squares = q1_anova$`Sum Sq`,
  Mean_Square = q1_anova$`Mean Sq`,
  F_value = q1_anova$`F value`,
  P_value = q1_anova$`Pr(>F)`,
  row.names = NULL
)

knitr::kable(
  q1_results,
  digits = 5,
  caption = "Full factorial ANOVA for silver powder density"
)
Full factorial ANOVA for silver powder density
Effect Df Sum_Squares Mean_Square F_value P_value
Ammonium 1 44.38891 44.38891 11.18026 0.01018
StirRate 1 70.68606 70.68606 17.80373 0.00292
Temperature 1 0.32776 0.32776 0.08255 0.78117
Ammonium:StirRate 1 28.11651 28.11651 7.08172 0.02875
Ammonium:Temperature 1 0.02176 0.02176 0.00548 0.94281
StirRate:Temperature 1 10.12831 10.12831 2.55102 0.14889
Ammonium:StirRate:Temperature 1 1.51906 1.51906 0.38261 0.55341
Residuals 8 31.76235 3.97029 NA NA

An effect is statistically significant when its p-value is below 0.05.

q1_significant <- subset(
  q1_results,
  !is.na(P_value) & P_value < 0.05
)

knitr::kable(
  q1_significant,
  digits = 5,
  caption = "Statistically significant effects at alpha = 0.05"
)
Statistically significant effects at alpha = 0.05
Effect Df Sum_Squares Mean_Square F_value P_value
1 Ammonium 1 44.38891 44.38891 11.18026 0.01018
2 StirRate 1 70.68606 70.68606 17.80373 0.00292
4 Ammonium:StirRate 1 28.11651 28.11651 7.08172 0.02875

The following code reports the conclusion for each effect directly from the computed ANOVA results.

q1_results$Conclusion <- ifelse(
  is.na(q1_results$P_value),
  "Not tested",
  ifelse(
    q1_results$P_value < 0.05,
    "Significant",
    "Not significant"
  )
)

knitr::kable(
  q1_results[, c("Effect", "P_value", "Conclusion")],
  digits = 5,
  caption = "Significance decisions for the full factorial model"
)
Significance decisions for the full factorial model
Effect P_value Conclusion
Ammonium 0.01018 Significant
StirRate 0.00292 Significant
Temperature 0.78117 Not significant
Ammonium:StirRate 0.02875 Significant
Ammonium:Temperature 0.94281 Not significant
StirRate:Temperature 0.14889 Not significant
Ammonium:StirRate:Temperature 0.55341 Not significant
Residuals NA Not tested

1.5 Treatment means

Calculate the mean powder density at each combination of factor levels.

q1_means <- aggregate(
  Density ~ Ammonium + StirRate + Temperature,
  data = powder,
  FUN = mean
)

knitr::kable(
  q1_means,
  digits = 4,
  caption = "Mean density by treatment combination"
)
Mean density by treatment combination
Ammonium StirRate Temperature Density
2 100 8 14.930
30 100 8 16.300
2 150 8 7.100
30 150 8 12.540
2 100 40 14.315
30 100 40 14.305
2 150 40 8.435
30 150 40 14.960

1.6 Interaction plot: ammonium and stir rate

The ammonium-by-stir-rate interaction is examined at each temperature level.

library(ggplot2)

q1_means_plot <- aggregate(
  Density ~ Ammonium + StirRate + Temperature,
  data = powder,
  FUN = mean
)

ggplot(
  q1_means_plot,
  aes(
    x = StirRate,
    y = Density,
    color = Ammonium,
    group = Ammonium
  )
) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 3) +
  facet_wrap(~ Temperature, labeller = label_both) +
  labs(
    title = "Ammonium concentration and stir-rate interaction",
    x = "Stir rate (RPM)",
    y = "Mean powder density",
    color = "Ammonium (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Nonparallel lines suggest that the effect of ammonium concentration changes with stir rate. The ANOVA p-value for the corresponding interaction determines whether the interaction is statistically significant at the 0.05 level.

The interaction plot should be interpreted together with the full factorial ANOVA, including the three-factor interaction.

2 Problem 2: Carbon-anode baked density

2.1 Introduction

This experiment investigates whether firing temperature and furnace position affect the baked density of a carbon anode.

The factors are:

  • Position: 1 and 2.
  • Temperature: 800°C, 825°C, and 850°C.

There are three observations per treatment combination, giving 18 observations.

The analysis considers three models:

  1. Position and temperature are both fixed effects.
  2. Position and temperature are both random effects.
  3. Position is fixed and temperature is random.

2.2 Data

The data below reproduce the observations in the assignment table.

anode <- data.frame(
  Position = factor(
    rep(c(1, 2), each = 9),
    levels = c(1, 2)
  ),

  Temperature = factor(
    rep(
      rep(c(800, 825, 850), each = 3),
      times = 2
    ),
    levels = c(800, 825, 850)
  ),

  Density = c(
    # Position 1
    570, 565, 583,       # 800 degrees
    1063, 1080, 1043,    # 825 degrees
    565, 510, 590,       # 850 degrees

    # Position 2
    528, 547, 521,       # 800 degrees
    988, 1026, 1004,     # 825 degrees
    526, 538, 532        # 850 degrees
  )
)

stopifnot(nrow(anode) == 18)

knitr::kable(
  anode,
  caption = "Carbon-anode baked-density data"
)
Carbon-anode baked-density data
Position Temperature Density
1 800 570
1 800 565
1 800 583
1 825 1063
1 825 1080
1 825 1043
1 850 565
1 850 510
1 850 590
2 800 528
2 800 547
2 800 521
2 825 988
2 825 1026
2 825 1004
2 850 526
2 850 538
2 850 532 
with(
  anode,
  table(Position, Temperature)
)
##         Temperature
## Position 800 825 850
##        1   3   3   3
##        2   3   3   3

2.3 Factorial model

The two-factor model is

\[ Y_{ijk} = \mu+P_i+T_j+(PT)_{ij}+\varepsilon_{ijk}, \]

where:

  • \(\mu\) is the overall mean.
  • \(P_i\) is the position effect.
  • \(T_j\) is the temperature effect.
  • \((PT)_{ij}\) is the position-by-temperature interaction.
  • \(\varepsilon_{ijk}\) is the experimental error.

The appropriate F-test denominator depends on the assumed fixed or random nature of the factors.

2.4 Preliminary ANOVA components

Fit the two-factor model and extract the sums of squares, degrees of freedom, and mean squares.

q2_model <- aov(
  Density ~ Position * Temperature,
  data = anode
)

q2_anova <- anova(q2_model)

q2_components <- data.frame(
  Effect = rownames(q2_anova),
  Df = q2_anova$Df,
  Sum_Squares = q2_anova$`Sum Sq`,
  Mean_Square = q2_anova$`Mean Sq`,
  F_value = q2_anova$`F value`,
  P_value = q2_anova$`Pr(>F)`,
  row.names = NULL
)

knitr::kable(
  q2_components,
  digits = 6,
  caption = "ANOVA components for the carbon-anode experiment"
)
ANOVA components for the carbon-anode experiment
Effect Df Sum_Squares Mean_Square F_value P_value
Position 1 7160.0556 7160.0556 15.998138 0.001762
Temperature 2 945342.1111 472671.0556 1056.117056 0.000000
Position:Temperature 2 818.1111 409.0556 0.913977 0.427110
Residuals 12 5370.6667 447.5556 NA NA

Store the mean squares and degrees of freedom for the fixed- and random-effects calculations.

SS <- q2_anova$`Sum Sq`
DF <- q2_anova$Df

names(SS) <- rownames(q2_anova)
names(DF) <- rownames(q2_anova)

MS <- SS / DF

MS
##             Position          Temperature Position:Temperature 
##            7160.0556          472671.0556             409.0556 
##            Residuals 
##             447.5556

2.5 Part (a): Both factors fixed

When both Position and Temperature are fixed effects, the two main effects and their interaction are tested against the residual mean square.

q2_fixed <- data.frame(
  Effect = rownames(q2_anova),
  Df = q2_anova$Df,
  F_value = q2_anova$`F value`,
  P_value = q2_anova$`Pr(>F)`,
  row.names = NULL
)

knitr::kable(
  q2_fixed,
  digits = 6,
  caption = "Part (a): Both factors fixed"
)
Part (a): Both factors fixed
Effect Df F_value P_value
Position 1 15.998138 0.001762
Temperature 2 1056.117056 0.000000
Position:Temperature 2 0.913977 0.427110
Residuals 12 NA NA

The null hypotheses are:

\[ H_{0P}: \text{there is no position effect}, \]

\[ H_{0T}: \text{there is no temperature effect}, \]

\[ H_{0PT}: \text{there is no position-by-temperature interaction}. \]

Use the p-values to determine which null hypotheses are rejected at the 0.05 significance level.

2.6 Part (b): Both factors random

When both factors are random effects in a balanced crossed design with replication, the main effects are tested against the position-by-temperature interaction mean square.

The interaction is tested against the residual mean square.

F_position_random <- MS["Position"] /
  MS["Position:Temperature"]

p_position_random <- pf(
  F_position_random,
  df1 = DF["Position"],
  df2 = DF["Position:Temperature"],
  lower.tail = FALSE
)

F_temperature_random <- MS["Temperature"] /
  MS["Position:Temperature"]

p_temperature_random <- pf(
  F_temperature_random,
  df1 = DF["Temperature"],
  df2 = DF["Position:Temperature"],
  lower.tail = FALSE
)

F_interaction <- MS["Position:Temperature"] /
  MS["Residuals"]

p_interaction <- pf(
  F_interaction,
  df1 = DF["Position:Temperature"],
  df2 = DF["Residuals"],
  lower.tail = FALSE
)

q2_random <- data.frame(
  Effect = c(
    "Position",
    "Temperature",
    "Position:Temperature"
  ),

  F_value = c(
    F_position_random,
    F_temperature_random,
    F_interaction
  ),

  Numerator_df = c(
    DF["Position"],
    DF["Temperature"],
    DF["Position:Temperature"]
  ),

  Denominator_df = c(
    DF["Position:Temperature"],
    DF["Position:Temperature"],
    DF["Residuals"]
  ),

  P_value = c(
    p_position_random,
    p_temperature_random,
    p_interaction
  ),

  row.names = NULL
)

knitr::kable(
  q2_random,
  digits = 6,
  caption = "Part (b): Both factors random"
)
Part (b): Both factors random
Effect F_value Numerator_df Denominator_df P_value
Position 17.503871 1 2 0.052658
Temperature 1155.517995 2 2 0.000865
Position:Temperature 0.913977 2 12 0.427110

The significance decision for each effect is determined by comparing its p-value with 0.05.

2.7 Part (c): Position fixed and temperature random

For the mixed-effects model, Position is fixed and Temperature is random.

In this balanced crossed design, the position and temperature main effects are tested against the interaction mean square. The interaction is tested against the residual mean square.

q2_mixed <- data.frame(
  Effect = c(
    "Position (fixed)",
    "Temperature (random)",
    "Position:Temperature"
  ),

  F_value = c(
    MS["Position"] / MS["Position:Temperature"],
    MS["Temperature"] / MS["Position:Temperature"],
    MS["Position:Temperature"] / MS["Residuals"]
  ),

  Numerator_df = c(
    DF["Position"],
    DF["Temperature"],
    DF["Position:Temperature"]
  ),

  Denominator_df = c(
    DF["Position:Temperature"],
    DF["Position:Temperature"],
    DF["Residuals"]
  ),

  P_value = c(
    pf(
      MS["Position"] / MS["Position:Temperature"],
      DF["Position"],
      DF["Position:Temperature"],
      lower.tail = FALSE
    ),

    pf(
      MS["Temperature"] / MS["Position:Temperature"],
      DF["Temperature"],
      DF["Position:Temperature"],
      lower.tail = FALSE
    ),

    pf(
      MS["Position:Temperature"] / MS["Residuals"],
      DF["Position:Temperature"],
      DF["Residuals"],
      lower.tail = FALSE
    )
  ),

  row.names = NULL
)

knitr::kable(
  q2_mixed,
  digits = 6,
  caption = "Part (c): Position fixed and temperature random"
)
Part (c): Position fixed and temperature random
Effect F_value Numerator_df Denominator_df P_value
Position (fixed) 17.503871 1 2 0.052658
Temperature (random) 1155.517995 2 2 0.000865
Position:Temperature 0.913977 2 12 0.427110

2.8 Part (d): Compare the three analyses

Compare the p-values from the three model specifications.

q2_comparison <- data.frame(
  Effect = c(
    "Position",
    "Temperature",
    "Position:Temperature"
  ),

  Both_Fixed = c(
    q2_fixed$P_value[
      q2_fixed$Effect == "Position"
    ],
    q2_fixed$P_value[
      q2_fixed$Effect == "Temperature"
    ],
    q2_fixed$P_value[
      q2_fixed$Effect == "Position:Temperature"
    ]
  ),

  Both_Random = c(
    p_position_random,
    p_temperature_random,
    p_interaction
  ),

  Position_Fixed_Temperature_Random = c(
    q2_mixed$P_value[1],
    q2_mixed$P_value[2],
    q2_mixed$P_value[3]
  )
)

knitr::kable(
  q2_comparison,
  digits = 6,
  caption = "Comparison of p-values under the three models"
)
Comparison of p-values under the three models
Effect Both_Fixed Both_Random Position_Fixed_Temperature_Random
Position Position 0.001762 0.052658 0.052658
Temperature Temperature 0.000000 0.000865 0.000865
Position:Temperature Position:Temperature 0.427110 0.427110 0.427110

The fixed-effects model tests the main effects against the residual mean square. The random-effects and mixed-effects models use the interaction mean square as the denominator for the main-effect tests in this balanced crossed design.

Changing the denominator changes the F-statistics and can change the p-values and significance decisions.

Parts (b) and (c) use the same denominators for the corresponding tests, so their p-values are equal for these data. The interaction also uses the residual mean square in all three analyses, giving the same interaction test.

The final conclusions should be based on the p-values in the comparison table.

2.9 Interaction plot

q2_means <- aggregate(
  Density ~ Position + Temperature,
  data = anode,
  FUN = mean
)

ggplot(
  q2_means,
  aes(
    x = Temperature,
    y = Density,
    color = Position,
    group = Position
  )
) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 3) +
  labs(
    title = "Carbon-anode density by temperature and position",
    x = "Firing temperature (degrees C)",
    y = "Mean baked density",
    color = "Furnace position"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Nonparallel lines suggest an interaction between furnace position and firing temperature. The interaction p-value from the ANOVA determines whether there is statistically significant evidence of that interaction.