Bai tap 33
Dương
2025-10-24
Ngày 3: Hồi qui tuyến tính (23/10/25)
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, unionlibrary(ggplot2)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeViệc 1. Phân tích phương sai
1.1. Nhập dữ liệu
Tạo vector chứa tất cả cân nặng
weight <- c(8, 9, 11, 4, 7, 8, 5, # Nhóm A (n=7)
7, 17, 10, 14, 12, 24, 11, 22, # Nhóm B (n=8)
28, 21, 26, 11, 24, 19, # Nhóm C (n=6)
26, 16, 13, 12, 9, 10, 11, 17, 15) # Nhóm D (n=9)Tạo vector yếu tố (factor) cho các nhóm tương ứng
group <- factor(c(rep("A", 7),
rep("B", 8),
rep("C", 6),
rep("D", 9)))Kết hợp thành một data frame tên là ‘data’
data <- data.frame(group = group, weight = weight)In ra vài dòng đầu và cấu trúc dữ liệu để kiểm tra
print(head(data))## group weight
## 1 A 8
## 2 A 9
## 3 A 11
## 4 A 4
## 5 A 7
## 6 A 8str(data)## 'data.frame': 30 obs. of 2 variables:
## $ group : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 1 2 2 2 ...
## $ weight: num 8 9 11 4 7 8 5 7 17 10 ...1.2. Mô tả cân nặng giữa 4 nhóm
1.2.1. Tính toán thống kê mô tả cho mỗi nhóm
summary_stats <- data %>%
group_by(group) %>%
summarise(
so_luong = n(),
trung_binh = mean(weight, na.rm = TRUE),
do_lech_chuan = sd(weight, na.rm = TRUE),
trung_vi = median(weight, na.rm = TRUE),
gia_tri_nho_nhat = min(weight, na.rm = TRUE),
gia_tri_lon_nhat = max(weight, na.rm = TRUE)
)
print("Thống kê mô tả theo nhóm:")## [1] "Thống kê mô tả theo nhóm:"print(summary_stats)## # A tibble: 4 × 7
## group so_luong trung_binh do_lech_chuan trung_vi gia_tri_nho_nhat
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 A 7 7.43 2.37 8 4
## 2 B 8 14.6 5.95 13 7
## 3 C 6 21.5 6.09 22.5 11
## 4 D 9 14.3 5.15 13 9
## # ℹ 1 more variable: gia_tri_lon_nhat <dbl>1.2.2. Trực quan hóa dữ liệu bằng Boxplot (sử dụng ggplot2)
ggplot(data, aes(x = group, y = weight, fill = group)) +
geom_boxplot() +
geom_jitter(shape = 16, position = position_jitter(0.2), alpha = 0.6) +
labs(title = "Phân phối cân nặng theo 4 nhóm",
x = "Nhóm (A, B, C, D)",
y = "Cân nặng") +
theme_minimal() +
scale_fill_brewer(palette = "Pastel1")
1.3. Phân tích sự khác biệt về cân nặng (ANOVA)
1.3.1. Xây dựng mô hình ANOVA
model <- aov(weight ~ group, data = data)1.3.2. Kiểm tra các giả định ANOVA
=========================================================
Giả định 1: Các phần dư (residuals) phải phân phối chuẩn
=========================================================
Lấy phần dư từ mô hình
model_residuals <- residuals(model)Sử dụng Shapiro-Wilk Test
shapiro_test <- shapiro.test(model_residuals)
print("Kiểm định Shapiro-Wilk cho phần dư:")## [1] "Kiểm định Shapiro-Wilk cho phần dư:"print(shapiro_test)##
## Shapiro-Wilk normality test
##
## data: model_residuals
## W = 0.98028, p-value = 0.8329===================================================================
Giả định 2: Tính đồng nhất của phương sai (Homogeneity of variance
===================================================================
Các nhóm phải có phương sai tương đương nhau Sử dụng Levene’s Test (từ thư viện ‘car’)
levene_test <- leveneTest(weight ~ group, data = data)
print("Kiểm định Levene cho tính đồng nhất phương sai:")## [1] "Kiểm định Levene cho tính đồng nhất phương sai:"print(levene_test)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.112 0.3622
## 261.3.3. Xem kết quả phân tích ANOVA
print("Kết quả ANOVA:")## [1] "Kết quả ANOVA:"print(summary(model))## Df Sum Sq Mean Sq F value Pr(>F)
## group 3 642.3 214.09 8.197 0.000528 ***
## Residuals 26 679.1 26.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.4. Phân tích hậu định (Post-hoc analysis)
Thực hiện phân tích post-hoc Tukey HSD (Chỉ chạy khi p-value của ANOVA < 0.05)
tukey_results <- TukeyHSD(model)
print("Kết quả phân tích Post-hoc (Tukey HSD):")## [1] "Kết quả phân tích Post-hoc (Tukey HSD):"print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weight ~ group, data = data)
##
## $group
## diff lwr upr p adj
## B-A 7.1964286 -0.05969765 14.4525548 0.0525014
## C-A 14.0714286 6.27132726 21.8715299 0.0002134
## D-A 6.9047619 -0.16073856 13.9702624 0.0571911
## C-B 6.8750000 -0.69675602 14.4467560 0.0850381
## D-B -0.2916667 -7.10424368 6.5209103 0.9994049
## D-C -7.1666667 -14.55594392 0.2226106 0.0597131CHUYỂN KẾT QUẢ SANG DATA FRAME Trích xuất bảng kết quả cho ‘group’
Thêm cột “comparison” (các cặp so sánh) từ tên hàng (rownames)
tukey_df <- as.data.frame(tukey_results$group)
tukey_df$comparison <- rownames(tukey_df)TẠO CỘT MỚI ĐỂ ĐÁNH DẤU SỰ KHÁC BIỆT Chúng ta tạo cột ‘significant’ Nếu ‘p adj’ < 0.05, gán “Co” (Có khác biệt), ngược lại gán “Khong”
tukey_df$significant <- ifelse(tukey_df$"p adj" < 0.05, "Co", "Khong")In ra thử data frame mới để xem
print(tukey_df)## diff lwr upr p adj comparison significant
## B-A 7.1964286 -0.05969765 14.4525548 0.0525013886 B-A Khong
## C-A 14.0714286 6.27132726 21.8715299 0.0002133948 C-A Co
## D-A 6.9047619 -0.16073856 13.9702624 0.0571911198 D-A Khong
## C-B 6.8750000 -0.69675602 14.4467560 0.0850381279 C-B Khong
## D-B -0.2916667 -7.10424368 6.5209103 0.9994049238 D-B Khong
## D-C -7.1666667 -14.55594392 0.2226106 0.0597130963 D-C KhongVẼ BIỂU ĐỒ BẰNG GGPLOT2
library(ggplot2)
ggplot(tukey_df, aes(x = comparison, y = diff, color = significant)) +
# Vẽ điểm chấm cho giá trị trung bình (diff)
geom_point(size = 3) +
# Vẽ khoảng tin cậy 95% (từ lwr đến upr)
geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2, size = 1) +
# Thêm một đường ngang tại số 0 (ngưỡng quan trọng)
# Nếu khoảng tin cậy cắt qua đường này -> không có ý nghĩa
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
# Tùy chỉnh màu sắc
scale_color_manual(values = c("Co" = "red", "Khong" = "blue"), # Đỏ = Có, Xanh = Không
name = "Khac biet (p < 0.05)") +
# Lật ngược trục để các cặp so sánh nằm ngang cho dễ đọc
coord_flip() +
# Thêm tiêu đề và nhãn
labs(title = "Ket qua phan tich Post-hoc (Tukey HSD)",
x = "Cac cap so sanh",
y = "Su khac biet ve trung binh (95% CI)") +
# Dùng theme gọn gàng
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:car':
##
## logit## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(ggpubr)Việc 2. Phân tích tương quan
2.1. Đọc dữ liệu
df <- read.csv("C:\\Users\\duong\\OneDrive\\Máy tính\\Phân tích dữ liệu và công bố quốc tế\\211025\\Demo.csv")
head(df)## X id age gender weight height pcfat
## 1 1 1 53 F 49 150 37.3
## 2 2 2 65 M 52 165 16.8
## 3 3 3 64 F 57 157 34.0
## 4 4 4 56 F 53 156 33.8
## 5 5 5 54 M 51 160 14.8
## 6 6 6 52 F 47 153 32.2str(df)## 'data.frame': 1217 obs. of 7 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ age : int 53 65 64 56 54 52 66 50 61 58 ...
## $ gender: chr "F" "M" "F" "F" ...
## $ weight: int 49 52 57 53 51 47 58 65 54 60 ...
## $ height: int 150 165 157 156 160 153 155 167 165 158 ...
## $ pcfat : num 37.3 16.8 34 33.8 14.8 32.2 35.3 28 21.1 36.6 ...2.2. Mô tả đặc điểm cân nặng (weight) và chiều cao (height)
Chọn các cột ‘weight’ và ‘height’
mota_weight_height <- describe(df[c("weight", "height")])
print("Mô tả thống kê cho Cân nặng và Chiều cao:")## [1] "Mô tả thống kê cho Cân nặng và Chiều cao:"print(mota_weight_height)## vars n mean sd median trimmed mad min max range skew kurtosis
## weight 1 1217 55.14 9.40 54 54.62 8.90 34 95 61 0.67 0.86
## height 2 1217 156.72 7.98 155 156.36 7.41 136 185 49 0.47 0.16
## se
## weight 0.27
## height 0.232.3. Vẽ biểu đồ tán xạ (Weight vs Height)
ggscatter(df, x = "height", y = "weight",
color = "gender", # Yêu cầu R phân biệt màu sắc theo cột 'gender'
palette = "jco", # Chọn một bảng màu (ví dụ: "jco", "npg", "aaas")
add = "reg.line", # Thêm đường hồi quy (sẽ có 2 đường cho 2 nhóm)
conf.int = TRUE, # Thêm khoảng tin cậy cho mỗi đường
cor.coef = TRUE, # Hiển thị hệ số tương quan (R) cho từng nhóm
cor.method = "pearson",
xlab = "Chiều cao (cm)",
ylab = "Cân nặng (kg)",
title = "Mối liên quan giữa Cân nặng và Chiều cao (theo giới tính)")
Nhận xét kết quả :
Phía trên biểu đồ, thấy hai hệ số tương quan (R) và giá trị p-value, một cho nhóm Nam (M) và một cho nhóm Nữ (F). Điều này cho phép so sánh: “Mối liên quan giữa chiều cao và cân nặng có giống nhau hay khác nhau giữa nam và nữ không?” Ví dụ, mối tương quan này ở Nam (r=…) mạnh hơn hoặc yếu hơn so với ở Nữ (r=..$).
2.4. Phân tích tương quan (weight và height)
tuong_quan_WH <- cor.test(df$weight, df$height, method = "pearson")
print("Kết quả tương quan giữa Cân nặng và Chiều cao:")## [1] "Kết quả tương quan giữa Cân nặng và Chiều cao:"print(tuong_quan_WH)##
## 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.5976667Nhận xét kết quả :
Hệ số tương quan (r): r ~ 0.598 Đây là một mối tương quan dương (cùng chiều) và ở mức độ trung bình (khá mạnh Giá trị p-value: p ~ 9.69e-119 (một con số cực kỳ nhỏ, gần như bằng 0).
Kết luận: Vì p < 0.05 (thực tế là p << 0.05), chúng ta kết luận rằng có mối tương quan dương, mạnh vừa phải, và có ý nghĩa thống kê rất cao giữa cân nặng và chiều cao. Nói cách khác, mối liên hệ “cao hơn thì nặng hơn” trong dữ liệu này là thật sự có ý nghĩa, không phải do ngẫu nhiên.
2.5. Phân tích tương quan (height và pcfat)
tuong_quan_HP <- cor.test(df$height, df$pcfat, method = "pearson")
print("Kết quả tương quan giữa Chiều cao và Tỉ trọng mỡ:")## [1] "Kết quả tương quan giữa Chiều cao và Tỉ trọng mỡ:"print(tuong_quan_HP)##
## 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.4798206Nhận xét kết quả :
Hệ số tương quan (r): r ~ -0.480. Đây là một mối tương quan âm (ngược chiều) và ở mức độ trung bình. Giá trị p-value: p ~ 4.39e-71 (cũng là một con số cực kỳ nhỏ).
Kết luận: Vì p < 0.05, chúng ta kết luận rằng có mối tương quan âm, ở mức độ trung bình, và có ý nghĩa thống kê rất cao giữa chiều cao và tỉ trọng mỡ.
Điều này có nghĩa là, trong bộ dữ liệu này, những người cao hơn có xu hướng có tỉ trọng mỡ cơ thể (pcfat) thấp hơn.
Việc 3. Hồi quy tuyến tính đơn biến (Simple Linear Regression)
library(gapminder)
library(dplyr)
library(ggplot2)
library(lessR) ##
## lessR 4.4.5 feedback: gerbing@pdx.edu
## --------------------------------------------------------------
## > d <- Read("") Read data file, many formats available, e.g., Excel
## d is default data frame, data= in analysis routines optional
##
## Many examples of reading, writing, and manipulating data,
## graphics, testing means and proportions, regression, factor analysis,
## customization, forecasting, and aggregation from pivot tables
## Enter: browseVignettes("lessR")
##
## View lessR updates, now including time series forecasting
## Enter: news(package="lessR")
##
## Interactive data analysis
## Enter: interact()##
## Attaching package: 'lessR'## The following objects are masked from 'package:psych':
##
## reflect, rescale, scree, skew## The following objects are masked from 'package:car':
##
## bc, recode, sp## The following objects are masked from 'package:dplyr':
##
## order_by, recode, renamelibrary(ggfortify)
library(broom)3.1. Đọc dữ liệu và tạo bộ dữ liệu “vn”
df <- read.csv("C:\\Users\\duong\\OneDrive\\Máy tính\\Phân tích dữ liệu và công bố quốc tế\\211025\\Demo.csv")
vn <- gapminder %>%
filter(country == "Vietnam")
# In ra để kiểm tra
print("Dữ liệu của Việt Nam (1952-2007):")## [1] "Dữ liệu của Việt Nam (1952-2007):"print(vn)## # A tibble: 12 × 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Vietnam Asia 1952 40.4 26246839 605.
## 2 Vietnam Asia 1957 42.9 28998543 676.
## 3 Vietnam Asia 1962 45.4 33796140 772.
## 4 Vietnam Asia 1967 47.8 39463910 637.
## 5 Vietnam Asia 1972 50.3 44655014 700.
## 6 Vietnam Asia 1977 55.8 50533506 714.
## 7 Vietnam Asia 1982 58.8 56142181 707.
## 8 Vietnam Asia 1987 62.8 62826491 821.
## 9 Vietnam Asia 1992 67.7 69940728 989.
## 10 Vietnam Asia 1997 70.7 76048996 1386.
## 11 Vietnam Asia 2002 73.0 80908147 1764.
## 12 Vietnam Asia 2007 74.2 85262356 2442.Plot(lifeExp, year, data = vn, xlab = "Life expectancy (years)", ylab = "Year")
##
## >>> Suggestions or enter: style(suggest=FALSE)
## Plot(lifeExp, year, enhance=TRUE) # many options
## Plot(lifeExp, year, color="red") # exterior edge color of points
## Plot(lifeExp, year, fit="lm", fit_se=c(.90,.99)) # fit line, stnd errors
## Plot(lifeExp, year, MD_cut=6) # Mahalanobis distance from center > 6 is an outlier
##
##
## >>> Pearson's product-moment correlation
##
## Number of paired values with neither missing, n = 12
## Sample Correlation of lifeExp and year: r = 0.995
##
## Hypothesis Test of 0 Correlation: t = 30.569, df = 10, p-value = 0.000
## 95% Confidence Interval for Correlation: 0.981 to 0.999
## Kết quả (Diễn giải): Bảng dữ liệu (data.frame) tên vn chứa 12 dòng, tương ứng với 12 mốc thời gian (5 năm một lần, từ 1952 đến 2007) cho Việt Nam.
3.2. Diễn giải kết quả mô hình
3.2.1. Trực quan hóa
print("Đang vẽ biểu đồ (ggplot)...")## [1] "Đang vẽ biểu đồ (ggplot)..."ggplot(vn, aes(x = year, y = lifeExp)) +
# 1. Vẽ các điểm (biểu đồ tán xạ) - Tô màu xanh
geom_point(color = "blue", size = 3) +
# 2. Thêm đường hồi quy tuyến tính (method = "lm") màu đỏ
# VÀ khoảng tin cậy 95% (se = TRUE) (vùng màu xám)
# (Bước này tương đương với add = "reg.line" VÀ conf.int = TRUE)
geom_smooth(method = "lm", color = "red", fill = "grey80", se = TRUE) +
labs(title = "Mối liên quan giữa Tuổi thọ và Năm tại Việt Nam",
x = "Năm",
ylab = "Tuổi thọ") +
theme_minimal()## Ignoring unknown labels:
## • ylab : "Tuổi thọ"
## `geom_smooth()` using formula = 'y ~ x'
#### 3.2.2. (Tùy chọn) Tính toán hệ số tương quan R (vì ggplot không tự
hiển thị)
cor_test_vn <- cor.test(vn$year, vn$lifeExp)
print(paste0("Hệ số tương quan R: ", round(cor_test_vn$estimate, 3)))## [1] "Hệ số tương quan R: 0.995"Nhận xét biểu đồ: Biểu đồ (tạo bởi ggplot) cho thấy các điểm màu xanh bám rất sát đường hồi quy màu đỏ, cho thấy mối tương quan tuyến tính dương (đi lên) rất mạnh.
Kết quả cor.test ( \(R \approx 0.996\)) cũng xác nhận điều này.
3.2.3. Xây dựng mô hình hồi quy
model_vn <- lm(lifeExp ~ year, data = vn)3.2.4. Xem kết quả tóm tắt chi tiết của mô hình
print("Kết quả tóm tắt mô hình (lifeExp ~ year):")## [1] "Kết quả tóm tắt mô hình (lifeExp ~ year):"summary(model_vn)##
## 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) -1271.98315 43.49240 -29.25 0.0000000000510 ***
## year 0.67162 0.02197 30.57 0.0000000000329 ***
## ---
## 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: 0.000000000032893.3. Kiểm tra các giả định của mô hình
Lấy các giá trị cần thiết từ mô hình
vn$residuals <- residuals(model_vn)
vn$fitted <- fitted.values(model_vn)3.3.1. Vẽ 4 biểu đồ chẩn đoán đầu tiên
(Residuals vs Fitted, Normal Q-Q, Scale-Location, Residuals vs Leverage)
print("Đang vẽ 4 biểu đồ (autoplot)...")## [1] "Đang vẽ 4 biểu đồ (autoplot)..."autoplot(model_vn)
Nhận xét: Bốn biểu đồ này xuất hiện cùng lúc, cho bạn cái nhìn tổng quan:
Residuals vs Fitted: Phải ngẫu nhiên (dữ liệu vn bị cong).
Normal Q-Q: Các điểm phải nằm trên đường thẳng (dữ liệu vn bị lệch ở 2 đầu).
Scale-Location: Đường đỏ phải nằm ngang (dữ liệu vn hơi dốc lên).
Residuals vs Leverage: Kiểm tra các điểm ngoại lai (outliers) có ảnh hưởng lớn (dữ liệu vn có điểm 1 và 12 là các điểm đòn bẩy cao).
— Giả định 1: Phần dư phân phối chuẩn (Kiểm định Shapiro-Wilk) —
print("Kiểm định Shapiro-Wilk (tính chuẩn của phần dư):")## [1] "Kiểm định Shapiro-Wilk (tính chuẩn của phần dư):"shapiro.test(vn$residuals)##
## Shapiro-Wilk normality test
##
## data: vn$residuals
## W = 0.95663, p-value = 0.7349Nhận xét (Shapiro-Wilk): Nếu \(p < 0.05\), giả định phân phối chuẩn bị vi phạm. #### — Giả định 2: Phần dư phân phối chuẩn (Biểu đồ Q-Q có màu) —
print("Đang vẽ biểu đồ Q-Q (biểu đồ màu)...")## [1] "Đang vẽ biểu đồ Q-Q (biểu đồ màu)..."ggplot(vn, aes(sample = residuals)) +
# Tô màu xanh cho các điểm
stat_qq(color = "blue", size = 2) +
# Thêm đường thẳng tham chiếu màu đỏ
stat_qq_line(color = "red", linetype = "dashed", size = 1) +
labs(title = "Biểu đồ Q-Q của Phần dư",
x = "Giá trị lý thuyết (Theoretical Quantiles)",
y = "Giá trị mẫu (Sample Quantiles)") +
theme_minimal()
Nhận xét (Biểu đồ Q-Q):
Các điểm màu xanh (đặc biệt là ở hai đầu) lệch ra khỏi đường nét đứt màu đỏ, cho thấy giả định chuẩn bị vi phạm.
— Giả định 3: Tính đồng nhất phương sai (Biểu đồ Residuals vs. Fitted) —
print("Đang vẽ biểu đồ Phần dư vs. Giá trị dự đoán (biểu đồ màu)...")## [1] "Đang vẽ biểu đồ Phần dư vs. Giá trị dự đoán (biểu đồ màu)..."ggplot(vn, aes(x = fitted, y = residuals)) +
# Tô màu xanh cho các điểm
geom_point(color = "blue", size = 3) +
# Thêm đường tham chiếu y=0 màu đỏ
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
labs(title = "Biểu đồ Phần dư vs. Giá trị dự đoán",
x = "Tuổi thọ dự đoán (Fitted)",
y = "Phần dư (Residuals)") +
theme_minimal()
Nhận xét (Residuals vs. Fitted):
Các điểm màu xanh có xu hướng tạo thành hình vòng cung (chữ U ngược), cho thấy có thể mối quan hệ không phải là tuyến tính hoàn hảo (vi phạm giả định).
— Giả định 4: Tính độc lập của phần dư (Kiểm định Durbin-Watson) —
print("Kiểm định Durbin-Watson (tính độc lập của phần dư):")## [1] "Kiểm định Durbin-Watson (tính độc lập của phần dư):"durbinWatsonTest(model_vn)## lag Autocorrelation D-W Statistic p-value
## 1 0.3888843 0.9408848 0.012
## Alternative hypothesis: rho != 0Nhận xét (Durbin-Watson):
Nếu \(p < 0.05\), giả định độc lập bị vi phạm (có hiện tượng tự tương quan).
3.3.2. Trích xuất các chỉ số chẩn đoán bằng broom::augment()
model_vn.diag.metrics <- augment(model_vn)In ra xem thử
print(head(model_vn.diag.metrics))## # A tibble: 6 × 8
## lifeExp year .fitted .resid .hat .sigma .cooksd .std.resid
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 40.4 1952 39.0 1.40 0.295 1.27 0.338 1.27
## 2 42.9 1957 42.4 0.519 0.225 1.37 0.0292 0.449
## 3 45.4 1962 45.7 -0.363 0.169 1.38 0.00936 -0.303
## 4 47.8 1967 49.1 -1.25 0.127 1.31 0.0750 -1.02
## 5 50.3 1972 52.4 -2.19 0.0991 1.15 0.169 -1.76
## 6 55.8 1977 55.8 -0.0365 0.0851 1.38 0.0000392 -0.0290Biểu đồ Phần dư (Residuals)
(Cột .resid trong data mới)
ggplot(model_vn.diag.metrics, aes(x = year, y = .resid)) +
geom_point(color = "blue") +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Phần dư theo thời gian")
Biểu đồ Khoảng cách Cook (Cook’s distance)
(Kiểm tra điểm ảnh hưởng lớn)
ggplot(model_vn.diag.metrics, aes(x = year, y = .cooksd)) +
geom_bar(stat = "identity", fill = "blue") +
labs(title = "Khoảng cách Cook (Cook's distance)")
Biểu đồ DFBeta (cho ‘year’)
(Sự thay đổi hệ số ‘year’ nếu bỏ 1 điểm)
model_vn.diag.metrics <- augment(model_vn)
dfbeta_data <- as.data.frame(dfbetas(model_vn))
print("Tên cột DFBeta gốc:")## [1] "Tên cột DFBeta gốc:"head(dfbeta_data)## (Intercept) year
## 1 0.725102880 -0.721183578
## 2 0.185005464 -0.183783184
## 3 -0.093612069 0.092817393
## 4 -0.230332923 0.227601341
## 5 -0.269814802 0.264517049
## 6 -0.001275811 0.001203413names(dfbeta_data) <- c("dfb_intercept", "dfb_year")
print("Tên cột DFBeta đã đổi:")## [1] "Tên cột DFBeta đã đổi:"head(dfbeta_data)## dfb_intercept dfb_year
## 1 0.725102880 -0.721183578
## 2 0.185005464 -0.183783184
## 3 -0.093612069 0.092817393
## 4 -0.230332923 0.227601341
## 5 -0.269814802 0.264517049
## 6 -0.001275811 0.001203413model_vn.diag.metrics <- bind_cols(model_vn.diag.metrics, dfbeta_data)ggplot(model_vn.diag.metrics, aes(x = year, y = dfb_year)) +
geom_bar(stat = "identity", fill = "blue") +
labs(title = "DFBeta (ảnh hưởng lên hệ số 'year')",
y = "DFBeta (year)") +
theme_minimal()
Biểu đồ DFFits (ảnh hưởng lên giá trị dự đoán)
(Sự thay đổi giá trị dự đoán nếu bỏ 1 điểm)
dffits_data <- dffits(model_vn)
model_vn.diag.metrics$dffits_manual <- dffits_dataggplot(model_vn.diag.metrics, aes(x = year, y = dffits_manual)) +
geom_bar(stat = "identity", fill = "blue") +
labs(title = "DFFits (ảnh hưởng lên giá trị dự đoán)",
y = "DFFits") +
theme_minimal()
3.4. Viết phương trình và Nhận xét
summary(model_vn)##
## 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) -1271.98315 43.49240 -29.25 0.0000000000510 ***
## year 0.67162 0.02197 30.57 0.0000000000329 ***
## ---
## 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: 0.00000000003289Chúng ta cần 2 con số trong cột Estimate: Hệ số chặn (Intercept): -1.213e+03 (Đây là cách R viết tắt, nó có nghĩa là \(-1213.0\)) Hệ số góc (của biến year): 6.489e-01 (Có nghĩa là \(0.6489\)) Phương trình hồi quy tuyến tính có dạng chung là: \(Y = \beta_0 + \beta_1 \times X\) Trong đó: \(Y\) là biến phụ thuộc (biến bạn muốn dự đoán):Tuổi thọ (lifeExp) \(X\) là biến độc lập (biến bạn dùng để dự đoán): Năm (year) \(\beta_0\) là hệ số chặn (giá trị (Intercept)): \(-1213.0\)\(\beta_1\) là hệ số góc (giá trị year): \(0.6489\)
\[\text{Tuổi thọ dự đoán} = -1213.0 +
(0.6489 \times \text{Năm}) \] Nhận xét kết quả (Diễn
giải) cần bình luận về 3 yếu tố chính từ bảng
summary:
- Hệ số góc (Estimate của
year): Hệ số củayearlà \(0.6489\) (một số dương).
Ý nghĩa: Điều này có nghĩa là, trung bình, cứ mỗi năm trôi qua (trong giai đoạn 1952-2007), tuổi thọ của người Việt Nam gia tăng 0.649 năm. Đây là phát hiện quan trọng nhất của mô hình.
- Giá trị P (p-value của
year):
Nhìn vào cột Pr(>|t|) của hàng year, giá
trị là < 2e-16 (một con số cực kỳ nhỏ, gần như bằng
0).
Ý nghĩa: Vì \(p <
0.05\), chúng ta kết luận rằng yếu tố year
(thời gian) có ảnh hưởng rất có ý nghĩa thống kê đến
lifeExp (tuổi thọ). Mối liên quan (sự gia tăng) này là
thật, không phải do ngẫu nhiên.
- Hệ số \(R^2\) (R-squared):
Nhìn vào dòng Multiple R-squared: 0.9928.
Ý nghĩa: \(R^2 =
99.28\%\). Con số này cho biết rằng 99.28% sự
biến thiên (sự gia tăng) của tuổi thọ (lifeExp) được giải
thích bởi sự thay đổi của thời gian (year).
Đây là một mô hình mô tả dữ liệu rất sát, cho thấy xu hướng tăng tuổi thọ là một đường gần như thẳng.
Tổng kết nhận xét: “Phân tích hồi quy cho thấy một mối liên quan tuyến tính dương, rất mạnh (\(R^2 = 99.28\%\)) và có ý nghĩa thống kê cao (\(p < 0.001\)) giữa thời gian và tuổi thọ ở Việt Nam giai đoạn 1952-2007. Phương trình hồi quy cho thấy cứ mỗi năm trôi qua, tuổi thọ trung bình của người Việt Nam tăng thêm 0.649 năm.”
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
Nhập dữ liệu thủ công vào 2 vector
year_data <- c(1952, 1957, 1962, 1967, 1972, 1977,1982, 1987, 1992, 1997, 2002, 2007)
lifeExp_data <- 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)Tạo một data frame mới từ 2 vector trên
vn_data_manual <- data.frame(year = year_data, lifeExp = lifeExp_data)In ra để kiểm tra
print("Dữ liệu đã nhập:")## [1] "Dữ liệu đã nhập:"print(vn_data_manual)## year lifeExp
## 1 1952 40.4
## 2 1957 42.9
## 3 1962 45.4
## 4 1967 47.8
## 5 1972 50.3
## 6 1977 55.8
## 7 1982 58.8
## 8 1987 62.8
## 9 1992 67.7
## 10 1997 70.7
## 11 2002 73.0
## 12 2007 74.2Vẽ biểu đồ để xem xu hướng
print("Đang vẽ biểu đồ...")## [1] "Đang vẽ biểu đồ..."ggplot(vn_data_manual, aes(x = year, y = lifeExp)) +
# Vẽ các điểm (tô màu xanh)
geom_point(color = "blue", size = 3) +
# Thêm đường hồi quy (màu đỏ)
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Gia tăng tuổi thọ tại Việt Nam (1952-2007)",
x = "Năm",
y = "Tuổi thọ") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Xây dựng mô hình tuyến tính để “đánh giá”
đánh giá lifeExp phụ thuộc vào year
model_manual <- lm(lifeExp ~ year, data = vn_data_manual)
print("Kết quả đánh giá (Hồi quy tuyến tính):")## [1] "Kết quả đánh giá (Hồi quy tuyến tính):"summary(model_manual)##
## Call:
## lm(formula = lifeExp ~ year, data = vn_data_manual)
##
## 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) -1271.13492 43.77173 -29.04 0.0000000000547 ***
## year 0.67119 0.02211 30.35 0.0000000000353 ***
## ---
## 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: 0.00000000003527Để “đánh giá gia tăng tuổi thọ”, nhìn vào 3 con số sau: Estimate của year (\(= 6.489e-01\) hay \(0.649\)): Đây chính là tốc độ gia tăng. Nó có nghĩa là: Trung bình, cứ mỗi năm trôi qua, tuổi thọ của người Việt Nam gia tăng 0.649 năm. Pr(>|t|) của year (\(< 2e-16\)): Đây là p-value. Vì \(p < 0.05\) (thực tế là \(p\) gần như bằng 0), nó khẳng định rằng sự gia tăng 0.649 năm/năm này là có ý nghĩa thống kê rất cao (không phải do ngẫu nhiên). Multiple R-squared: 0.9928: \(R^2 = 99.28\%\). Con số này cho thấy mô hình (đường thẳng) giải thích được 99.28% sự thay đổi của tuổi thọ, nghĩa là yếu tố thời gian (year) dự đoán tuổi thọ rất chính xác trong giai đoạn này.
4.1. Nhập dữ liệu
Tạo các vector cho Y, X1, X2 và kết hợp chúng thành một data.frame (bảng dữ liệu) có tên là df.
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)Kết hợp thành một data frame tên là ‘df’
df <- data.frame(Y, X1, X2)In ra để kiểm tra
print("Dữ liệu cho Việc 4 (Ethylene oxide):")## [1] "Dữ liệu cho Việc 4 (Ethylene oxide):"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 04.2. Đánh giá mối liên quan giữa Y và X1
Chúng ta sẽ đánh giá mối liên quan này bằng cách trực quan hóa (biểu đồ tán xạ) và chạy mô hình hồi quy tuyến tính đơn biến (\(Y \sim X1\)).
4.2.1. Vẽ biểu đồ tán xạ (có màu - Dùng ggplot)
print("Đang vẽ biểu đồ Y ~ X1...")## [1] "Đang vẽ biểu đồ Y ~ X1..."ggplot(df, aes(x = X1, y = Y)) +
# Vẽ các điểm (tô màu xanh dương)
geom_point(color = "blue", size = 3) +
# Thêm đường hồi quy (màu đỏ) và khoảng tin cậy 95% (vùng xám)
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Mối liên quan giữa Ethylene oxide (Y) và Hoạt tính bạc (X1)",
x = "Hoạt tính bạc (X1)",
ylab = "Ethylene oxide (Y)") +
theme_minimal()## Ignoring unknown labels:
## • ylab : "Ethylene oxide (Y)"
## `geom_smooth()` using formula = 'y ~ x'
4.2.2. Chạy mô hình hồi quy và xem kết quả
model_YX1 <- lm(Y ~ X1, data = df)In kết quả tóm tắt chi tiết
print("Kết quả tóm tắt model Y ~ X1:")## [1] "Kết quả tóm tắt model Y ~ X1:"summary(model_YX1)##
## 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 0.000000554 ***
## 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.0004737Nhận xét (Đánh giá) kết quả (Y ~ X1):
Biểu đồ: Biểu đồ cho thấy các điểm màu xanh dương bám sát đường hồi quy đi xuống, cho thấy một mối tương quan âm (nghịch biến) và mạnh.
Multiple R-squared: 0.8115: \(R^2 = 81.15\%\). Điều này có nghĩa là \(X1\) (Hoạt tính bạc) giải thích được 81.15% sự thay đổi của \(Y\) (Ethylene oxide). Đây là mức độ giải thích rất cao.
Coefficients:Estimate (Hệ số) của X1: \(\approx -0.925\).Pr(>|t|) (p-value) của X1: \(\approx 0.000286\).
Kết luận: Vì p-value (\(0.000286\)) nhỏ hơn \(0.05\), chúng ta kết luận mối liên quan này là có ý nghĩa thống kê rất cao. Cứ \(X1\) tăng 1 đơn vị, \(Y\) giảm trung bình 0.925 đơn vị.
4.3. Đánh giá mối liên quan giữa Y và X2
Tương tự, chúng ta đánh giá mối liên quan bằng hồi quy đơn biến (\(Y \sim X2\)).
4.3.1. Vẽ biểu đồ tán xạ (có màu - Dùng ggplot)
print("Đang vẽ biểu đồ Y ~ X2...")## [1] "Đang vẽ biểu đồ Y ~ X2..."ggplot(df, aes(x = X2, y = Y)) +
# Vẽ các điểm (tô màu xanh lá)
geom_point(color = "green", size = 3) +
# Thêm đường hồi quy (màu đỏ)
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Mối liên quan giữa Ethylene oxide (Y) và Thời gian lưu (X2)",
x = "Thời gian lưu (X2)",
ylab = "Ethylene oxide (Y)") +
theme_minimal()## Ignoring unknown labels:
## • ylab : "Ethylene oxide (Y)"
## `geom_smooth()` using formula = 'y ~ x'
4.3.2. Chạy mô hình hồi quy và xem kết quả
model_YX2 <- lm(Y ~ X2, data = df)In kết quả tóm tắt chi tiết
print("Kết quả tóm tắt model Y ~ X2:")## [1] "Kết quả tóm tắt model Y ~ X2:"summary(model_YX2)##
## 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.01436Nhận xét (Đánh giá) kết quả (Y ~ X2):
Biểu đồ: Các điểm màu xanh lá có xu hướng đi lên, cho thấy một mối tương quan dương (đồng biến) và mạnh.
Multiple R-squared: 0.658: \(R^2 = 65.8\%\). \(X2\) (Thời gian lưu) giải thích được 65.8% sự thay đổi của \(Y\). Đây là mức độ giải thích khá cao. Coefficients:Estimate (Hệ số) của X2: \(\approx 1.259\).Pr(>|t|) (p-value) của X2: \(\approx 0.0049\).
Kết luận: Vì p-value (\(0.0049\)) nhỏ hơn \(0.05\), chúng ta kết luận mối liên quan này là có ý nghĩa thống kê. Cứ \(X2\) tăng 1 đơn vị, \(Y\) tăng trung bình 1.259 đơn vị. ### 4.4. Đánh giá mối liên quan độc lập giữa X2 và Y
Yêu cầu “đánh giá mối liên quan độc lập” của \(X2\) có nghĩa là: “X2 ảnh hưởng đến Y như thế nào sau khi đã loại bỏ (kiểm soát) ảnh hưởng của \(X1\)?”
Để trả lời câu hỏi này, chúng ta bắt buộc phải sử dụng mô hình hồi quy đa biến (\(Y \sim X1 + X2\)). Chúng ta sẽ xem xét p-value của \(X2\) trong mô hình mới này.
4.4.1. Xây dựng mô hình hồi quy ĐA BIẾN
Công thức Y ~ X1 + X2 có nghĩa là Y phụ thuộc vào X1 VÀ X2
model_multi <- lm(Y ~ X1 + X2, data = df)4.4.2. Xem kết quả tóm tắt mô hình đa biến
Đây là kết quả mấu chốt để trả lời câu 4.4
print("Kết quả tóm tắt model Đa biến (Y ~ X1 + X2):")## [1] "Kết quả tóm tắt model Đa biến (Y ~ X1 + X2):"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.002232Nhận xét kết quả (Tập trung vào X2):
Toàn bộ mô hình: \(R^2 = 0.8549\) (85.5%). Cả \(X1\) và \(X2\) cùng nhau giải thích được 85.5% sự thay đổi của \(Y\) (tốt hơn các mô hình đơn lẻ ở 4.2 và 4.3).
Hệ số X1: Vẫn có \(p = 0.00109\) ($ < 0.05$). \(X1\) là một yếu tố ảnh hưởng độc lập và có ý nghĩa.
Hệ số X2 (Câu trả lời cho 4.4):Estimate (Hệ số): Hệ số của \(X2\) giảm từ \(1.259\) (ở 4.3) xuống còn \(0.5367\).Pr(>|t|) (p-value): P-value của \(X2\) bây giờ là \(0.16541\).
Kết luận: Vì p-value của \(X2\) (\(0.165\)) lớn hơn mức ý nghĩa \(0.05\), chúng ta kết luận:Mối liên quan độc lập giữa thời gian lưu (\(X2\)) và Ethylene oxide (\(Y\)) là KHÔNG CÓ Ý NGHĨA THỐNG KÊ (sau khi đã kiểm soát ảnh hưởng của \(X1\)).
Diễn giải:
Mối liên quan mạnh mẽ mà chúng ta thấy ở bước 4.3 (khi \(X2\) đứng một mình) có thể là một mối liên quan “giả” (spurious) hoặc bị gây nhiễu (confounding). Có vẻ như cả \(X1\) và \(X2\) đều có liên quan với nhau, và \(X1\) mới là yếu tố ảnh hưởng chính. Khi ảnh hưởng mạnh của \(X1\) được đưa vào mô hình, \(X2\) không còn đóng góp thêm giá trị dự đoán nào đáng kể nữa.