# ==========================================
# 3. TIỀN XỬ LÝ SỐ LIỆU
# ==========================================
# 3.1 Đọc dữ liệu
printer_data <- read.csv("D:/HCMUT/HK252/BTL XSTK/data.csv")
head(printer_data, 3)
## layer_height wall_thickness infill_density infill_pattern nozzle_temperature
## 1 0.02 8 90 grid 220
## 2 0.02 7 90 honeycomb 225
## 3 0.02 1 80 grid 230
## bed_temperature print_speed material fan_speed roughness tension_strenght
## 1 60 40 abs 0 25 18
## 2 65 40 abs 25 32 16
## 3 70 40 abs 50 40 8
## elongation
## 1 1.2
## 2 1.4
## 3 0.8
# 3.2 Làm sạch dữ liệu
new_DF <- printer_data[,c("layer_height","wall_thickness","infill_density",
"infill_pattern","nozzle_temperature","bed_temperature","print_speed",
"material","fan_speed","tension_strenght")]
head(new_DF, 3)
## layer_height wall_thickness infill_density infill_pattern nozzle_temperature
## 1 0.02 8 90 grid 220
## 2 0.02 7 90 honeycomb 225
## 3 0.02 1 80 grid 230
## bed_temperature print_speed material fan_speed tension_strenght
## 1 60 40 abs 0 18
## 2 65 40 abs 25 16
## 3 70 40 abs 50 8
library(questionr)
freq.na(new_DF)
## missing %
## layer_height 0 0
## wall_thickness 0 0
## infill_density 0 0
## infill_pattern 0 0
## nozzle_temperature 0 0
## bed_temperature 0 0
## print_speed 0 0
## material 0 0
## fan_speed 0 0
## tension_strenght 0 0
# ==========================================
# 4. THỐNG KÊ MÔ TẢ
# ==========================================
new_function <- function(x){
c(n=length(x),xtb=mean(x),sd=sd(x),Q1=quantile(x,probs=0.25),
Q2 = median(x),Q3=quantile(x,probs=0.75),min=min(x),max=max(x))
}
continous_data <- new_DF[,c("layer_height","wall_thickness","infill_density",
"nozzle_temperature","bed_temperature","print_speed",
"fan_speed","tension_strenght")]
apply(continous_data, 2, new_function)
## layer_height wall_thickness infill_density nozzle_temperature
## n 50.00000000 50.000000 50.00000 50.00000
## xtb 0.10600000 5.220000 53.40000 221.50000
## sd 0.06439673 2.922747 25.36348 14.82035
## Q1.25% 0.06000000 3.000000 40.00000 210.00000
## Q2 0.10000000 5.000000 50.00000 220.00000
## Q3.75% 0.15000000 7.000000 80.00000 230.00000
## min 0.02000000 1.000000 10.00000 200.00000
## max 0.20000000 10.000000 90.00000 250.00000
## bed_temperature print_speed fan_speed tension_strenght
## n 50.000000 50.0000 50.00000 50.000000
## xtb 70.000000 64.0000 50.00000 20.080000
## sd 7.142857 29.6923 35.71429 8.925634
## Q1.25% 65.000000 40.0000 25.00000 12.000000
## Q2 70.000000 60.0000 50.00000 19.000000
## Q3.75% 75.000000 60.0000 75.00000 27.000000
## min 60.000000 40.0000 0.00000 4.000000
## max 80.000000 120.0000 100.00000 37.000000
table(new_DF$infill_pattern)
##
## grid honeycomb
## 25 25
table(new_DF$material)
##
## abs pla
## 25 25
hist(new_DF$tension_strenght, xlab="tension_strenght",
main="Histogram of tension_strenght", col="blue", labels=T, ylim=c(0,15))
boxplot(tension_strenght~infill_pattern, new_DF, col=c("red","blue"),
main="tension_strenght and infill_pattern")
boxplot(tension_strenght~material, new_DF, col=c("red","blue"),
main="tension_strenght and material")
par(mfrow = c(2,4))
plot(new_DF$layer_height, new_DF$tension_strenght, xlab="layer_height",
ylab="tension_strenght", main="layer_height & tension_strenght")
plot(new_DF$wall_thickness, new_DF$tension_strenght, xlab="wall_thickness",
ylab="tension_strenght", main="wall_thickness & tension_strenght")
plot(new_DF$infill_density, new_DF$tension_strenght, xlab="infill_density",
ylab="tension_strenght", main="infill_density & tension_strenght")
plot(new_DF$nozzle_temperature, new_DF$tension_strenght, xlab="nozzle_temperature",
ylab="tension_strenght", main="nozzle_temperature & tension_strenght")
plot(new_DF$bed_temperature, new_DF$tension_strenght, xlab="bed_temperature",
ylab="tension_strenght", main="bed_temperature & tension_strenght")
plot(new_DF$print_speed, new_DF$tension_strenght, xlab="print_speed",
ylab="tension_strenght", main="print_speed & tension_strenght")
plot(new_DF$fan_speed, new_DF$tension_strenght, xlab="fan_speed",
ylab="tension_strenght", main="fan_speed & tension_strenght")
# Khôi phục khung đồ thị
par(mfrow = c(1,1))
library(corrplot)
## corrplot 0.95 loaded
corrplot(cor(continous_data), method = "number")
# ==========================================
# 5.1 BÀI TOÁN MỘT MẪU
# ==========================================
n <- length(new_DF$tension_strenght)
xtb <- mean(new_DF$tension_strenght)
s <- sd(new_DF$tension_strenght)
data.frame(n, xtb, s)
## n xtb s
## 1 50 20.08 8.925634
qqnorm(new_DF$tension_strenght)
qqline(new_DF$tension_strenght)
shapiro.test(new_DF$tension_strenght)
##
## Shapiro-Wilk normality test
##
## data: new_DF$tension_strenght
## W = 0.9566, p-value = 0.06404
epsilon <- qt(p=.05/2, df=n-1, lower.tail=FALSE)*s/sqrt(n)
print(epsilon)
## [1] 2.536637
data.frame(u1 = xtb-epsilon, u2 = xtb+epsilon)
## u1 u2
## 1 17.54336 22.61664
# ==========================================
# 5.2 BÀI TOÁN HAI MẪU
# ==========================================
abs_data <- subset(new_DF, material=="abs")
qqnorm(abs_data$tension_strenght)
qqline(abs_data$tension_strenght)
shapiro.test(abs_data$tension_strenght)
##
## Shapiro-Wilk normality test
##
## data: abs_data$tension_strenght
## W = 0.90604, p-value = 0.02489
pla_data <- subset(new_DF, material=="pla")
qqnorm(pla_data$tension_strenght)
qqline(pla_data$tension_strenght)
shapiro.test(pla_data$tension_strenght)
##
## Shapiro-Wilk normality test
##
## data: pla_data$tension_strenght
## W = 0.92466, p-value = 0.06547
var.test(tension_strenght~material, new_DF, alternative="greater")
##
## F test to compare two variances
##
## data: tension_strenght by material
## F = 1.3045, num df = 24, denom df = 24, p-value = 0.26
## alternative hypothesis: true ratio of variances is greater than 1
## 95 percent confidence interval:
## 0.6575797 Inf
## sample estimates:
## ratio of variances
## 1.30448
t.test(tension_strenght~material, new_DF, var.equal = T)
##
## Two Sample t-test
##
## data: tension_strenght by material
## t = -2.0972, df = 48, p-value = 0.04126
## alternative hypothesis: true difference in means between group abs and group pla is not equal to 0
## 95 percent confidence interval:
## -10.028585 -0.211415
## sample estimates:
## mean in group abs mean in group pla
## 17.52 22.64
# ==========================================
# 5.3 ANOVA MỘT NHÂN TỐ (ONE-WAY ANOVA)
# ==========================================
# 1. Chia nhóm nozzle temperature (Đã gán thẳng thành Factor)
new_DF$nozzle_temperature_2 <- factor(
ifelse(new_DF$nozzle_temperature < 220, "Group_1",
ifelse(new_DF$nozzle_temperature > 230, "Group_3", "Group_2")),
levels = c("Group_1", "Group_2", "Group_3")
)
# 2. Boxplot
boxplot(
tension_strenght ~ nozzle_temperature_2,
data = new_DF,
main = "Tensile strength vs Nozzle temperature",
xlab = "Nozzle temperature group",
ylab = "Tensile strength",
col = "lightblue"
)
# 3. QQ-plot cho từng nhóm (Hiển thị ngang 3 biểu đồ)
par(mfrow = c(1,3))
by(
new_DF$tension_strenght,
new_DF$nozzle_temperature_2,
function(x){
qqnorm(x, main = "Normal Q-Q Plot")
qqline(x)
}
)
## new_DF$nozzle_temperature_2: Group_1
## NULL
## ------------------------------------------------------------
## new_DF$nozzle_temperature_2: Group_2
## NULL
## ------------------------------------------------------------
## new_DF$nozzle_temperature_2: Group_3
## NULL
par(mfrow = c(1,1)) # Trả lại khung mặc định
# 4. Kiểm định Shapiro–Wilk
by(
new_DF$tension_strenght,
new_DF$nozzle_temperature_2,
shapiro.test
)
## new_DF$nozzle_temperature_2: Group_1
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.91815, p-value = 0.09128
##
## ------------------------------------------------------------
## new_DF$nozzle_temperature_2: Group_2
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95684, p-value = 0.4829
##
## ------------------------------------------------------------
## new_DF$nozzle_temperature_2: Group_3
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.87204, p-value = 0.1056
# 5. Kiểm định Levene
library(car)
## Loading required package: carData
leveneTest(tension_strenght ~ nozzle_temperature_2, data = new_DF)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.9575 0.3912
## 47
# 6. One-way ANOVA
ANOVA_model <- aov(tension_strenght ~ nozzle_temperature_2, data = new_DF)
summary(ANOVA_model)
## Df Sum Sq Mean Sq F value Pr(>F)
## nozzle_temperature_2 2 688 343.9 5.026 0.0105 *
## Residuals 47 3216 68.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 7. Tính F_value và F_critical (miền bác bỏ)
anova_table <- summary(ANOVA_model)[[1]]
# Lưu ý: Tên biến trong bảng anova có thể có khoảng trắng, dùng vị trí dòng/cột để lấy giá trị cho an toàn
F_value <- anova_table[1, "F value"]
df1 <- anova_table[1, "Df"] # k - 1
df2 <- anova_table["Residuals", "Df"] # n - k
alpha <- 0.05
F_critical <- qf(1 - alpha, df1, df2)
cat("F_value =", F_value, "\n")
## F_value = 5.025912
cat("F_critical =", F_critical, "\n")
## F_critical = 3.195056
if (F_value > F_critical) {
cat("Kết luận: Bác bỏ H0 (F_value > F_critical)\n")
} else {
cat("Kết luận: Chưa đủ cơ sở bác bỏ H0 (F_value <= F_critical)\n")
}
## Kết luận: Bác bỏ H0 (F_value > F_critical)
# 8. Tukey HSD (Áp dụng căn chỉnh lề để nhãn trục y không bị mất)
TukeyHSD(ANOVA_model)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = tension_strenght ~ nozzle_temperature_2, data = new_DF)
##
## $nozzle_temperature_2
## diff lwr upr p adj
## Group_2-Group_1 -3.10 -9.430522 3.2305221 0.4678313
## Group_3-Group_1 -10.15 -17.903275 -2.3967255 0.0074607
## Group_3-Group_2 -7.05 -14.803275 0.7032745 0.0814970
par(mar = c(5, 8, 4, 2) + 0.1)
plot(TukeyHSD(ANOVA_model), las = 1)
par(mar = c(5, 4, 4, 2) + 0.1)
# ==========================================
# 5.4 HỒI QUY TUYẾN TÍNH ĐA BIẾN
# ==========================================
# 1. Xây dựng mô hình 1 (Đầy đủ các biến)
model_1 <- lm(tension_strenght ~ layer_height + wall_thickness + infill_density +
infill_pattern + nozzle_temperature + bed_temperature + print_speed + material,
data = new_DF)
summary(model_1)
##
## Call:
## lm(formula = tension_strenght ~ layer_height + wall_thickness +
## infill_density + infill_pattern + nozzle_temperature + bed_temperature +
## print_speed + material, data = new_DF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4101 -3.8491 0.0338 3.9073 13.5086
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 171.62809 54.21779 3.166 0.00292 **
## layer_height 55.59721 12.78771 4.348 8.87e-05 ***
## wall_thickness 1.06871 0.31942 3.346 0.00176 **
## infill_density 0.16286 0.03415 4.769 2.35e-05 ***
## infill_patternhoneycomb -1.14271 1.64581 -0.694 0.49140
## nozzle_temperature -1.04681 0.36901 -2.837 0.00705 **
## bed_temperature 1.00534 0.47426 2.120 0.04012 *
## print_speed -0.01559 0.03006 -0.519 0.60679
## materialpla -17.30508 8.51530 -2.032 0.04864 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.58 on 41 degrees of freedom
## Multiple R-squared: 0.673, Adjusted R-squared: 0.6092
## F-statistic: 10.55 on 8 and 41 DF, p-value: 6.91e-08
# 2. Xây dựng mô hình 2 (Đã loại bỏ các biến không có ý nghĩa thống kê)
model_2 <- lm(tension_strenght ~ layer_height + wall_thickness + infill_density +
nozzle_temperature + bed_temperature + material,
data = new_DF)
summary(model_2)
##
## Call:
## lm(formula = tension_strenght ~ layer_height + wall_thickness +
## infill_density + nozzle_temperature + bed_temperature + material,
## data = new_DF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4717 -4.1695 0.0914 3.8586 13.9581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.02463 53.19880 3.140 0.003055 **
## layer_height 56.36695 12.44845 4.528 4.67e-05 ***
## wall_thickness 1.11169 0.28024 3.967 0.000271 ***
## infill_density 0.16643 0.03339 4.985 1.06e-05 ***
## nozzle_temperature -1.02890 0.36321 -2.833 0.006996 **
## bed_temperature 0.98346 0.46685 2.107 0.041026 *
## materialpla -17.10365 8.39202 -2.038 0.047722 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.501 on 43 degrees of freedom
## Multiple R-squared: 0.6667, Adjusted R-squared: 0.6201
## F-statistic: 14.33 on 6 and 43 DF, p-value: 6.783e-09
# ==========================================
# CÁC KIỂM ĐỊNH GIẢ ĐỊNH CỦA MÔ HÌNH HỒI QUY
# ==========================================
# A. Kiểm định sai số có kỳ vọng bằng 0 và tính tuyến tính
t.test(model_2$residuals)
##
## One Sample t-test
##
## data: model_2$residuals
## t = -7.5559e-16, df = 49, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.464546 1.464546
## sample estimates:
## mean of x
## -5.506598e-16
plot(model_2, which = 1) # Đồ thị Residuals vs Fitted
# B. Kiểm định phương sai của các sai số là hằng số (Đồng phương sai)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
bptest(model_2)
##
## studentized Breusch-Pagan test
##
## data: model_2
## BP = 6.6773, df = 6, p-value = 0.3517
plot(model_2, which = 3) # Đồ thị Scale-Location
# C. Kiểm định sai số tuân theo phân phối chuẩn
shapiro.test(model_2$residuals)
##
## Shapiro-Wilk normality test
##
## data: model_2$residuals
## W = 0.97458, p-value = 0.3517
plot(model_2, which = 2) # Đồ thị Normal Q-Q
# D. Kiểm định các sai số độc lập với nhau
dwtest(model_2)
##
## Durbin-Watson test
##
## data: model_2
## DW = 1.5005, p-value = 0.01575
## alternative hypothesis: true autocorrelation is greater than 0
# Đồ thị kiểm tra tính độc lập của sai số theo thứ tự quan sát
plot(residuals(model_2), type = "o",
main = "Residuals vs Order",
xlab = "Observation Order",
ylab = "Residuals")
# ==========================================
# KIỂM TRA ĐIỂM NGOẠI LAI VÀ ĐA CỘNG TUYẾN
# ==========================================
# E. Kiểm tra giá trị ngoại lai (Outlier)
library(car)
outlierTest(model_2)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 11 2.910173 0.0057559 0.2878
plot(model_2, which = 5) # Đồ thị Residuals vs Leverage
# F. Kiểm tra đa cộng tuyến
library(psych)
##
## Attaching package: 'psych'
## The following object is masked from 'package:car':
##
## logit
## The following object is masked from 'package:questionr':
##
## describe
data_numeric <- model_2$model[, sapply(model_2$model, is.numeric)]
cortest.bartlett(data_numeric)
## R was not square, finding R from data
## $chisq
## [1] 75.23686
##
## $p.value
## [1] 5.130664e-10
##
## $df
## [1] 15
# Thiết lập khoảng cách lề dưới rộng ra để hiển thị tên biến không bị cắt
par(mar = c(8, 4, 4, 2) + 0.1)
# Đồ thị Barplot cho VIF
barplot(car::vif(model_2),
main = "Chỉ số VIF của các biến",
col = "steelblue",
las = 2) # las = 2 giúp xoay dọc tên biến để dễ đọc
# Thêm đường ngưỡng báo động (10)
abline(h = 10, col = "red", lwd = 2, lty = 2)
# Đồ thị Plot dạng đường (Line plot) cho VIF
vif_values <- vif(model_2)
plot(vif_values, type = "b",
main = "Chỉ số VIF của các biến",
xlab = "", ylab = "Giá trị VIF", xaxt = "n")
axis(1, at = 1:length(vif_values), labels = names(vif_values), las = 2)
abline(h = 10, col = "red", lty = 2)
# Trả lại lề mặc định sau khi vẽ xong
par(mar = c(5, 4, 4, 2) + 0.1)