library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
wine = read_excel("C:/Users/user/Downloads/Wine_Quality.xlsx")
head(wine)
## # A tibble: 6 × 7
## Clarity Aroma Body Flavor Oakiness Quality Region
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 3.3 2.8 3.1 4.1 9.8 1
## 2 1 4.4 4.9 3.5 3.9 12.6 1
## 3 1 3.9 5.3 4.8 4.7 11.9 1
## 4 1 3.9 2.6 3.1 3.6 11.1 1
## 5 1 5.6 5.1 5.5 5.1 13.3 1
## 6 1 4.6 4.7 5 4.1 12.8 1
str(wine)
## tibble [38 × 7] (S3: tbl_df/tbl/data.frame)
## $ Clarity : num [1:38] 1 1 1 1 1 1 1 1 1 1 ...
## $ Aroma : num [1:38] 3.3 4.4 3.9 3.9 5.6 4.6 4.8 5.3 4.3 4.3 ...
## $ Body : num [1:38] 2.8 4.9 5.3 2.6 5.1 4.7 4.8 4.5 4.3 3.9 ...
## $ Flavor : num [1:38] 3.1 3.5 4.8 3.1 5.5 5 4.8 4.3 3.9 4.7 ...
## $ Oakiness: num [1:38] 4.1 3.9 4.7 3.6 5.1 4.1 3.3 5.2 2.9 3.9 ...
## $ Quality : num [1:38] 9.8 12.6 11.9 11.1 13.3 12.8 12.8 12 13.6 13.9 ...
## $ Region : num [1:38] 1 1 1 1 1 1 1 1 3 1 ...
summary(wine)
## Clarity Aroma Body Flavor
## Min. :0.5000 Min. :3.300 Min. :2.600 Min. :2.900
## 1st Qu.:0.8250 1st Qu.:4.125 1st Qu.:4.150 1st Qu.:4.225
## Median :1.0000 Median :4.650 Median :4.750 Median :4.800
## Mean :0.9237 Mean :4.847 Mean :4.684 Mean :4.768
## 3rd Qu.:1.0000 3rd Qu.:5.450 3rd Qu.:5.375 3rd Qu.:5.500
## Max. :1.0000 Max. :7.700 Max. :6.600 Max. :7.000
## Oakiness Quality Region
## Min. :2.900 Min. : 7.90 Min. :1.000
## 1st Qu.:3.700 1st Qu.:11.15 1st Qu.:1.000
## Median :4.100 Median :12.45 Median :2.000
## Mean :4.255 Mean :12.44 Mean :1.868
## 3rd Qu.:4.775 3rd Qu.:13.75 3rd Qu.:3.000
## Max. :6.000 Max. :16.10 Max. :3.000
## Matriks scatterplot semua variabel sekaligus
pairs(wine, main = "Matriks Scatterplot - Wine Quality Data",
pch = 19, col = "steelblue")
## Scatterplot Quality vs masing-masing prediktor + garis regresi
par(mfrow = c(2, 3))
prediktor <- c("Clarity", "Aroma", "Body", "Flavor", "Oakiness")
for (p in prediktor) {
plot(wine[[p]], wine$Quality,
xlab = p, ylab = "Quality",
main = paste("Quality vs", p,
"\n(r =", round(cor(wine[[p]], wine$Quality), 3), ")"),
pch = 19, col = "darkslategray")
abline(lm(wine$Quality ~ wine[[p]]), col = "red", lwd = 2)
}
par(mfrow = c(1, 1))
## Matriks korelasi
korelasi <- cor(wine)
round(korelasi, 3)
## Clarity Aroma Body Flavor Oakiness Quality Region
## Clarity 1.000 0.062 -0.308 -0.085 0.183 0.028 -0.020
## Aroma 0.062 1.000 0.549 0.737 0.202 0.707 0.620
## Body -0.308 0.549 1.000 0.647 0.152 0.549 0.435
## Flavor -0.085 0.737 0.647 1.000 0.180 0.790 0.506
## Oakiness 0.183 0.202 0.152 0.180 1.000 -0.047 -0.060
## Quality 0.028 0.707 0.549 0.790 -0.047 1.000 0.507
## Region -0.020 0.620 0.435 0.506 -0.060 0.507 1.000
## Heatmap korelasi (opsional, butuh package corrplot)
# install.packages("corrplot")
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
corrplot(korelasi, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", number.cex = 0.8,
title = "Matriks Korelasi Antarvariabel", mar = c(0,0,2,0))
model <- lm(Quality ~ Clarity + Aroma + Body + Flavor + Oakiness, data = wine)
summary(model)
##
## Call:
## lm(formula = Quality ~ Clarity + Aroma + Body + Flavor + Oakiness,
## data = wine)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.85552 -0.57448 -0.07092 0.67275 1.68093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9969 2.2318 1.791 0.082775 .
## Clarity 2.3395 1.7348 1.349 0.186958
## Aroma 0.4826 0.2724 1.771 0.086058 .
## Body 0.2732 0.3326 0.821 0.417503
## Flavor 1.1683 0.3045 3.837 0.000552 ***
## Oakiness -0.6840 0.2712 -2.522 0.016833 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.163 on 32 degrees of freedom
## Multiple R-squared: 0.7206, Adjusted R-squared: 0.6769
## F-statistic: 16.51 on 5 and 32 DF, p-value: 4.703e-08
coef(model)
## (Intercept) Clarity Aroma Body Flavor Oakiness
## 3.9968648 2.3394535 0.4825505 0.2731612 1.1683238 -0.6840102
## Persamaan model:
## Quality_hat = 4.2147 + 0.6998*Clarity + 0.6295*Aroma + 0.7058*Body
## + 0.7317*Flavor - 0.5295*Oakiness
## Selang kepercayaan koefisien
confint(model, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) -0.54910206 8.5428317
## Clarity -1.19427368 5.8731807
## Aroma -0.07240642 1.0375075
## Body -0.40424262 0.9505650
## Flavor 0.54811681 1.7885307
## Oakiness -1.23641174 -0.1316086
## Tabel ANAVA overall (regresi vs residual)
anova_overall <- anova(model)
print(anova_overall)
## Analysis of Variance Table
##
## Response: Quality
## Df Sum Sq Mean Sq F value Pr(>F)
## Clarity 1 0.125 0.125 0.0926 0.7628120
## Aroma 1 77.353 77.353 57.2351 1.286e-08 ***
## Body 1 6.414 6.414 4.7461 0.0368417 *
## Flavor 1 19.050 19.050 14.0953 0.0006946 ***
## Oakiness 1 8.598 8.598 6.3616 0.0168327 *
## Residuals 32 43.248 1.352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Ringkas jadi tabel Regresi - Residual - Total
SSR <- sum(anova_overall$"Sum Sq"[1:5]) # jumlah kuadrat regresi
SSE <- anova_overall$"Sum Sq"[6] # jumlah kuadrat residual
SST <- SSR + SSE
dfR <- sum(anova_overall$Df[1:5])
dfE <- anova_overall$Df[6]
dfT <- dfR + dfE
MSR <- SSR / dfR
MSE <- SSE / dfE
Fhitung <- MSR / MSE
pvalue <- 1 - pf(Fhitung, dfR, dfE)
tabel_anava <- data.frame(
Sumber = c("Regresi", "Residual", "Total"),
db = c(dfR, dfE, dfT),
JK = c(SSR, SSE, SST),
KT = c(MSR, MSE, NA),
Fhitung = c(Fhitung, NA, NA),
pvalue = c(pvalue, NA, NA)
)
print(tabel_anava, digits = 5)
## Sumber db JK KT Fhitung pvalue
## 1 Regresi 5 111.540 22.3081 16.506 4.703e-08
## 2 Residual 32 43.248 1.3515 NA NA
## 3 Total 37 154.788 NA NA NA
## Nilai F tabel pada alpha = 5%
alpha <- 0.05
F_tabel <- qf(1 - alpha, dfR, dfE)
cat("F_tabel(alpha=0.05, db1=", dfR, ", db2=", dfE, ") =", F_tabel, "\n")
## F_tabel(alpha=0.05, db1= 5 , db2= 32 ) = 2.512255
cat("F_hitung =", Fhitung, "\n")
## F_hitung = 16.50616
cat("Keputusan:", ifelse(Fhitung > F_tabel,
"Tolak H0 -> model signifikan",
"Gagal tolak H0 -> model tidak signifikan"), "\n")
## Keputusan: Tolak H0 -> model signifikan
## R-squared dan Adjusted R-squared
cat("R-squared :", summary(model)$r.squared, "\n")
## R-squared : 0.7205992
cat("Adj R-squared :", summary(model)$adj.r.squared, "\n")
## Adj R-squared : 0.6769428