#Việc 1. Phân tích phương sai
#1.1 Giả sử bạn có dữ liệu về cân nặng của 4 nhóm A, B, C và D như sau:
#Nhóm A (n= 7): 8, 9, 11, 4, 7, 8, 5
#Nhóm B (n= 8): 7, 17, 10, 14, 12, 24, 11, 22
#Nhóm C (n= 6): 28, 21, 26, 11, 24, 19
#Nhóm D (n= 9): 26, 16, 13, 12, 9, 10, 11, 17, 15
#Nhập dữ liệu vào R và đặt tên dữ liệu là data.
# Tạo vector cân nặng cho từng nhóm
group_A <- c(8, 9, 11, 4, 7, 8, 5)
group_B <- c(7, 17, 10, 14, 12, 24, 11, 22)
group_C <- c(28, 21, 26, 11, 24, 19)
group_D <- c(26, 16, 13, 12, 9, 10, 11, 17, 15)
# Gộp lại thành 1 vector duy nhất
weight <- c(group_A, group_B, group_C, group_D)
# Tạo vector nhóm tương ứng
group <- c(rep("A", length(group_A)),
rep("B", length(group_B)),
rep("C", length(group_C)),
rep("D", length(group_D)))
# Tạo data frame
data <- data.frame(
group = factor(group),
weight = weight
)
# Xem dữ liệu
print(data)## group weight
## 1 A 8
## 2 A 9
## 3 A 11
## 4 A 4
## 5 A 7
## 6 A 8
## 7 A 5
## 8 B 7
## 9 B 17
## 10 B 10
## 11 B 14
## 12 B 12
## 13 B 24
## 14 B 11
## 15 B 22
## 16 C 28
## 17 C 21
## 18 C 26
## 19 C 11
## 20 C 24
## 21 C 19
## 22 D 26
## 23 D 16
## 24 D 13
## 25 D 12
## 26 D 9
## 27 D 10
## 28 D 11
## 29 D 17
## 30 D 15#1.2 Mô tả cân nặng giữa 4 nhóm.
# Nếu chưa có thì cài và gọi thư viện
# Tóm tắt cân nặng theo nhóm
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondata %>%
group_by(group) %>%
summarise(
count = n(),
mean = mean(weight),
median = median(weight),
sd = sd(weight),
min = min(weight),
max = max(weight),
Q1 = quantile(weight, 0.25),
Q3 = quantile(weight, 0.75)
)## # A tibble: 4 × 9
## group count mean median sd min max Q1 Q3
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A 7 7.43 8 2.37 4 11 6 8.5
## 2 B 8 14.6 13 5.95 7 24 10.8 18.2
## 3 C 6 21.5 22.5 6.09 11 28 19.5 25.5
## 4 D 9 14.3 13 5.15 9 26 11 16#Mô tả bằng biểu đồ (boxplot)
boxplot(weight ~ group, data = data,
main = "Boxplot of Weight by Group",
xlab = "Group", ylab = "Weight",
col = c("skyblue", "lightgreen", "salmon", "khaki"))
#1.3 Phân tích sự khác biệt về cân nặng giữa 4 nhóm. Diễn giải kết quả.
#a. Kiểm tra phân phối chuẩn từng nhóm (Shapiro-Wilk test):
by(data$weight, data$group, shapiro.test)## data$group: A
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96875, p-value = 0.8893
##
## ------------------------------------------------------------
## data$group: B
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.93765, p-value = 0.5881
##
## ------------------------------------------------------------
## data$group: C
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.93365, p-value = 0.6086
##
## ------------------------------------------------------------
## data$group: D
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.86649, p-value = 0.1127#Ketluan: Dữ liệu KHÔNG phân phối chuẩn vì p-vaue<0.05
# → Dùng Kruskal-Wallis test (thay thế ANOVA)
kruskal.test(weight ~ group, data = data)##
## Kruskal-Wallis rank sum test
##
## data: weight by group
## Kruskal-Wallis chi-squared = 15.35, df = 3, p-value = 0.001541 #Ketluan: p-value = 0.001541 < 0.05 → Có sự khác biệt có ý nghĩa thống kê về cân nặng giữa ít nhất 2 trong 4 nhóm.
#Kết quả này không nói nhóm nào khác biệt với nhóm nào → bạn cần phân tích tiếp post-hoc test. #1.4 Thực hiện phân tích hậu định (post-hoc analysis) để xác định cụ thể nhóm có khác biệt về cân nặng. Diễn giải kết quả.
#So sánh cặp đôi (post-hoc)
library(FSA)## ## FSA v0.10.0. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book. dunnTest(weight ~ group, data = data, method = "bonferroni")## Dunn (1964) Kruskal-Wallis multiple comparison## p-values adjusted with the Bonferroni method.## Comparison Z P.unadj P.adj
## 1 A - B -2.45056032 0.0142634064 0.0855804382
## 2 A - C -3.83109442 0.0001275745 0.0007654472
## 3 B - C -1.59822234 0.1099935038 0.6599610229
## 4 A - D -2.57941358 0.0098968218 0.0593809311
## 5 B - D -0.06506684 0.9481207798 1.0000000000
## 6 C - D 1.57770190 0.1146341103 0.6878046619 #Ketluan: Kết quả kiểm định hậu kiểm cho thấy chỉ có nhóm A và nhóm C có sự khác biệt có ý nghĩa thống kê về cân nặng (p.adj = 0.00076 < 0.05).
#Các cặp còn lại không có sự khác biệt đáng kể sau khi điều chỉnh Bonferroni.# Việc 2. Phân tích tương quan
# 2.1 Đọc dữ liệu "Demo data.csv" vào R và gọi dữ liệu là "df"
df = read.csv("D:\\2025-10-Tap huan_Phan tich du lieu va cong bo Quoc te danh cho giang vien\\Demo.csv")#2.2 Mô tả đặc điểm cân nặng (weight) và chiều cao (height).
# Tạo bảng mô tả cho cân nặng và chiều cao toàn bộ dữ liệu
library(gtsummary)
df %>%
tbl_summary(by = gender, statistic = list(all_continuous() ~ "{mean} ± {sd}"))| Characteristic | F N = 8621 |
M N = 3551 |
|---|---|---|
| X | 632 ± 359 | 553 ± 326 |
| id | 638 ± 362 | 558 ± 329 |
| age | 49 ± 16 | 44 ± 19 |
| weight | 52 ± 8 | 62 ± 10 |
| height | 153 ± 6 | 165 ± 7 |
| pcfat | 35 ± 5 | 24 ± 6 |
| 1 Mean ± SD | ||
#2.3 Vẽ biểu đồ tán xạ đánh giá mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.
plot(df$height, df$weight,
main = "Biểu đồ tán xạ giữa Chiều cao và Cân nặng",
xlab = "Chiều cao (cm)",
ylab = "Cân nặng (kg)",
pch = 16, col = "blue")
abline(lm(weight ~ height, data = df), col = "red", lwd = 2)
library(ggplot2)
ggplot(df, aes(x = height, y = weight)) +
geom_point(color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "red") +
labs(title = "Biểu đồ tán xạ giữa Chiều cao và Cân nặng",
x = "Chiều cao (cm)",
y = "Cân nặng (kg)") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
#2.4 Tiến hành phân tích tương quan định lượng mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.
#Bước 1: Kiểm tra phân phối chuẩn
#Trước khi chọn loại tương quan, hãy kiểm tra xem 2 biến có phân phối chuẩn không (dùng shapiro.test()):
shapiro.test(df$weight)##
## Shapiro-Wilk normality test
##
## data: df$weight
## W = 0.974, p-value = 5.471e-14 shapiro.test(df$height)##
## Shapiro-Wilk normality test
##
## data: df$height
## W = 0.9827, p-value = 7.179e-11#Bước 2: Tính hệ số tương quan
#Pearson (tương quan tuyến tính):
cor.test(df$weight, df$height, method = "pearson")##
## Pearson's product-moment correlation
##
## data: df$weight and df$height
## t = 25.984, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5602911 0.6326135
## sample estimates:
## cor
## 0.5976667#Ketluan:Phân tích tương quan Pearson cho thấy có mối tương quan tuyến tính dương trung bình giữa cân nặng và chiều cao (r = 0.598, p < 0.001).
#Điều này có nghĩa là chiều cao càng lớn thì cân nặng có xu hướng càng tăng, và mối liên hệ này là có ý nghĩa thống kê.#2.5 Tiến hành phân tích tương quan định lượng mối liên quan giữa chiều cao (height) và tỉ trọng mỡ (pcfat). Nhận xét kết quả.
#Bước 1: Kiểm tra phân phối chuẩn của các biến
#Dùng shapiro.test() để kiểm tra xem dữ liệu có phân phối chuẩn không:
shapiro.test(df$height)##
## Shapiro-Wilk normality test
##
## data: df$height
## W = 0.9827, p-value = 7.179e-11 shapiro.test(df$pcfat)##
## Shapiro-Wilk normality test
##
## data: df$pcfat
## W = 0.98284, p-value = 8.193e-11#Bước 2: Phân tích tương quan
#a. Nếu phân phối chuẩn → Pearson:
cor.test(df$height, df$pcfat, method = "pearson") ##
## Pearson's product-moment correlation
##
## data: df$height and df$pcfat
## t = -19.063, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5219407 -0.4353664
## sample estimates:
## cor
## -0.4798206#Ketluan: Kết quả kiểm định tương quan Pearson cho thấy có mối tương quan tuyến tính âm ở mức trung bình giữa chiều cao và tỉ trọng mỡ cơ thể (r = –0.480, p < 0.001).
#Điều này có nghĩa là những người có chiều cao cao hơn thường có tỉ lệ mỡ cơ thể thấp hơn, và mối quan hệ này có ý nghĩa thống kê. #Việc 3. Hồi qui tuyến tính
#3.1 Đọc dữ liệu "gapminder" vào R từ gói lệnh gapminder. Tạo bộ dữ liệu “vn” gồm dữ liệu của Việt Nam.
library(gapminder)
vn <- subset(gapminder, country == "Vietnam")#3.2 Thực hiện phân tích hồi qui tuyến tính để xác định xem mỗi năm, trong thời gian từ 1952-2007, tuổi thọ của người Việt Nam gia tăng như thế nào?
#Bước 1: Hồi quy tuyến tính
model_year_life <- lm(lifeExp ~ year, data = vn)
summary(model_year_life)##
## Call:
## lm(formula = lifeExp ~ year, data = vn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1884 -0.5840 0.1335 0.7396 1.7873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.272e+03 4.349e+01 -29.25 5.10e-11 ***
## year 6.716e-01 2.197e-02 30.57 3.29e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared: 0.9894, Adjusted R-squared: 0.9884
## F-statistic: 934.5 on 1 and 10 DF, p-value: 3.289e-11#Ve bieu do
plot(vn$year, vn$lifeExp,
main = "Tuổi thọ người Việt Nam qua các năm (1952–2007)",
xlab = "Năm", ylab = "Tuổi thọ",
pch = 16, col = "blue")
abline(model_year_life, col = "red", lwd = 2)
#Ketluan: Kết quả hồi quy tuyến tính cho thấy hệ số của biến year là 0.6716 (p < 0.001), nghĩa là trong giai đoạn từ 1952 đến 2007, tuổi thọ trung bình của người Việt Nam tăng khoảng 0.67 tuổi mỗi năm — và sự gia tăng này là có ý nghĩa thống kê.
#Giá trị R-squared = 0.9894 cho thấy mô hình giải thích tới 98.94% sự biến thiên của tuổi thọ theo thời gian, chứng tỏ mối quan hệ tuyến tính rất mạnh giữa năm và tuổi thọ. #3.3 Kiểm tra các giả định của mô hình hồi qui tuyến tính
#1. Vẽ các biểu đồ chẩn đoán
# Mô hình bạn đã chạy
model_year_life <- lm(lifeExp ~ year, data = vn)
# Vẽ 4 biểu đồ chẩn đoán cơ bản
par(mfrow = c(2, 2))
plot(model_year_life) 
# Giải thích từng biểu đồ:
#Residuals vs Fitted → Kiểm tra tính tuyến tính và phương sai không đổi
#Mong muốn: phân tán ngẫu nhiên quanh đường 0
#Normal Q-Q → Kiểm tra phân phối chuẩn của phần dư
# Mong muốn: các điểm nằm trên đường chéo
#Scale-Location (hoặc Spread-Location) → Kiểm tra phương sai không đổi
#Mong muốn: các điểm phân tán đều (không hình phễu)
#Residuals vs Leverage → Kiểm tra điểm ảnh hưởng quá lớn
#Không nên có điểm nào quá xa (ảnh hưởng mạnh)
# Kiểm tra chuẩn hóa phần dư (normality)
res <- residuals(model_year_life)
shapiro.test(res) # Kiểm định Shapiro-Wilk##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.95663, p-value = 0.7349 # Ketluan: Nếu p-value > 0.05 → phần dư phân phối chuẩn (thỏa giả định)
#Kiểm tra phương sai đồng nhất (homoscedasticity)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Kiểm định Breusch-Pagan
bptest(model_year_life)##
## studentized Breusch-Pagan test
##
## data: model_year_life
## BP = 0.22628, df = 1, p-value = 0.6343#Ketluan: Nếu p-value > 0.05 → không có hiện tượng phương sai thay đổi → thỏa giả định
# Kiểm tra tự tương quan phần dư
dwtest(model_year_life)##
## Durbin-Watson test
##
## data: model_year_life
## DW = 0.94088, p-value = 0.005797
## alternative hypothesis: true autocorrelation is greater than 0#Ketluan: Nếu p-value > 0.05 → không có tự tương quan → thỏa giả định
#Ketluan: kiểm tra giả định không có tự tương quan trong mô hình hồi quy tuyến tín, dùng cac mo hinh ARMIh.#3.4 Viết phương trình đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007. Nhận xét kết quả.
# Hồi quy tuyến tính: Tuổi thọ theo năm
model_year_life <- lm(lifeExp ~ year, data = vn)
# Xem tóm tắt kết quả mô hình
summary(model_year_life)##
## Call:
## lm(formula = lifeExp ~ year, data = vn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1884 -0.5840 0.1335 0.7396 1.7873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.272e+03 4.349e+01 -29.25 5.10e-11 ***
## year 6.716e-01 2.197e-02 30.57 3.29e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared: 0.9894, Adjusted R-squared: 0.9884
## F-statistic: 934.5 on 1 and 10 DF, p-value: 3.289e-11 # Lấy hệ số hồi quy
intercept <- coef(model_year_life)[1]
slope <- coef(model_year_life)[2]
# In phương trình hồi quy
cat("Phương trình hồi quy tuổi thọ theo năm là:\n")## Phương trình hồi quy tuổi thọ theo năm là: cat("lifeExp =", round(intercept, 3), "+", round(slope, 4), "* year\n")## lifeExp = -1271.983 + 0.6716 * year # Viết nhận xét
cat("\n✅ Nhận xét:\n")##
## ✅ Nhận xét: cat("Trong giai đoạn 1952–2007, tuổi thọ trung bình của người Việt Nam tăng khoảng",
round(slope, 2), "tuổi mỗi năm.\n")## Trong giai đoạn 1952–2007, tuổi thọ trung bình của người Việt Nam tăng khoảng 0.67 tuổi mỗi năm. cat("Hệ số R-squared cao cho thấy mối liên hệ rất mạnh giữa năm và tuổi thọ.\n")## Hệ số R-squared cao cho thấy mối liên hệ rất mạnh giữa năm và tuổi thọ. cat("Lưu ý: Kiểm định Durbin-Watson cho thấy có thể có tự tương quan phần dư trong mô hình.\n") ## Lưu ý: Kiểm định Durbin-Watson cho thấy có thể có tự tương quan phần dư trong mô hình.#3.5 Sử dụng ChatGPT viết code để đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007 bằng dữ liệu sau:
# year: 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, 2002, 2007
#lifeExp: 40.4, 42.9, 45.4, 47.8, 50.3, 55.8, 58.8, 62.8, 67.7, 70.7, 73.0, 74.2
# Tạo dữ liệu
year <- c(1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, 2002, 2007)
lifeExp <- c(40.4, 42.9, 45.4, 47.8, 50.3, 55.8, 58.8, 62.8, 67.7, 70.7, 73.0, 74.2)
# Kết hợp thành data frame
vn_data <- data.frame(year, lifeExp)
# Hồi quy tuyến tính
model <- lm(lifeExp ~ year, data = vn_data)
# Tóm tắt kết quả
summary(model)##
## Call:
## lm(formula = lifeExp ~ year, data = vn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1494 -0.5944 0.1387 0.7324 1.8268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.271e+03 4.377e+01 -29.04 5.47e-11 ***
## year 6.712e-01 2.211e-02 30.35 3.53e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.322 on 10 degrees of freedom
## Multiple R-squared: 0.9893, Adjusted R-squared: 0.9882
## F-statistic: 921.4 on 1 and 10 DF, p-value: 3.527e-11 # Trích xuất hệ số
intercept <- coef(model)[1]
slope <- coef(model)[2]
# In phương trình hồi quy
cat("📌 Phương trình hồi quy tuổi thọ theo năm:\n")## 📌 Phương trình hồi quy tuổi thọ theo năm: cat("lifeExp =", round(intercept, 3), "+", round(slope, 4), "* year\n")## lifeExp = -1271.135 + 0.6712 * year # Nhận xét đơn giản
cat("\n✅ Nhận xét:\n")##
## ✅ Nhận xét: cat("Tuổi thọ trung bình của người Việt Nam tăng khoảng", round(slope, 2),
"tuổi mỗi năm trong giai đoạn 1952–2007.\n")## Tuổi thọ trung bình của người Việt Nam tăng khoảng 0.67 tuổi mỗi năm trong giai đoạn 1952–2007. plot(vn_data$year, vn_data$lifeExp,
type = "p",
pch = 16,
col = "blue",
main = "Gia tăng tuổi thọ người Việt Nam (1952–2007)",
xlab = "Năm", ylab = "Tuổi thọ (năm)")
abline(model, col = "red", lwd = 2)
#Việc 4. Hồi qui tuyến tính đa biến
#4.1 Giả sử bạn có dữ liệu về ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1) và thời gian lưu (X2) như sau:
#Y: 12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6
# X1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
# X2: 7, 4, 4, 6, 4, 2, 1, 1, 1, 0
# Nhập dữ liệu vào R và đặt tên tập dữ liệu là df.
# Tạo các biến
Y <- c(12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6)
X1 <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
X2 <- c(7, 4, 4, 6, 4, 2, 1, 1, 1, 0)
# Tạo data frame
df <- data.frame(Y, X1, X2)
# Xem thử dữ liệu
print(df)## Y X1 X2
## 1 12.1 0 7
## 2 11.9 1 4
## 3 10.2 2 4
## 4 8.0 3 6
## 5 7.7 4 4
## 6 5.3 5 2
## 7 7.9 6 1
## 8 7.8 7 1
## 9 5.5 8 1
## 10 2.6 9 0#42 Đánh giá mối liên quan giữa ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1).
# Hồi quy tuyến tính đơn biến: Y theo X1
model_X1 <- lm(Y ~ X1, data = df)
# Tóm tắt kết quả mô hình
summary(model_X1)##
## Call:
## lm(formula = Y ~ X1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1606 -1.0735 0.1742 0.8621 2.0970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.8545 0.8283 14.312 5.54e-07 ***
## X1 -0.8788 0.1552 -5.664 0.000474 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.409 on 8 degrees of freedom
## Multiple R-squared: 0.8004, Adjusted R-squared: 0.7755
## F-statistic: 32.08 on 1 and 8 DF, p-value: 0.0004737 # Biểu đồ tán xạ + đường hồi quy
plot(df$X1, df$Y,
main = "Mối liên quan giữa Y và X1",
xlab = "Hoạt tính bạc (X1)",
ylab = "Ethylene Oxide (Y)",
pch = 16, col = "blue")
abline(model_X1, col = "red", lwd = 2)
# Trích xuất hệ số
intercept <- coef(model_X1)[1]
slope <- coef(model_X1)[2]
r2 <- summary(model_X1)$r.squared
p_val <- summary(model_X1)$coefficients[2,4]
# In phương trình và nhận xét
cat("📌 Phương trình hồi quy:\n")## 📌 Phương trình hồi quy: cat("Y =", round(intercept, 4), "+", round(slope, 4), "* X1\n\n")## Y = 11.8545 + -0.8788 * X1 cat("📝 Nhận xét:\n")## 📝 Nhận xét: cat("- Khi hoạt tính bạc (X1) tăng 1 đơn vị, ethylene oxide (Y) trung bình giảm khoảng",
abs(round(slope, 2)), "đơn vị.\n")## - Khi hoạt tính bạc (X1) tăng 1 đơn vị, ethylene oxide (Y) trung bình giảm khoảng 0.88 đơn vị. cat("- P-value =", signif(p_val, 3), "→", ifelse(p_val < 0.05,
"Ảnh hưởng của X1 có ý nghĩa thống kê.\n",
"Không có ý nghĩa thống kê.\n"))## - P-value = 0.000474 → Ảnh hưởng của X1 có ý nghĩa thống kê. cat("- R-squared =", round(r2 * 100, 1), "% → X1 giải thích được",
round(r2 * 100, 1), "% biến thiên của Y.\n") ## - R-squared = 80 % → X1 giải thích được 80 % biến thiên của Y.#4.3 Đánh giá mối liên quan giữa ethylene oxide (Y) và thời gian lưu (X2).
# Hồi quy đơn biến: Y theo X2
model_X2 <- lm(Y ~ X2, data = df)
# Tóm tắt mô hình
summary(model_X2)##
## Call:
## lm(formula = Y ~ X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.702 -1.533 -0.034 1.667 3.066
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0980 1.1222 4.543 0.00189 **
## X2 0.9340 0.2999 3.114 0.01436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.121 on 8 degrees of freedom
## Multiple R-squared: 0.548, Adjusted R-squared: 0.4915
## F-statistic: 9.698 on 1 and 8 DF, p-value: 0.01436 # Trích xuất hệ số và thông tin
intercept <- coef(model_X2)[1]
slope <- coef(model_X2)[2]
r2 <- summary(model_X2)$r.squared
p_val <- summary(model_X2)$coefficients[2, 4]
# In phương trình và diễn giải
cat("📌 Phương trình hồi quy:\n")## 📌 Phương trình hồi quy: cat("Y =", round(intercept, 4), "+", round(slope, 4), "* X2\n\n")## Y = 5.098 + 0.934 * X2 cat("📝 Nhận xét:\n")## 📝 Nhận xét: cat("- Khi thời gian lưu (X2) tăng 1 đơn vị, lượng ethylene oxide (Y) trung bình thay đổi khoảng",
round(slope, 2), "đơn vị.\n")## - Khi thời gian lưu (X2) tăng 1 đơn vị, lượng ethylene oxide (Y) trung bình thay đổi khoảng 0.93 đơn vị. cat("- P-value =", signif(p_val, 3), "→", ifelse(p_val < 0.05,
"Ảnh hưởng của X2 có ý nghĩa thống kê.\n",
"Không có ý nghĩa thống kê.\n"))## - P-value = 0.0144 → Ảnh hưởng của X2 có ý nghĩa thống kê. cat("- R-squared =", round(r2 * 100, 1), "% → X2 giải thích được",
round(r2 * 100, 1), "% biến thiên của Y.\n")## - R-squared = 54.8 % → X2 giải thích được 54.8 % biến thiên của Y.#4.4 Bạn muốn đánh giá mối liên quan độc lập giữa thời gian lưu (X2) và ethylene oxide (Y). Nhận xét kết quả.
# Mô hình hồi quy đa biến
model_multi <- lm(Y ~ X1 + X2, data = df)
# Tóm tắt kết quả
summary(model_multi)##
## Call:
## lm(formula = Y ~ X1 + X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.46078 -0.33384 0.00026 0.81856 1.98476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7076 2.9785 4.938 0.00168 **
## X1 -1.2042 0.3614 -3.332 0.01255 *
## X2 -0.4629 0.4642 -0.997 0.35187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.41 on 7 degrees of freedom
## Multiple R-squared: 0.8252, Adjusted R-squared: 0.7753
## F-statistic: 16.53 on 2 and 7 DF, p-value: 0.002232 cat("Kết quả hồi quy đa biến:\n")## Kết quả hồi quy đa biến: cat("Phương trình hồi quy: Y =", round(14.7076, 3),
"-", round(1.2042, 3), "* X1",
"-", round(0.4629, 3), "* X2\n\n")## Phương trình hồi quy: Y = 14.708 - 1.204 * X1 - 0.463 * X2 cat("Hệ số hồi quy:\n")## Hệ số hồi quy: cat("- Hoạt tính bạc (X1): Estimate = -1.2042, p-value = 0.01255\n")## - Hoạt tính bạc (X1): Estimate = -1.2042, p-value = 0.01255 cat(" -> Có ảnh hưởng nghịch chiều và có ý nghĩa thống kê đến ethylene oxide (Y).\n\n")## -> Có ảnh hưởng nghịch chiều và có ý nghĩa thống kê đến ethylene oxide (Y). cat("- Thời gian lưu (X2): Estimate = -0.4629, p-value = 0.35187\n")## - Thời gian lưu (X2): Estimate = -0.4629, p-value = 0.35187 cat(" -> Không có ảnh hưởng độc lập đáng kể đến Y khi đã điều chỉnh X1.\n\n")## -> Không có ảnh hưởng độc lập đáng kể đến Y khi đã điều chỉnh X1. cat("Độ phù hợp mô hình:\n")## Độ phù hợp mô hình: cat("- R-squared =", 0.8252, "\n")## - R-squared = 0.8252 cat("- P-value của F-test =", 0.002232, "\n")## - P-value của F-test = 0.002232 cat("=> Mô hình hồi quy tổng thể có ý nghĩa thống kê và giải thích được 82.5% biến thiên của Y.\n")## => Mô hình hồi quy tổng thể có ý nghĩa thống kê và giải thích được 82.5% biến thiên của Y.