# Memuat data
espresso <- read.csv("D:/Kuliah/Semester 3/Rancangan Percobaan/Proyek G/espresso2.csv")
espresso
## foamIndx trt_id tempC prssBar
## 1 94.21 1 75 15
## 2 120.39 1 75 15
## 3 119.72 1 75 15
## 4 120.45 1 75 15
## 5 132.68 1 75 15
## 6 121.32 1 75 15
## 7 88.85 1 75 15
## 8 100.33 1 75 15
## 9 122.65 1 75 15
## 10 134.84 2 75 20
## 11 154.17 2 75 20
## 12 138.15 2 75 20
## 13 109.89 2 75 20
## 14 178.94 2 75 20
## 15 110.92 2 75 20
## 16 169.65 2 75 20
## 17 110.27 2 75 20
## 18 108.18 2 75 20
## 19 95.84 3 85 15
## 20 98.26 3 85 15
## 21 102.74 3 85 15
## 22 102.00 3 85 15
## 23 94.37 3 85 15
## 24 88.05 3 85 15
## 25 97.54 3 85 15
## 26 128.96 3 85 15
## 27 113.84 3 85 15
## 28 143.38 4 85 20
## 29 137.78 4 85 20
## 30 85.06 4 85 20
## 31 112.60 4 85 20
## 32 99.79 4 85 20
## 33 102.44 4 85 20
## 34 116.04 4 85 20
## 35 94.14 4 85 20
## 36 132.97 4 85 20
## 37 75.63 5 95 15
## 38 78.38 5 95 15
## 39 113.45 5 95 15
## 40 104.62 5 95 15
## 41 92.47 5 95 15
## 42 96.92 5 95 15
## 43 80.92 5 95 15
## 44 95.74 5 95 15
## 45 81.78 5 95 15
## 46 120.10 6 95 20
## 47 136.54 6 95 20
## 48 118.50 6 95 20
## 49 95.86 6 95 20
## 50 144.83 6 95 20
## 51 79.89 6 95 20
## 52 102.04 6 95 20
## 53 140.31 6 95 20
## 54 120.33 6 95 20
# Memastikan kolom 'tempC' dianggap sebagai faktor
espresso$tempC <- as.factor(espresso$tempC)
espresso$trt_id <- as.factor(espresso$trt_id)
# Melakukan analisis variansi
model <- aov(foamIndx ~ trt_id + tempC, data = espresso)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## trt_id 5 9848 1969.7 5.402 0.000509 ***
## Residuals 48 17501 364.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji lanjut Tukey untuk perlakuan
TukeyHSD(model, which = "trt_id")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = foamIndx ~ trt_id + tempC, data = espresso)
##
## $trt_id
## diff lwr upr p adj
## 2-1 21.60111 -5.114020 48.316242 0.1768428
## 3-1 -11.00000 -37.715131 15.715131 0.8239626
## 4-1 0.40000 -26.315131 27.115131 1.0000000
## 5-1 -22.29889 -49.014020 4.416242 0.1513853
## 6-1 4.20000 -22.515131 30.915131 0.9970720
## 3-2 -32.60111 -59.316242 -5.885980 0.0086879
## 4-2 -21.20111 -47.916242 5.514020 0.1928126
## 5-2 -43.90000 -70.615131 -17.184869 0.0001705
## 6-2 -17.40111 -44.116242 9.314020 0.3952510
## 4-3 11.40000 -15.315131 38.115131 0.8013930
## 5-3 -11.29889 -38.014020 15.416242 0.8072165
## 6-3 15.20000 -11.515131 41.915131 0.5457826
## 5-4 -22.69889 -49.414020 4.016242 0.1381232
## 6-4 3.80000 -22.915131 30.515131 0.9981822
## 6-5 26.49889 -0.216242 53.214020 0.0530291
# Pengecekan asumsi
par(mfrow = c(2, 2))
plot(model)
# Uji normalitas residual
shapiro.test(residuals(model))
##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.98577, p-value = 0.7674
# Uji kesamaan variansi antar perlakuan
bartlett.test(foamIndx ~ trt_id, data = espresso)
##
## Bartlett test of homogeneity of variances
##
## data: foamIndx by trt_id
## Bartlett's K-squared = 7.9006, df = 5, p-value = 0.1618
# Uji kesamaan variansi antar kelompok
bartlett.test(foamIndx ~ tempC, data = espresso)
##
## Bartlett test of homogeneity of variances
##
## data: foamIndx by tempC
## Bartlett's K-squared = 1.8234, df = 2, p-value = 0.4018
#1. Boxplot Respon Berdasarkan Perlakuan
library(ggplot2)
ggplot(espresso, aes(x = trt_id, y = foamIndx, fill = trt_id)) +
geom_boxplot() +
theme_minimal() +
labs(
title = "Boxplot of Foam Index by Treatment",
x = "Treatment",
y = "Foam Index"
) +
scale_fill_brewer(palette = "Set3")
#2. Boxplot Respon Berdasarkan Kelompok
ggplot(espresso, aes(x = tempC, y = foamIndx, fill = tempC)) +
geom_boxplot() +
theme_minimal() +
labs(
title = "Boxplot of Foam Index by Temperature Group",
x = "Temperature (°C)",
y = "Foam Index"
) +
scale_fill_brewer(palette = "Set2")
#3. Scatter Plot untuk Hubungan Suhu dan Indeks Busa
ggplot(espresso, aes(x = as.numeric(tempC), y = foamIndx, color = trt_id)) +
geom_point(size = 3) +
theme_minimal() +
labs(
title = "Scatter Plot of Foam Index by Temperature",
x = "Temperature (°C)",
y = "Foam Index",
color = "Treatment"
) +
scale_color_brewer(palette = "Dark2")
#histogram respon
ggplot(espresso, aes(x = foamIndx)) +
geom_histogram(binwidth = 10, fill = "steelblue", color = "black") +
theme_minimal() +
labs(
title = "Histogram of Foam Index",
x = "Foam Index",
y = "Frequency"
)
#interaksi plot untuk perlakuan dan kelompok
interaction.plot(
x.factor = espresso$tempC,
trace.factor = espresso$trt_id,
response = espresso$foamIndx,
col = rainbow(length(levels(espresso$trt_id))),
type = "b",
legend = TRUE,
xlab = "Temperature (°C)",
ylab = "Foam Index",
main = "Interaction Plot"
)
# Asumsi Normalitas
par(mfrow = c(1, 1)) # Reset layout
plot(model, 2) # QQ plot untuk residual
shapiro.test(residuals(model)) # Uji Shapiro-Wilk
##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.98577, p-value = 0.7674
# Asumsi Kesamaan Variansi
# Variansi antar perlakuan
bartlett.test(foamIndx ~ trt_id, data = espresso)
##
## Bartlett test of homogeneity of variances
##
## data: foamIndx by trt_id
## Bartlett's K-squared = 7.9006, df = 5, p-value = 0.1618
# Variansi antar kelompok
bartlett.test(foamIndx ~ tempC, data = espresso)
##
## Bartlett test of homogeneity of variances
##
## data: foamIndx by tempC
## Bartlett's K-squared = 1.8234, df = 2, p-value = 0.4018
# Visualisasi Residuals
library(ggplot2)
# Residuals vs Perlakuan
p1 <- ggplot(espresso, aes(x = trt_id, y = residuals(model))) +
geom_point(size = 3, position = "jitter") +
labs(x = "Jenis Perlakuan (trt_id)", y = "Residual") +
theme_bw()
p1
# Residuals vs Kelompok
p2 <- ggplot(espresso, aes(x = tempC, y = residuals(model))) +
geom_point(size = 3, position = "jitter") +
labs(x = "Kelompok Suhu (tempC)", y = "Residual") +
theme_bw()
p2
# Gabung plot residual
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.4.1
ggarrange(p1, p2, labels = c("A", "B"), ncol = 2, nrow = 1)
# Asumsi Aditivitas
interaction.plot(
x.factor = espresso$tempC,
trace.factor = espresso$trt_id,
response = espresso$foamIndx,
col = 1:length(unique(espresso$trt_id)),
type = "b",
legend = TRUE,
xlab = "Kelompok Suhu (tempC)",
ylab = "Rata-rata Foam Index",
main = "Interaction Plot"
)
