MrXuan_Bao2
Dinh Phuong
2025-10-17
## Them font chu
library(showtext)## Loading required package: sysfonts## Loading required package: showtextdbfont_add("Times New Roman", "times.ttf") # Đường dẫn tới file font Times New Roman
font_add("Arial", "arial.ttf")
font_add("Times New Roman Bold", "timesbd.ttf")
showtext_auto()
# 1
library(ggplot2)
library(viridis)## Loading required package: viridisLitelibrary(ggthemes)QHTN (bố trí theo kiểu Box-wilson design) với các yếu tố lần lượt là Gocluon_R, KcTam_e, HsMasat_μ.
library(rsm)
# Box-Wilson Design bbd
QHTN <- bbd(
k = 3, # 3 yếu tố
n0 = 3, # 3 điểm trung tâm
block = FALSE, # No blocks
randomize = FALSE, # Not randomized
coding = list(
x1 ~ (KcTam_e - 25)/10,
x2 ~ (Gocluon_R - 5)/2,
x3 ~ (HsMasat_μ - 0.15)/0.05
)
)
print(QHTN)## run.order std.order KcTam_e Gocluon_R HsMasat_μ
## 1 1 1 15 3 0.15
## 2 2 2 35 3 0.15
## 3 3 3 15 7 0.15
## 4 4 4 35 7 0.15
## 5 5 5 15 5 0.10
## 6 6 6 35 5 0.10
## 7 7 7 15 5 0.20
## 8 8 8 35 5 0.20
## 9 9 9 25 3 0.10
## 10 10 10 25 7 0.10
## 11 11 11 25 3 0.20
## 12 12 12 25 7 0.20
## 13 13 13 25 5 0.15
## 14 14 14 25 5 0.15
## 15 15 15 25 5 0.15
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (KcTam_e - 25)/10
## x2 ~ (Gocluon_R - 5)/2
## x3 ~ (HsMasat_μ - 0.15)/0.05Bước 2: Thí nghiệm thực tế và đưa dữ liệu vào R
## Ở đây mình chạy số liệu mô phỏng về LucF và MucdoBD
library(readxl)
data <- read_excel("~/1.Detai_Caohoc/MrXuan/Bao2/DataQHTN_ECAPPC_211025.xlsx",
sheet = "BoxBehnken", col_types = c("numeric",
"numeric", "numeric", "numeric",
"numeric"))
head(data)## # A tibble: 6 × 5
## KCTam_e GocLuon_R HesoMS_f LucF_Y1 CiBD_Y2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 15 3 0.15 124. 0.499
## 2 35 3 0.15 131. 0.340
## 3 15 7 0.15 116. 0.485
## 4 35 7 0.15 117. 0.578
## 5 15 5 0.1 111. 0.544
## 6 35 5 0.1 113. 0.457# Kiểm tra số dòng của QHTN và data_excel
if (nrow(data) != nrow(QHTN)) {
stop("Số dòng của dữ liệu Excel không khớp với số dòng của QHTN!")}
# Gắn cột y1 và y2 từ data vào QHTN
QHTN$LucF_Y1 <- data$LucF_Y1
QHTN$CiBD_Y2 <- data$CiBD_Y2
# In kết quả để kiểm tra
print(QHTN)## run.order std.order KcTam_e Gocluon_R HsMasat_μ LucF_Y1 CiBD_Y2
## 1 1 1 15 3 0.15 123.7659 0.4994021
## 2 2 2 35 3 0.15 130.7823 0.3400517
## 3 3 3 15 7 0.15 116.1473 0.4849841
## 4 4 4 35 7 0.15 116.7619 0.5778713
## 5 5 5 15 5 0.10 110.9146 0.5444828
## 6 6 6 35 5 0.10 113.0822 0.4565588
## 7 7 7 15 5 0.20 133.1021 0.6148683
## 8 8 8 35 5 0.20 133.7514 0.5768149
## 9 9 9 25 3 0.10 121.8491 0.3216295
## 10 10 10 25 7 0.10 107.1278 0.4771377
## 11 11 11 25 3 0.20 133.2195 0.3620474
## 12 12 12 25 7 0.20 130.8387 0.5030914
## 13 13 13 25 5 0.15 124.2060 0.4579854
## 14 14 14 25 5 0.15 125.6337 0.4834678
## 15 15 15 25 5 0.15 127.0610 0.5089501
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (KcTam_e - 25)/10
## x2 ~ (Gocluon_R - 5)/2
## x3 ~ (HsMasat_μ - 0.15)/0.05Bước 3: Tính toán Hàm hôi quy & model tối ưu hóa #1: Hàm hồi quy đủ hệ số (Full model)
ket_qua_model <- rsm((LucF_Y1) ~ SO(x1, x2, x3),
data = QHTN)
summary(ket_qua_model)##
## Call:
## rsm(formula = (LucF_Y1) ~ SO(x1, x2, x3), data = QHTN)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.63356 1.06778 117.6587 8.411e-10 ***
## x1 1.30600 0.65388 1.9973 0.1022901
## x2 -4.84264 0.65388 -7.4060 0.0007064 ***
## x3 9.74223 0.65388 14.8991 2.465e-05 ***
## x1:x2 -1.60043 0.92472 -1.7307 0.1440592
## x1:x3 -0.37958 0.92472 -0.4105 0.6984543
## x2:x3 3.08516 0.92472 3.3363 0.0206360 *
## x1^2 -2.15771 0.96248 -2.2418 0.0750438 .
## x2^2 -1.61149 0.96248 -1.6743 0.1549243
## x3^2 -0.76328 0.96248 -0.7930 0.4637111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9837, Adjusted R-squared: 0.9545
## F-statistic: 33.63 on 9 and 5 DF, p-value: 0.000602
##
## Analysis of Variance Table
##
## Response: (LucF_Y1)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 960.54 320.18 93.6077 8.194e-05
## TWI(x1, x2, x3) 3 48.89 16.30 4.7649 0.06283
## PQ(x1, x2, x3) 3 25.86 8.62 2.5204 0.17199
## Residuals 5 17.10 3.42
## Lack of fit 3 13.03 4.34 2.1309 0.33523
## Pure error 2 4.08 2.04
##
## Stationary point of response surface:
## x1 x2 x3
## 1.735140 -3.565064 -1.254566
##
## Stationary point in original units:
## KcTam_e Gocluon_R HsMasat_μ
## 42.35139736 -2.13012879 0.08727168
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.5655325 -1.9867145 -3.1112960
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2338790 0.8179920 -0.5255375
## x2 -0.6178802 -0.2922928 -0.7299240
## x3 -0.7506828 0.4954331 0.4370600# pred.r_squared
PRESS <- sum((residuals(ket_qua_model)/(1-lm.influence(ket_qua_model)$hat))^2)
SST <- sum((ket_qua_model$fitted.values + residuals(ket_qua_model) - mean(ket_qua_model$model[,1]))^2)
predR2 <- 1 - PRESS/SST
#cat("Predicted R-squared:", round(predR2, 4), "\n")
cat(sprintf("Std. Dev.: %.4f Mean: %.4f C.V.%%: %.2f%% Predicted R-squared: %.4f\n", sigma(ket_qua_model), mean(ket_qua_model$model[,1]), sigma(ket_qua_model)/mean(ket_qua_model$model[,1])*100, round(predR2, 4)))## Std. Dev.: 1.8494 Mean: 123.2162 C.V.%: 1.50% Predicted R-squared: 0.7932#2: Hàm hồi quy (đã kiểm định Studien Ttest và Fisher Ftest)
#Hoặc
ket_qua_model <- rsm((LucF_Y1) ~ FO(x1, x2, x3)+TWI(x1,x2,x3)+PQ(x1,x2,x3),
data = QHTN)
summary(ket_qua_model)##
## Call:
## rsm(formula = (LucF_Y1) ~ FO(x1, x2, x3) + TWI(x1, x2, x3) +
## PQ(x1, x2, x3), data = QHTN)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.63356 1.06778 117.6587 8.411e-10 ***
## x1 1.30600 0.65388 1.9973 0.1022901
## x2 -4.84264 0.65388 -7.4060 0.0007064 ***
## x3 9.74223 0.65388 14.8991 2.465e-05 ***
## x1:x2 -1.60043 0.92472 -1.7307 0.1440592
## x1:x3 -0.37958 0.92472 -0.4105 0.6984543
## x2:x3 3.08516 0.92472 3.3363 0.0206360 *
## x1^2 -2.15771 0.96248 -2.2418 0.0750438 .
## x2^2 -1.61149 0.96248 -1.6743 0.1549243
## x3^2 -0.76328 0.96248 -0.7930 0.4637111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9837, Adjusted R-squared: 0.9545
## F-statistic: 33.63 on 9 and 5 DF, p-value: 0.000602
##
## Analysis of Variance Table
##
## Response: (LucF_Y1)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 960.54 320.18 93.6077 8.194e-05
## TWI(x1, x2, x3) 3 48.89 16.30 4.7649 0.06283
## PQ(x1, x2, x3) 3 25.86 8.62 2.5204 0.17199
## Residuals 5 17.10 3.42
## Lack of fit 3 13.03 4.34 2.1309 0.33523
## Pure error 2 4.08 2.04
##
## Stationary point of response surface:
## x1 x2 x3
## 1.735140 -3.565064 -1.254566
##
## Stationary point in original units:
## KcTam_e Gocluon_R HsMasat_μ
## 42.35139736 -2.13012879 0.08727168
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.5655325 -1.9867145 -3.1112960
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2338790 0.8179920 -0.5255375
## x2 -0.6178802 -0.2922928 -0.7299240
## x3 -0.7506828 0.4954331 0.4370600# Thực hiện stepwise regression dựa trên AIC
# Sử dụng direction = "both" để thử cả forward và backward selection
ket_qua_model_opt <- step(ket_qua_model, direction = "both", trace = 1)## Start: AIC=21.97
## (LucF_Y1) ~ FO(x1, x2, x3) + TWI(x1, x2, x3) + PQ(x1, x2, x3)
##
## Df Sum of Sq RSS AIC
## <none> 17.10 21.967
## - PQ(x1, x2, x3) 3 25.86 42.96 29.785
## - TWI(x1, x2, x3) 3 48.89 66.00 36.223
## - FO(x1, x2, x3) 3 960.54 977.65 76.656# Hiển thị tóm tắt mô hình tối ưu
summary(ket_qua_model_opt)##
## Call:
## rsm(formula = (LucF_Y1) ~ FO(x1, x2, x3) + TWI(x1, x2, x3) +
## PQ(x1, x2, x3), data = QHTN)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.63356 1.06778 117.6587 8.411e-10 ***
## x1 1.30600 0.65388 1.9973 0.1022901
## x2 -4.84264 0.65388 -7.4060 0.0007064 ***
## x3 9.74223 0.65388 14.8991 2.465e-05 ***
## x1:x2 -1.60043 0.92472 -1.7307 0.1440592
## x1:x3 -0.37958 0.92472 -0.4105 0.6984543
## x2:x3 3.08516 0.92472 3.3363 0.0206360 *
## x1^2 -2.15771 0.96248 -2.2418 0.0750438 .
## x2^2 -1.61149 0.96248 -1.6743 0.1549243
## x3^2 -0.76328 0.96248 -0.7930 0.4637111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9837, Adjusted R-squared: 0.9545
## F-statistic: 33.63 on 9 and 5 DF, p-value: 0.000602
##
## Analysis of Variance Table
##
## Response: (LucF_Y1)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 960.54 320.18 93.6077 8.194e-05
## TWI(x1, x2, x3) 3 48.89 16.30 4.7649 0.06283
## PQ(x1, x2, x3) 3 25.86 8.62 2.5204 0.17199
## Residuals 5 17.10 3.42
## Lack of fit 3 13.03 4.34 2.1309 0.33523
## Pure error 2 4.08 2.04
##
## Stationary point of response surface:
## x1 x2 x3
## 1.735140 -3.565064 -1.254566
##
## Stationary point in original units:
## KcTam_e Gocluon_R HsMasat_μ
## 42.35139736 -2.13012879 0.08727168
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.5655325 -1.9867145 -3.1112960
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2338790 0.8179920 -0.5255375
## x2 -0.6178802 -0.2922928 -0.7299240
## x3 -0.7506828 0.4954331 0.4370600# Lấy phương trình hồi quy từ mô hình tối ưu
formula_opt <- formula(ket_qua_model_opt)
cat("Phương trình hồi quy tối ưu:\n", as.character(formula_opt)[3], "\n")## Phương trình hồi quy tối ưu:
## FO(x1, x2, x3) + TWI(x1, x2, x3) + PQ(x1, x2, x3)# Kiểm tra các hệ số có ý nghĩa (p-value < 0.05)
coef_summary <- summary(ket_qua_model_opt)$coefficients
significant_coefs <- coef_summary[coef_summary[, "Pr(>|t|)"] < 0.05, , drop = FALSE]
print("Các hệ số có ý nghĩa thống kê (p-value < 0.05):")## [1] "Các hệ số có ý nghĩa thống kê (p-value < 0.05):"print(significant_coefs)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.633561 1.0677792 117.658747 8.410982e-10
## x2 -4.842645 0.6538785 -7.406031 7.064389e-04
## x3 9.742233 0.6538785 14.899148 2.464679e-05
## x2:x3 3.085157 0.9247239 3.336301 2.063595e-02par(mfrow = c(2, 3)) # 2 x 3 pictures on one plot
persp(
ket_qua_model, # Our model
~ x1 + x2 + x3, # A formula to obtain the 6 possible graphs
col = topo.colors(100), # Color palette
contours = "colors", # Include contours with the same color palette
# xlabs = c("Nhiệt độ (°C)",
# "Glucse (g/L)",
# "Peptone (g/L)",
# "Thời gian (giờ)" ),
# zlab = "Mật độ tế bào (log CFU/mL)",
expand = 1,
cex = 0.8
)## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter
## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter
## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter
Đồ thị tương quan giữa x1 và x2
persp(
ket_qua_model, # Our model
~ x2 + x3, # A formula to obtain the 6 possible graphs
col = topo.colors(100), # Color palette
contours = "colors", # Include contours with the same color palette
xlabs = c("\n\nNhiệt độ (°C)",
"\n\nGlucose (g/L)"),
# zlab = "\n\nMật độ tế bào (log CFU/mL)", # Dùng \n để cách hàng
cex.lab = 0.9,
expand = 0.8,
# theta = 30, phi = 30 # Xoay đồ thị,
ticktype = "detailed",
nticks = 4
)## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "nticks" is not a graphical parameterpar(xpd = NA, srt = 95) ## disable clipping and set string rotation
text(-0.28, 0.02,"Mật độ tế bào (log CFU/mL)")
################################################################################
# Biểu đồ 7
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(viridis)
library(showtext)
# Font chữ
font_add("Times New Roman", "times.ttf")
showtext_auto()
# --- Tạo lưới x2, x3 ở dạng mã hoá ---
x2_seq_coded <- seq(min(QHTN$x2), max(QHTN$x2), length.out = 80)
x3_seq_coded <- seq(min(QHTN$x3), max(QHTN$x3), length.out = 80)
R_seq_coded <- seq(min(data$GocLuon_R), max(data$GocLuon_R), length.out = 80)
f_seq_coded <- seq(min(data$HesoMS_f), max(data$HesoMS_f), length.out = 80)
grid <- expand.grid(x1 = 0, x3 = x2_seq_coded, x2 = x3_seq_coded)
grid_real <-expand.grid(x1 = 25, x2 = R_seq_coded, x3 = f_seq_coded)
grid$F_pred <- predict(ket_qua_model, newdata = grid)
# --- Chuyển thành ma trận z ---
z_mat <- matrix(grid$F_pred, nrow = length(x2_seq_coded), ncol = length(x3_seq_coded))
# --- Tick marks tự động theo giá trị thực ---
x2_ticks <- pretty(grid_real$x2)
x3_ticks <- pretty(grid_real$x3)
z_ticks <- c(110,115,120,125,130,135) # chỉ hiện 110->135
# format 2 chữ số thập phân
x3_ticktext <- sprintf("%.2f", x3_ticks)
# --- Surface 3D ---
fig <- plot_ly(
x = unique(grid_real$x2),
y = unique(grid_real$x3),
z = z_mat,
type = "surface",
colors = viridis(120, option = "turbo"),
showscale = FALSE,
colorbar = list(
len = 0.3, # ← chiều CAO thanh màu (0–1)
thickness = 25, # ← độ RỘNG thanh màu (pixels)
x = 0.85, # ← vị trí ngang (1 = mép phải đồ thị)
y = 0.4, # ← vị trí dọc (0=bottom, 1=top)
yanchor = "middle", # giữ tâm thanh
# ---- ÉP plotly phải hiển thị đầy đủ text ----
tickmode = "array",
tickvals = z_ticks,
ticktext = z_ticks,
tickfont = list(family = "Times New Roman", size = 12),
# ---- Title cho colorbar ----
title = list(
text = "Force F (ton)",
font = list(family = "Times New Roman", size = 12),
side = "top", # ← Thêm dòng này: vị trí title (top/bottom)
standoff = 20 # ← Thêm dòng này: khoảng cách từ title đến colorbar (pixels)
)),
contours = list(
x = list(show=TRUE, start = min(R_seq_coded), end = max(R_seq_coded), size = 0.2),
y = list(show=TRUE, start = min(f_seq_coded), end = max(f_seq_coded), size = 0.005)))
# --- Layout: tick marks, labels, font, box ---
fig <- fig %>% layout(
scene = list(
xaxis = list(
title = "Corner radius R (mm)",
tickmode = "array",
tickvals = x2_ticks,
ticktext = x2_ticks, # hiển thị giá trị thực
tickangle = -30, # xoay tick-label
showline = TRUE, linewidth = 2, mirror = TRUE,
ticks = "outside", # tick marks
ticklen = 5, # chiều dài tick
titlefont = list(family="Times New Roman", size=12),
tickfont = list(family="Times New Roman", size=12)
),
yaxis = list(
title = "Friction coefficient μ",
tickmode = "array",
tickvals = x3_ticks,
ticktext = x3_ticktext,
tickangle = -20,
showline = TRUE, linewidth = 2, mirror = TRUE,
ticks = "outside", ticklen = 5,
titlefont = list(family="Times New Roman", size=12),
tickfont = list(family="Times New Roman", size=12)
),
zaxis = list(
title = "Force F (ton)",
tickmode = "array",
tickvals = z_ticks, # chỉ hiện các giá trị mong muốn
ticktext = z_ticks,
showline = TRUE, linewidth = 2, mirror = TRUE,
ticks = "outside", ticklen = 5,
titlefont = list(family="Times New Roman", size=12),
tickfont = list(family="Times New Roman", size=12)
),
aspectmode = "cube", # hộp chữ nhật
bgcolor = "white",
# --- Thêm fix camera ---
camera = list(
eye = list(
x = 3 * cos(-250*pi/360) * cos(50*pi/360), # scale 2 để vừa khít
y = 3 * sin(-250*pi/360) * cos(50*pi/360),
z = 3 * sin(50*pi/360)
)
)
),
font = list(family = "Times New Roman"),
dragmode = FALSE, # không cho xoay/pan/zoom
margin = list(
l = 0, # left margin (tương ứng với 4 trong mar)
r = 2, # right margin (tương ứng với 0.5 trong mar)
t = 0, # top margin (tương ứng với 1 trong mar)
b = 2 # bottom margin (tương ứng với 3.5 trong mar)
)
#annotations = list(
#list(
#text = "Slice at Distance e = 25 (mm)",
#xref = "paper",
#yref = "paper",
#x = 0.4,
#y = 0.08,
#xanchor = "left",
#yanchor = "top",
#font = list(family = "Times New Roman", size = 12),
#showarrow = FALSE
#bgcolor = "white",
#bordercolor = "black",
#borderwidth = 1
)
figĐồ thị tương quan giữa x1 và x2 (#2)
# Load thư viện cần thiết
library(rsm)
library(plotly)
# Tạo dữ liệu lưới để dự đoán y dựa trên x1 và x2, giữ x3 cố định
x2_range <- seq(min(QHTN$x2), max(QHTN$x2), length.out = 200)
x3_range <- seq(min(QHTN$x3), max(QHTN$x3), length.out = 200)
x1_fixed <- 0 # Giữ x3 cố định tại -1
# Tạo lưới dữ liệu
grid <- expand.grid(x1 = x1_fixed, x2 = x2_range, x3 = x3_range)
y_pred <- predict(ket_qua_model, newdata = grid)
# Chuyển đổi y_pred thành ma trận để vẽ surface plot
y_matrix <- matrix(y_pred, nrow = length(x2_range), ncol = length(x3_range))
# Cách 1 Tạo các điểm chia cho trục x (10 khoảng)
#x_ticks <- seq(-1.215, 1.215, length.out = 11) # 11 điểm để tạo 10 khoảng
# Cách 2 Xác định các giá trị nhãn cụ thể cho trục x
x_ticks <- c(-1, -0.5, 0, 0.5, 1)
# Tùy chỉnh nhãn cho trục x
x_ticktext <- c("-1", "-0.5", "0", "0.5", "1")
# Vẽ biểu đồ 3D bằng plotly
plot_ly(x = x2_range, y = x3_range, z = y_matrix, type = "surface") %>%
layout(
title = "Ảnh hưởng của x1 và x2 đến y khi x3 = 0",
scene = list(
xaxis = list(title = "x1",range = c(0, 1), # Thiết lập phạm vi trục x
tickvals = x_ticks, # Các điểm chia trên trục x
ticktext = x_ticktext),
yaxis = list(title = "x2",range = c(-1, 1), # Thiết lập phạm vi trục x
tickvals = x_ticks, # Các điểm chia trên trục x
ticktext = x_ticktext),
zaxis = list(title = "y", range = c(min(y_pred), max(y_pred))),
camera = list(eye = list(x = -1.5, y = -1.5, z = 1.5))
)
)Đồ thị đường đồng mức
par(mfrow = c(2,3)) # 2 x 3 pictures on one plot
contour(
ket_qua_model, # Our model
~ x1 + x2 + x3, # A formula to obtain the 6 possible graphs
image = TRUE, # If image = TRUE, apply color to each contour
)
Tìm điểm tối ưu theo Central Composite Design chuẩn (CCD)
opt_point <- summary(ket_qua_model)$canonical$xs
# opt_point
op_point_ru <- code2val(
opt_point, # Optimal point in coded units
codings = codings(QHTN) # Formulas to convert to factor units
)
op_point_ru## KcTam_e Gocluon_R HsMasat_μ
## 42.35139736 -2.13012879 0.08727168Tìm điểm tối ưu theo QH Trực giao Cấp 2 (BBD)
# Lấy điểm tối ưu trong đơn vị mã hóa (coded units) với giới hạn
# Sử dụng optim để tìm điểm tối ưu lớn nhất trong khoảng (-1, 1)
objective <- function(x) {
-predict(ket_qua_model, newdata = data.frame(x1 = x[1], x2 = x[2], x3 = x[3]))
}
# Tìm điểm tối ưu lớn nhất với giới hạn
opt_point <- optim(
par = c(0, 0, 0), # Điểm khởi đầu
fn = objective, # Hàm mục tiêu (lấy âm để tìm max)
method = "L-BFGS-B", # Phương pháp với ràng buộc
lower = c(-1, -1, -1), # Giới hạn dưới
upper = c(1, 1, 1) # Giới hạn trên
)$par
cat("Điểm tối ưu trong đơn vị mã hóa (coded units):", round(opt_point, 3), "\n")## Điểm tối ưu trong đơn vị mã hóa (coded units): 0.511 -0.799 1# Chuyển đổi sang đơn vị thực tế (factor units)
op_point_ru <- code2val(
opt_point, # Optimal point in coded units
codings = codings(QHTN) # Formulas to convert to factor units
)
cat("Điểm tối ưu trong đơn vị thực tế:", round(op_point_ru, 3), "\n")## Điểm tối ưu trong đơn vị thực tế: 0.511 -0.799 1# Giá trị y tại điểm tối ưu
y_opt <- predict(ket_qua_model, newdata = data.frame(x1 = opt_point[1], x2 = opt_point[2], x3 = opt_point[3]))
cat("Giá trị y tại điểm tối ưu:", round(y_opt, 2), "\n")## Giá trị y tại điểm tối ưu: 135.55Dự đoán kết quả đầu ra từ điểm tối ưu
opt_point_df <- data.frame( # predict() needs a data frame with the points
x1 = opt_point[1], # to be predicted
x2 = opt_point[2],
x3 = opt_point[3]
)
best_response <- predict(
ket_qua_model, # Our model
opt_point_df # Data frame with points to be predicted
)
names(best_response) <- "Luc F"
best_response## Luc F
## 135.5514