============================================================
library(tidyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(GAD)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(agricolae)
# ============================================================
# PROBLEM 3.7 (c,d,e,f) — Portland Cement
# ============================================================
cement <- data.frame(
Technique = factor(rep(1:4, each = 4)),
Strength = c(3129, 3000, 2865, 2890,
3200, 3300, 2975, 3150,
2800, 2900, 2985, 3050,
2600, 2700, 2600, 2765)
)
# Fit model
cement_aov <- aov(Strength ~ Technique, data = cement)
summary(cement_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Technique 3 489740 163247 12.73 0.000489 ***
## Residuals 12 153908 12826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (c) Fisher's LSD
LSD_cement <- LSD.test(cement_aov, "Technique", group = FALSE, console = TRUE)
##
## Study: cement_aov ~ "Technique"
##
## LSD t Test for Strength
##
## Mean Square Error: 12825.69
##
## Technique, means and individual ( 95 %) CI
##
## Strength std r se LCL UCL Min Max Q25 Q50
## 1 2971.00 120.55704 4 56.62528 2847.624 3094.376 2865 3129 2883.75 2945.0
## 2 3156.25 135.97641 4 56.62528 3032.874 3279.626 2975 3300 3106.25 3175.0
## 3 2933.75 108.27242 4 56.62528 2810.374 3057.126 2800 3050 2875.00 2942.5
## 4 2666.25 80.97067 4 56.62528 2542.874 2789.626 2600 2765 2600.00 2650.0
## Q75
## 1 3032.25
## 2 3225.00
## 3 3001.25
## 4 2716.25
##
## Alpha: 0.05 ; DF Error: 12
## Critical Value of t: 2.178813
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## 1 - 2 -185.25 0.0392 * -359.72984 -10.77016
## 1 - 3 37.25 0.6501 -137.22984 211.72984
## 1 - 4 304.75 0.0025 ** 130.27016 479.22984
## 2 - 3 222.50 0.0167 * 48.02016 396.97984
## 2 - 4 490.00 0.0001 *** 315.52016 664.47984
## 3 - 4 267.50 0.0059 ** 93.02016 441.97984
# (d) Normal probability plot of residuals
par(mfrow = c(2, 2))
plot(cement_aov)
par(mfrow = c(1, 1))
shapiro.test(residuals(cement_aov))
##
## Shapiro-Wilk normality test
##
## data: residuals(cement_aov)
## W = 0.97046, p-value = 0.846
# (e) Residuals vs predicted
plot(fitted(cement_aov), residuals(cement_aov),
xlab = "Predicted", ylab = "Residuals",
main = "Residuals vs Predicted: Cement")
abline(h = 0, col = "red", lty = 2)
# (f) Scatter plot of results
plot(cement$Technique, cement$Strength,
xlab = "Mixing Technique", ylab = "Tensile Strength (lb/in^2)",
main = "Scatter Plot of Cement Tensile Strength",
col = "steelblue", pch = 19)
# ============================================================
# PROBLEM 3.10 (b,c) — Cotton Content
# ============================================================
cotton <- data.frame(
Cotton = factor(rep(c(15, 20, 25, 30, 35), each = 5)),
Strength = c(7,7,15,11,9,
12,17,12,18,18,
14,19,19,18,18,
19,25,22,19,23,
7,10,11,15,11)
)
cotton_aov <- aov(Strength ~ Cotton, data = cotton)
summary(cotton_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Cotton 4 475.8 118.94 14.76 9.13e-06 ***
## Residuals 20 161.2 8.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (b) Fisher's LSD
LSD_cotton <- LSD.test(cotton_aov, "Cotton", group = FALSE, console = TRUE)
##
## Study: cotton_aov ~ "Cotton"
##
## LSD t Test for Strength
##
## Mean Square Error: 8.06
##
## Cotton, means and individual ( 95 %) CI
##
## Strength std r se LCL UCL Min Max Q25 Q50 Q75
## 15 9.8 3.346640 5 1.269646 7.151566 12.44843 7 15 7 9 11
## 20 15.4 3.130495 5 1.269646 12.751566 18.04843 12 18 12 17 18
## 25 17.6 2.073644 5 1.269646 14.951566 20.24843 14 19 18 18 19
## 30 21.6 2.607681 5 1.269646 18.951566 24.24843 19 25 19 22 23
## 35 10.8 2.863564 5 1.269646 8.151566 13.44843 7 15 10 11 11
##
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 2.085963
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## 15 - 20 -5.6 0.0054 ** -9.3454518 -1.8545482
## 15 - 25 -7.8 0.0003 *** -11.5454518 -4.0545482
## 15 - 30 -11.8 0.0000 *** -15.5454518 -8.0545482
## 15 - 35 -1.0 0.5838 -4.7454518 2.7454518
## 20 - 25 -2.2 0.2347 -5.9454518 1.5454518
## 20 - 30 -6.2 0.0025 ** -9.9454518 -2.4545482
## 20 - 35 4.6 0.0186 * 0.8545482 8.3454518
## 25 - 30 -4.0 0.0375 * -7.7454518 -0.2545482
## 25 - 35 6.8 0.0012 ** 3.0545482 10.5454518
## 30 - 35 10.8 0.0000 *** 7.0545482 14.5454518
# (c) Residual analysis
par(mfrow = c(2, 2))
plot(cotton_aov)
par(mfrow = c(1, 1))
shapiro.test(residuals(cotton_aov))
##
## Shapiro-Wilk normality test
##
## data: residuals(cotton_aov)
## W = 0.94387, p-value = 0.1818
# ============================================================
# PROBLEM 3.20 (a,b) — Rodding of Concrete
# ============================================================
rodding <- data.frame(
Rodding = factor(rep(c(10, 15, 20, 25), each = 3)),
Strength = c(1530,1530,1440,
1610,1650,1500,
1560,1730,1530,
1500,1490,1510)
)
rodding_aov <- aov(Strength ~ Rodding, data = rodding)
summary(rodding_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Rodding 3 28633 9544 1.865 0.214
## Residuals 8 40933 5117
# ============================================================
# PROBLEM 3.23 — Insulating Fluid Life
# ============================================================
fluid <- data.frame(
Fluid = factor(rep(1:4, each = 6)),
Life = c(17.6,18.9,16.3,17.4,20.1,21.6,
16.9,15.3,18.6,17.1,19.5,20.3,
21.4,23.6,19.4,18.5,20.5,22.3,
19.3,21.1,16.9,17.5,18.3,19.8)
)
fluid_aov <- aov(Life ~ Fluid, data = fluid)
summary(fluid_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fluid 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
# (b) Which fluid to select?
fluid %>% group_by(Fluid) %>% summarise(Mean = mean(Life))
## # A tibble: 4 × 2
## Fluid Mean
## <fct> <dbl>
## 1 1 18.6
## 2 2 18.0
## 3 3 21.0
## 4 4 18.8
# (c) Residuals
par(mfrow = c(2,2)); plot(fluid_aov); par(mfrow = c(1,1))
shapiro.test(residuals(fluid_aov))
##
## Shapiro-Wilk normality test
##
## data: residuals(fluid_aov)
## W = 0.95671, p-value = 0.376
# ============================================================
# PROBLEM 3.28 — Insulating Materials Failure Times
# ============================================================
materials <- data.frame(
Material = factor(rep(1:5, each = 4)),
Failure = c(110,157,194,178,
1,2,4,18,
880,1256,5276,4355,
495,7040,5307,10050,
7,5,29,2)
)
materials_aov <- aov(Failure ~ Material, data = materials)
summary(materials_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Material 4 103191489 25797872 6.191 0.00379 **
## Residuals 15 62505657 4167044
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (b) Residual plots
par(mfrow = c(2,2)); plot(materials_aov); par(mfrow = c(1,1))
# (c) Log transform (since variance is highly unequal)
materials$logFailure <- log(materials$Failure)
materials_log_aov <- aov(logFailure ~ Material, data = materials)
summary(materials_log_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Material 4 165.06 41.26 37.66 1.18e-07 ***
## Residuals 15 16.44 1.10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ============================================================
# PROBLEM 3.29 — Semiconductor Particle Counts
# ============================================================
particle <- data.frame(
Method = factor(rep(1:3, each = 5)),
Count = c(31,10,21,4,1,
62,40,24,30,35,
53,27,120,97,68)
)
particle_aov <- aov(Count ~ Method, data = particle)
summary(particle_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Method 2 8964 4482 7.914 0.00643 **
## Residuals 12 6796 566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (b) Residuals
par(mfrow = c(2,2)); plot(particle_aov); par(mfrow = c(1,1))
# (c) Kruskal-Wallis
kruskal.test(Count ~ Method, data = particle)
##
## Kruskal-Wallis rank sum test
##
## data: Count by Method
## Kruskal-Wallis chi-squared = 8.54, df = 2, p-value = 0.01398
# ============================================================
# PROBLEM 3.44 & 3.45 — Sample Size for Four Means
# ============================================================
# Means: 50, 60, 50, 60; sigma^2 = 25
# Use power.anova.test
power.anova.test(
groups = 4,
n = NULL,
between.var = var(c(50, 60, 50, 60)),
within.var = 25,
sig.level = 0.05,
power = 0.90
)
##
## Balanced one-way analysis of variance power calculation
##
## groups = 4
## n = 4.658128
## between.var = 33.33333
## within.var = 25
## sig.level = 0.05
## power = 0.9
##
## NOTE: n is number in each group
# 3.45(a) sigma^2 = 36
power.anova.test(groups = 4, n = NULL,
between.var = var(c(50,60,50,60)),
within.var = 36, sig.level = 0.05, power = 0.90)
##
## Balanced one-way analysis of variance power calculation
##
## groups = 4
## n = 6.180885
## between.var = 33.33333
## within.var = 36
## sig.level = 0.05
## power = 0.9
##
## NOTE: n is number in each group
# 3.45(b) sigma^2 = 49
power.anova.test(groups = 4, n = NULL,
between.var = var(c(50,60,50,60)),
within.var = 49, sig.level = 0.05, power = 0.90)
##
## Balanced one-way analysis of variance power calculation
##
## groups = 4
## n = 7.998751
## between.var = 33.33333
## within.var = 49
## sig.level = 0.05
## power = 0.9
##
## NOTE: n is number in each group
# ============================================================
# PROBLEM 3.51 & 3.52 — Kruskal-Wallis on Problem 3.23
# ============================================================
kruskal.test(Life ~ Fluid, data = fluid)
##
## Kruskal-Wallis rank sum test
##
## data: Life by Fluid
## Kruskal-Wallis chi-squared = 6.2177, df = 3, p-value = 0.1015
# ============================================================
# PROBLEM 4.3 — RCBD: Chemical on Cloth (fixed blocks)
# ============================================================
cloth <- data.frame(
Chemical = factor(rep(1:4, each = 5)),
Bolt = factor(rep(1:5, times = 4)),
Strength = c(73,68,74,71,67,
73,67,75,72,70,
75,68,78,73,68,
73,71,75,75,69)
)
cloth$Chemical_GAD <- as.fixed(cloth$Chemical)
cloth$Bolt_GAD <- as.fixed(cloth$Bolt) # 4.16 asks for fixed bolts
cloth_model <- lm(Strength ~ Chemical_GAD + Bolt_GAD, data = cloth)
gad(cloth_model)
## $anova
## Analysis of Variance Table
##
## Response: Strength
## Df Sum Sq Mean Sq F value Pr(>F)
## Chemical_GAD 3 12.95 4.317 2.3761 0.1211
## Bolt_GAD 4 157.00 39.250 21.6055 2.059e-05 ***
## Residuals 12 21.80 1.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ============================================================
# PROBLEM 4.16 — RCBD: Estimate model parameters
# ============================================================
# With fixed chemical and fixed bolts:
cloth_model_fixed <- aov(Strength ~ Chemical + Bolt, data = cloth)
summary(cloth_model_fixed)
## Df Sum Sq Mean Sq F value Pr(>F)
## Chemical 3 12.95 4.32 2.376 0.121
## Bolt 4 157.00 39.25 21.606 2.06e-05 ***
## Residuals 12 21.80 1.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.tables(cloth_model_fixed, type = "effects")
## Tables of effects
##
## Chemical
## Chemical
## 1 2 3 4
## -1.15 -0.35 0.65 0.85
##
## Bolt
## Bolt
## 1 2 3 4 5
## 1.75 -3.25 3.75 1.00 -3.25
# ============================================================
# PROBLEM 4.22 — Latin Square Design
# ============================================================
# Treatments: A, B, C, D, E
# Rows: Batch (1-5); Columns: Day (1-5)
lat <- data.frame(
Batch = factor(rep(1:5, each = 5)),
Day = factor(rep(1:5, times = 5)),
Treat = factor(c("A","B","D","C","E",
"C","E","A","D","B",
"B","A","C","E","D",
"D","C","E","B","A",
"E","D","B","A","C")),
Resp = c(8,7,1,7,3,
11,2,7,3,8,
4,9,10,1,5,
6,8,6,6,10,
4,2,3,8,8)
)
lat_model <- aov(Resp ~ Batch + Day + Treat, data = lat)
summary(lat_model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Batch 4 15.44 3.86 1.235 0.347618
## Day 4 12.24 3.06 0.979 0.455014
## Treat 4 141.44 35.36 11.309 0.000488 ***
## Residuals 12 37.52 3.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Residuals for Latin square
par(mfrow = c(2,2)); plot(lat_model); par(mfrow = c(1,1))
