# ==========================================
# 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)


# ==========================================
# TÌM MÔ HÌNH HỒI QUY TỐI ƯU CHO ROUGHNESS & ELONGATION (CÁCH 2)
# ==========================================
# 1. Mô hình dự báo Độ nhám (roughness)
model_full_roughness <- lm(roughness ~ layer_height + wall_thickness + infill_density + 
                             infill_pattern + nozzle_temperature + bed_temperature + 
                             print_speed + material, data = printer_data)

best_model_roughness <- step(model_full_roughness)
## Start:  AIC=372.48
## roughness ~ layer_height + wall_thickness + infill_density + 
##     infill_pattern + nozzle_temperature + bed_temperature + print_speed + 
##     material
## 
##                      Df Sum of Sq    RSS    AIC
## - infill_pattern      1         0  59968 370.48
## - infill_density      1        48  60015 370.52
## - wall_thickness      1      1663  61630 371.84
## <none>                             59968 372.48
## - print_speed         1     14538  74505 381.33
## - bed_temperature     1     36006  95973 393.99
## - material            1     38246  98213 395.14
## - nozzle_temperature  1     51832 111800 401.62
## - layer_height        1    306817 366785 461.03
## 
## Step:  AIC=370.48
## roughness ~ layer_height + wall_thickness + infill_density + 
##     nozzle_temperature + bed_temperature + print_speed + material
## 
##                      Df Sum of Sq    RSS    AIC
## - infill_density      1        48  60016 368.52
## - wall_thickness      1      1697  61665 369.87
## <none>                             59968 370.48
## - print_speed         1     14580  74548 379.36
## - bed_temperature     1     36115  96083 392.05
## - material            1     38274  98241 393.16
## - nozzle_temperature  1     51990 111957 399.69
## - layer_height        1    307226 367193 459.08
## 
## Step:  AIC=368.52
## roughness ~ layer_height + wall_thickness + nozzle_temperature + 
##     bed_temperature + print_speed + material
## 
##                      Df Sum of Sq    RSS    AIC
## - wall_thickness      1      1651  61667 367.87
## <none>                             60016 368.52
## - print_speed         1     14651  74667 377.44
## - bed_temperature     1     37456  97472 390.76
## - material            1     39029  99045 391.57
## - nozzle_temperature  1     54242 114258 398.71
## - layer_height        1    307270 367285 457.09
## 
## Step:  AIC=367.87
## roughness ~ layer_height + nozzle_temperature + bed_temperature + 
##     print_speed + material
## 
##                      Df Sum of Sq    RSS    AIC
## <none>                             61667 367.87
## - print_speed         1     13210  74877 375.58
## - bed_temperature     1     36693  98360 389.22
## - material            1     38454 100121 390.11
## - nozzle_temperature  1     53228 114895 396.99
## - layer_height        1    314770 376437 456.32
# Xem kết quả phương trình tối ưu của roughness
summary(best_model_roughness)
## 
## Call:
## lm(formula = roughness ~ layer_height + nozzle_temperature + 
##     bed_temperature + print_speed + material, data = printer_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -74.084 -26.500  -1.662  22.585  92.356 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -2310.7356   353.2009  -6.542 5.38e-08 ***
## layer_height        1246.5353    83.1780  14.986  < 2e-16 ***
## nozzle_temperature    14.7774     2.3979   6.163 1.95e-07 ***
## bed_temperature      -15.8078     3.0895  -5.117 6.55e-06 ***
## print_speed            0.5538     0.1804   3.070  0.00366 ** 
## materialpla          294.1610    56.1586   5.238 4.38e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37.44 on 44 degrees of freedom
## Multiple R-squared:  0.8717, Adjusted R-squared:  0.8571 
## F-statistic: 59.78 on 5 and 44 DF,  p-value: < 2.2e-16
# 2. Mô hình dự báo Độ giãn dài (elongation)
model_full_elongation <- lm(elongation ~ layer_height + wall_thickness + infill_density + 
                              infill_pattern + nozzle_temperature + bed_temperature + 
                              print_speed + material, data = printer_data)

best_model_elongation <- step(model_full_elongation)
## Start:  AIC=-69.64
## elongation ~ layer_height + wall_thickness + infill_density + 
##     infill_pattern + nozzle_temperature + bed_temperature + print_speed + 
##     material
## 
##                      Df Sum of Sq     RSS     AIC
## - infill_pattern      1    0.0446  8.7088 -71.384
## - wall_thickness      1    0.2900  8.9542 -69.995
## <none>                             8.6643 -69.641
## - print_speed         1    0.3924  9.0567 -69.426
## - material            1    1.3661 10.0304 -64.320
## - bed_temperature     1    1.5036 10.1678 -63.640
## - infill_density      1    2.6119 11.2762 -58.467
## - nozzle_temperature  1    2.7487 11.4130 -57.864
## - layer_height        1    7.7548 16.4191 -39.679
## 
## Step:  AIC=-71.38
## elongation ~ layer_height + wall_thickness + infill_density + 
##     nozzle_temperature + bed_temperature + print_speed + material
## 
##                      Df Sum of Sq     RSS     AIC
## - wall_thickness      1    0.2624  8.9713 -71.900
## <none>                             8.7088 -71.384
## - print_speed         1    0.4092  9.1180 -71.089
## - material            1    1.3548 10.0636 -66.155
## - bed_temperature     1    1.4804 10.1893 -65.534
## - nozzle_temperature  1    2.7193 11.4281 -59.797
## - infill_density      1    2.7397 11.4485 -59.708
## - layer_height        1    7.7218 16.4306 -41.644
## 
## Step:  AIC=-71.9
## elongation ~ layer_height + infill_density + nozzle_temperature + 
##     bed_temperature + print_speed + material
## 
##                      Df Sum of Sq     RSS     AIC
## <none>                             8.9713 -71.900
## - print_speed         1    0.9175  9.8887 -69.031
## - material            1    1.4352 10.4065 -66.480
## - bed_temperature     1    1.5981 10.5694 -65.703
## - nozzle_temperature  1    2.9072 11.8785 -59.865
## - infill_density      1    3.0158 11.9871 -59.410
## - layer_height        1    7.4871 16.4584 -43.559
# Xem kết quả phương trình tối ưu của elongation
summary(best_model_elongation)
## 
## Call:
## lm(formula = elongation ~ layer_height + infill_density + nozzle_temperature + 
##     bed_temperature + print_speed + material, data = printer_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.12952 -0.32192  0.05021  0.30777  0.83809 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        19.023103   4.384019   4.339 8.51e-05 ***
## layer_height        6.079479   1.014849   5.991 3.77e-07 ***
## infill_density      0.010456   0.002750   3.802 0.000448 ***
## nozzle_temperature -0.112096   0.030029  -3.733 0.000551 ***
## bed_temperature     0.106920   0.038632   2.768 0.008295 ** 
## print_speed        -0.004639   0.002212  -2.097 0.041909 *  
## materialpla        -1.824195   0.695511  -2.623 0.012016 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4568 on 43 degrees of freedom
## Multiple R-squared:  0.7053, Adjusted R-squared:  0.6642 
## F-statistic: 17.15 on 6 and 43 DF,  p-value: 5.335e-10