“This analysis explores data from a reliability engineering study designed to select the optimal power source for industrial sensors used in harsh environments (e.g., downhole drilling or autoclaves). The experiment investigates the Service Life (in hours) of the batteries based on two factors:
Here is the modified and formatted text, ready for your R Markdown file:
1: Standard Li-ion (Consumer Grade)2: Industrial Li-MnO\(_2\) (Lithium Manganese Dioxide)3: High-Temp Li-SOCl\(_2\) (Lithium Thionyl Chloride)The objective is to characterize the reliability of these power sources across the full thermal range. A critical specific aim is to investigate the presence of an interaction effect: to determine if the ‘High-Temp’ chemistry successfully maintains its voltage stability at 125°C, a point where standard consumer-grade batteries are expected to suffer catastrophic capacity loss.
First, we organize the raw data into a structured data.frame. This is a best practice in R that makes the data easier to manage and the code more readable. We create columns for the response variable life and the two factors, material and temperature, ensuring they are treated as categorical variables (factors) for the analysis.
## Response variable: battery life
life <- c(130,155,74,180, 34,40,80,75, 20,70,82,58,
150,188,159,126, 136,122,106,115, 25,70,58,45,
138,110,168,160, 174,120,150,139, 96,104,82,60)
## Create the data frame
battery_df <- data.frame(
life = life,
material = factor(rep(1:3, each = 12)),
temperature = factor(rep(rep(c(15, 70, 125), each = 4), 3))
)
## Preview the data
battery_dfBefore fitting a formal model, we visualize the data to get an intuition for the relationships between the factors and the response.
Boxplots are excellent for examining the distribution of battery life for each level of our factors independently. This gives us a preliminary look at the main effects—the individual impact of material type and temperature.
library(ggplot2)
## Boxplot for Material Type
ggplot(battery_df, aes(x = material, y = life, fill = material)) +
geom_boxplot() +
labs(title = "Battery Life by Material Type", x = "Material Type", y = "Life (hours)") +
theme_minimal() +
theme(legend.position = "none")
## Boxplot for Temperature
ggplot(battery_df, aes(x = temperature, y = life, fill = temperature)) +
geom_boxplot() +
labs(title = "Battery Life by Temperature", x = "Temperature (°C)", y = "Life (hours)") +
theme_minimal() +
theme(legend.position = "none")
The most crucial plot for a factorial experiment is the interaction plot. It displays the mean battery life for each combination of material and temperature. If the lines are parallel, it suggests there is no interaction. If the lines are not parallel (i.e., they cross or diverge), it indicates that the effect of temperature on battery life is different for each material type, signaling a likely interaction.
ggplot(battery_df, aes(x = temperature, y = life, group = material, color = material)) +
stat_summary(fun = mean, geom = "line", size = 1) +
stat_summary(fun = mean, geom = "point", size = 3) +
labs(
title = "Interaction Plot: Material Type and Temperature",
x = "Temperature (°C)",
y = "Average Battery Life (hours)",
color = "Material Type"
) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
The interaction plot reveals a strong interaction effect between Material Type and Temperature, meaning the performance profile is distinct for each chemistry:
Summary: There is no single “best” material overall. The optimal choice depends entirely on the environment: Material 2 for cold applications and Material 3 for any environment exceeding room temperature.
We now fit a linear model to formally test the significance of the main effects and the interaction term. The model life ~ material * temperature is shorthand for life ~ material + temperature + material:temperature. We use a sum-to-zero contrast (contr.sum) for balanced interpretation of the effects. The ANOVA table will tell us if the variation caused by our factors is statistically significant compared to the random variation in the data.
## Fit the full factorial model
battery_fit <- lm(life ~ material * temperature,
data = battery_df,
contrasts = list(material = contr.sum, temperature = contr.sum))
summary(battery_fit)
Call:
lm(formula = life ~ material * temperature, data = battery_df,
contrasts = list(material = contr.sum, temperature = contr.sum))
Residuals:
Min 1Q Median 3Q Max
-60.750 -14.625 1.375 17.938 45.250
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 105.528 4.331 24.367 < 2e-16 ***
material1 -22.361 6.125 -3.651 0.00111 **
material2 2.806 6.125 0.458 0.65057
temperature1 39.306 6.125 6.418 7.1e-07 ***
temperature2 2.056 6.125 0.336 0.73975
material1:temperature1 12.278 8.662 1.417 0.16778
material2:temperature1 8.111 8.662 0.936 0.35735
material1:temperature2 -27.972 8.662 -3.229 0.00325 **
material2:temperature2 9.361 8.662 1.081 0.28936
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25.98 on 27 degrees of freedom
Multiple R-squared: 0.7652, Adjusted R-squared: 0.6956
F-statistic: 11 on 8 and 27 DF, p-value: 9.426e-07ANOVA Results
## Generate the ANOVA table
knitr::kable(anova(battery_fit), digits=3)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| material | 2 | 10683.722 | 5341.861 | 7.911 | 0.002 |
| temperature | 2 | 39118.722 | 19559.361 | 28.968 | 0.000 |
| material:temperature | 4 | 9613.778 | 2403.444 | 3.560 | 0.019 |
| Residuals | 27 | 18230.750 | 675.213 | NA | NA |
The ANOVA table shows very small p-values (Pr(>F)) for material, temperature, and, most importantly, the material:temperature interaction. This confirms our visual inspection: all effects are statistically significant. Because the interaction is significant, our interpretation should focus on the interaction itself rather than the main effects in isolation.
The validity of our ANOVA results depends on the model’s residuals meeting certain assumptions (normality, constant variance, independence). We check these with diagnostic plots.
## Extract standardized residuals and fitted values
battery_fit_diag <- data.frame(
residuals = rstandard(battery_fit),
fitted = fitted.values(battery_fit)
)
## Normal Q-Q Plot
p1 <- ggplot(battery_fit_diag, aes(sample = residuals)) +
stat_qq() +
stat_qq_line() +
labs(title = "Normal Q-Q Plot", x = "Theoretical Quantiles", y = "Standardized Residuals") +
theme_minimal()
## Residuals vs. Fitted Plot
p2 <- ggplot(battery_fit_diag, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Residuals vs. Fitted Values", x = "Fitted Values", y = "Standardized Residuals") +
theme_minimal()
p1
p2
The Normal Q-Q plot shows the points falling roughly along the line, suggesting the normality assumption is met. The Residuals vs. Fitted plot shows a random scatter of points around the zero line, indicating that the variance is reasonably constant. The model assumptions appear to be satisfied.
Since the interaction is significant, we must compare the means of the nine specific treatment combinations (3 materials × 3 temperatures). Simply comparing the average effect of Material 1 vs. Material 2 would be misleading, as that difference depends on the temperature.
Tukey’s Honest Significant Difference (HSD) test is a post-hoc test that compares all possible pairs of means while controlling the family-wise error rate. We apply it to an aov model object. The output for the material:temperature interaction shows which specific combinations are significantly different from one another.
## Fit the model using aov() for Tukey's test
battery_aov <- aov(life ~ material * temperature, data = battery_df)
## Perform Tukey's HSD test
TukeyHSD(battery_aov) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = life ~ material * temperature, data = battery_df)
$material
diff lwr upr p adj
2-1 25.16667 -1.135677 51.46901 0.0627571
3-1 41.91667 15.614323 68.21901 0.0014162
3-2 16.75000 -9.552344 43.05234 0.2717815
$temperature
diff lwr upr p adj
70-15 -37.25000 -63.55234 -10.94766 0.0043788
125-15 -80.66667 -106.96901 -54.36432 0.0000001
125-70 -43.41667 -69.71901 -17.11432 0.0009787
$`material:temperature`
diff lwr upr p adj
2:15-1:15 21.00 -40.823184 82.823184 0.9616404
3:15-1:15 9.25 -52.573184 71.073184 0.9998527
1:70-1:15 -77.50 -139.323184 -15.676816 0.0065212
2:70-1:15 -15.00 -76.823184 46.823184 0.9953182
3:70-1:15 11.00 -50.823184 72.823184 0.9994703
1:125-1:15 -77.25 -139.073184 -15.426816 0.0067471
2:125-1:15 -85.25 -147.073184 -23.426816 0.0022351
3:125-1:15 -49.25 -111.073184 12.573184 0.2016535
3:15-2:15 -11.75 -73.573184 50.073184 0.9991463
1:70-2:15 -98.50 -160.323184 -36.676816 0.0003449
2:70-2:15 -36.00 -97.823184 25.823184 0.5819453
3:70-2:15 -10.00 -71.823184 51.823184 0.9997369
1:125-2:15 -98.25 -160.073184 -36.426816 0.0003574
2:125-2:15 -106.25 -168.073184 -44.426816 0.0001152
3:125-2:15 -70.25 -132.073184 -8.426816 0.0172076
1:70-3:15 -86.75 -148.573184 -24.926816 0.0018119
2:70-3:15 -24.25 -86.073184 37.573184 0.9165175
3:70-3:15 1.75 -60.073184 63.573184 1.0000000
1:125-3:15 -86.50 -148.323184 -24.676816 0.0018765
2:125-3:15 -94.50 -156.323184 -32.676816 0.0006078
3:125-3:15 -58.50 -120.323184 3.323184 0.0742711
2:70-1:70 62.50 0.676816 124.323184 0.0460388
3:70-1:70 88.50 26.676816 150.323184 0.0014173
1:125-1:70 0.25 -61.573184 62.073184 1.0000000
2:125-1:70 -7.75 -69.573184 54.073184 0.9999614
3:125-1:70 28.25 -33.573184 90.073184 0.8281938
3:70-2:70 26.00 -35.823184 87.823184 0.8822881
1:125-2:70 -62.25 -124.073184 -0.426816 0.0474675
2:125-2:70 -70.25 -132.073184 -8.426816 0.0172076
3:125-2:70 -34.25 -96.073184 27.573184 0.6420441
1:125-3:70 -88.25 -150.073184 -26.426816 0.0014679
2:125-3:70 -96.25 -158.073184 -34.426816 0.0004744
3:125-3:70 -60.25 -122.073184 1.573184 0.0604247
2:125-1:125 -8.00 -69.823184 53.823184 0.9999508
3:125-1:125 28.00 -33.823184 89.823184 0.8347331
3:125-2:125 36.00 -25.823184 97.823184 0.5819453The Fisher’s Least Significant Difference (LSD) method is another option for pairwise comparisons. To test the interaction means, we must specify both factors in the trt argument.
library(agricolae)
## Perform LSD test on the interaction term
lsd_results <- LSD.test(battery_aov, trt = c("material", "temperature"),
p.adj = "none", group = FALSE)
## Print the comparison table
print(lsd_results$comparison) difference pvalue signif. LCL UCL
1:125 - 1:15 -77.25 0.0003 *** -114.950479 -39.549521
1:125 - 1:70 0.25 0.9892 -37.450479 37.950479
1:125 - 2:125 8.00 0.6667 -29.700479 45.700479
1:125 - 2:15 -98.25 0.0000 *** -135.950479 -60.549521
1:125 - 2:70 -62.25 0.0022 ** -99.950479 -24.549521
1:125 - 3:125 -28.00 0.1392 -65.700479 9.700479
1:125 - 3:15 -86.50 0.0001 *** -124.200479 -48.799521
1:125 - 3:70 -88.25 0.0001 *** -125.950479 -50.549521
1:15 - 1:70 77.50 0.0002 *** 39.799521 115.200479
1:15 - 2:125 85.25 0.0001 *** 47.549521 122.950479
1:15 - 2:15 -21.00 0.2631 -58.700479 16.700479
1:15 - 2:70 15.00 0.4214 -22.700479 52.700479
1:15 - 3:125 49.25 0.0124 * 11.549521 86.950479
1:15 - 3:15 -9.25 0.6187 -46.950479 28.450479
1:15 - 3:70 -11.00 0.5544 -48.700479 26.700479
1:70 - 2:125 7.75 0.6765 -29.950479 45.450479
1:70 - 2:15 -98.50 0.0000 *** -136.200479 -60.799521
1:70 - 2:70 -62.50 0.0021 ** -100.200479 -24.799521
1:70 - 3:125 -28.25 0.1358 -65.950479 9.450479
1:70 - 3:15 -86.75 0.0001 *** -124.450479 -49.049521
1:70 - 3:70 -88.50 0.0000 *** -126.200479 -50.799521
2:125 - 2:15 -106.25 0.0000 *** -143.950479 -68.549521
2:125 - 2:70 -70.25 0.0007 *** -107.950479 -32.549521
2:125 - 3:125 -36.00 0.0605 . -73.700479 1.700479
2:125 - 3:15 -94.50 0.0000 *** -132.200479 -56.799521
2:125 - 3:70 -96.25 0.0000 *** -133.950479 -58.549521
2:15 - 2:70 36.00 0.0605 . -1.700479 73.700479
2:15 - 3:125 70.25 0.0007 *** 32.549521 107.950479
2:15 - 3:15 11.75 0.5279 -25.950479 49.450479
2:15 - 3:70 10.00 0.5907 -27.700479 47.700479
2:70 - 3:125 34.25 0.0732 . -3.450479 71.950479
2:70 - 3:15 -24.25 0.1980 -61.950479 13.450479
2:70 - 3:70 -26.00 0.1685 -63.700479 11.700479
3:125 - 3:15 -58.50 0.0036 ** -96.200479 -20.799521
3:125 - 3:70 -60.25 0.0029 ** -97.950479 -22.549521
3:15 - 3:70 -1.75 0.9248 -39.450479 35.950479The results from both Tukey’s HSD and Fisher’s LSD provide detailed p-values for comparing pairs of treatment combinations, allowing us to make specific conclusions, such as “at 125°C, Material 3 has a significantly longer life than Materials 1 and 2.”
ANOEN Table
library(ANOEN)
print(anoen(battery_fit), format = "gt")| Analysis of Entropy (ANOEN) | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Model Family: gaussian | ||||||||||
| Term |
Model Complexity
|
Entropy Info
|
Entropy-based \(R^2\)
|
Significance
|
||||||
| \(p\) | \(\Delta p\) | \(\text{df}_c\) | \(d_{\text{adj}}\) | \(\hat I_H\) | Partial | Comp. | Cum. | \(\chi^2\) | \(F_{Seq}\) | |
| material | 3.0000 | 2.0000 | 2.1183 | 10.4533 | 0.0892 | 0.0853 | 0.0853 | 0.0853 | 0.0775 | 0.0869 |
| temperature | 5.0000 | 2.0000 | 2.2507 | 9.6383 | 0.8150 | 0.5574 | 0.5098 | 0.5951 | 0.0000 | 0.0000 |
| material:temperature | 9.0000 | 4.0000 | 4.9734 | 9.3529 | 0.2854 | 0.2483 | 0.1005 | 0.6956 | 0.0092 | 0.0186 |
| Model | 9.0000 | 8.0000 | 9.3424 | 9.3529 | 1.1896 | 0.6956 | 0.6956 | 0.6956 | 0.0000 | 0.0000 |
ANOEN Pie Chart
Figure 8.1 shows the component \(R^2_\text{adj}\). This pie chart shows us the percentage of contributions of each term in the total variance of battery life.
graph.compR2(anoen(battery_fit))