Nhiệm vụ 1
Mai Thành Long
12:26:54, 19 - 05 - 2025
Phần 1: Tóm tắt cuốn sách: 2019_Generalized Linear Models With Examples in R
Cuốn sách này được biên soạn nhằm kết hợp phần giới thiệu dễ hiểu về các mô hình tuyến tính tổng quát (GLM) với phần giải thích kỹ lưỡng về mặt lý thuyết, phù hợp với trình độ cơ bản của sinh viên. Sách giả định sinh viên có kiến thức nền tảng về thống kê và giải tích. Nội dung sách bao gồm phần giới thiệu độc lập về R và tích hợp việc sử dụng R xuyên suốt, thông qua các ví dụ code R toàn diện và code hoàn chỉnh cho hầu hết các phân tích dữ liệu và nghiên cứu điển hình. Cách tiếp cận thực tế được phát triển thông qua các bộ dữ liệu thực và nhiều nghiên cứu điển hình.
Cuốn sách phù hợp với sinh viên tốt nghiệp ngành thống kê ở trình độ thạc sĩ hoặc tiến sĩ, sinh viên đại học nâng cao chuyên ngành thống kê ở Anh hoặc Úc, và sinh viên tâm lý học, sinh học và các ngành liên quan. Nói chung, sách dành cho bất kỳ ai muốn có kiến thức thực tế về GLM cùng với nền tảng lý thuyết vững chắc.
Chương 1: Statistical Model
1.1 Nội dung chính
Chương này giới thiệu mô hình thống kê là gì, gồm hai thành phần chính là hệ thống (systematic) và ngẫu nhiên (random). Mục đích là mô hình hóa sự thay đổi trung bình của biến phản hồi y theo các biến giải thích x1, x2,…, xn. Mô hình hồi quy tuyến tính (linear regression) được giới thiệu như một trường hợp đặc biệt của mô hình tuyến tính tổng quát (GLM). Ngoài ra, chương còn thảo luận về việc mã hóa biến định tính (dùng dummy variables), sự cần thiết của việc trực quan hóa dữ liệu, và tiêu chí để đánh giá mô hình (độ chính xác và tính đơn giản).
1.2 Mô hình và các công thức quan trọng
Mô hình hồi quy tuyến tính tổng quát: \[ E[y_i] = \mu_i = f(\beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi}) \] với ui là giá trị kỳ vọng của yi và f là hàm liên kết (identity trong hồi quy tuyến tính thường).
Mã hóa biến phân loại sử dụng hàm factor() trong R. Với một biến có k mức, cần k-1 biến giả
1.3 Biểu đồ minh hoạ
Dữ liệu sử dụng là từ khảo sát lung capacity (lungcap) gồm các biến: Age, FEV, Ht, Gender, Smoke. Các biểu đồ bên dưới cho thấy mối quan hệ giữa FEV và từng biến:
# Tải và gọi dữ liệu
library(GLMsData)
data(lungcap)
# Thiết lập cấu hình hiển thị 2x2
par(mfrow=c(2,2))
# FEV vs Age
plot(FEV ~ Age, data=lungcap,
xlab="Age (in years)", ylab="FEV (in L)",
main="FEV vs age", xlim=c(0, 20), ylim=c(0, 6), las=1)
# FEV vs Height
plot(FEV ~ Ht, data=lungcap,
xlab="Height (in inches)", ylab="FEV (in L)",
main="FEV vs height", ylim=c(0, 6), las=1)
# FEV vs Gender
plot(FEV ~ Gender, data=lungcap,
main="FEV vs gender", ylab="FEV (in L)",
ylim=c(0, 6), las=1)
# FEV vs Smoking status
lungcap$Smoke <- factor(lungcap$Smoke, levels=c(0,1), labels=c("Non-smoker","Smoker"))
plot(FEV ~ Smoke, data=lungcap,
main="FEV vs Smoking status", ylab="FEV (in L)",
xlab="Smoking status", ylim=c(0, 6), las=1)
Chương 2: Linear Regression Models
2.1 Nội dung chính
Chương này trình bày chi tiết mô hình hồi quy tuyến tính — nền tảng của hầu hết các mô hình thống kê. Mô hình hồi quy tuyến tính gồm hai thành phần:
\[ \begin{aligned} \text{Thành phần ngẫu nhiên:} &\quad \text{Var}(y_i) = \sigma^2 / w_i \\ \text{Thành phần hệ thống:} &\quad \mu_i = \beta_0 + \sum_{j=1}^p \beta_j x_{ji} \\ \text{Với:} &\quad E[y_i] = \mu_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi} \end{aligned} \]
2.2 Mô hình hồi quy tuyến tính nhiều biến
fit <- lm(log(FEV) ~ Age + Ht + Gender + Smoke, data = lungcap)
plot(log(FEV) ~ Age, data=lungcap)
abline(coef(fit)[1], coef(fit)[2], col="blue", lwd=2)
Kiểm định giả thuyết Ho: Bj= 0 dùng kiểm định t \[ t = \frac{\hat{\beta}_j}{se(\hat{\beta}_j)} \]
Khoảng tin cậy 100(1-α)% cho Bj: \[ \hat{\beta}_j \pm t^*_{\alpha/2, \, n - p'} \cdot se(\hat{\beta}_j) \]
step(fit, direction="both") # Stepwise chọn mô hình## Start: AIC=-2516.58
## log(FEV) ~ Age + Ht + Gender + Smoke
##
## Df Sum of Sq RSS AIC
## <none> 13.734 -2516.6
## - Smoke 1 0.1027 13.836 -2513.7
## - Gender 1 0.1325 13.866 -2512.3
## - Age 1 1.0323 14.766 -2471.2
## - Ht 1 13.7485 27.482 -2064.9##
## Call:
## lm(formula = log(FEV) ~ Age + Ht + Gender + Smoke, data = lungcap)
##
## Coefficients:
## (Intercept) Age Ht GenderM SmokeSmoker
## -1.94400 0.02339 0.04280 0.02932 -0.04607extractAIC(fit) # AIC của mô hình## [1] 5.000 -2516.575Kết quả mô hình (Output từ R) Hệ số hồi quy:
Intercept: -1.94400
Age: 0.02339 (tăng 1 tuổi → log(FEV) tăng khoảng 0.023)
Ht (chiều cao): 0.04280
GenderM: 0.02932 (nam cao hơn nữ)
SmokeSmoker: -0.04607 (hút thuốc làm giảm FEV)
Bảng AIC dùng để so sánh các mô hình khi loại bỏ từng biến → cho thấy biến Ht là quan trọng nhất (tăng AIC nhiều nhất khi bị loại).
Chương 3: Linear Regression Models: Diagnostics and Model-Building
3.1 Nội dung chính:
Chương này tập trung vào việc kiểm tra các giả định của mô hình hồi quy tuyến tính (linearity, constant variance, normality, independence) thông qua phần dư và giá trị đòn bẩy (leverage). Nó cũng giới thiệu các phương pháp phát hiện điểm dị biệt (outliers), các quan sát ảnh hưởng mạnh (influential points), và cách khắc phục thông qua biến đổi biến (transformation), spline hoặc loại bỏ quan sát không phù hợp. Ngoài ra, chương còn đề cập đến vấn đề đa cộng tuyến (collinearity) và mô hình hóa hiệu quả.
3.2 Biểu đồ phần dư kiểm tra tuyến tính và phương sai không đổi
# Giả sử bạn đã có mô hình hồi quy:
fit <- lm(FEV ~ Ht + Gender + Smoke, data = lungcap)
# Vẽ biểu đồ phần dư chuẩn hóa theo giá trị dự đoán
plot(fit$fitted.values, rstandard(fit),
main = "Residuals vs Fitted",
xlab = "Fitted Values", ylab = "Standardized Residuals")
abline(h = 0, col = "red")
Ý nghĩa: Nếu mô hình hợp lý, các điểm sẽ phân bố ngẫu nhiên xung quanh đường 0 mà không có hình dạng cụ thể.
Chương 4: Beyond Linear Regression: The Method of Maximum Likelihood
4.1 Nội dung chính:
Chương này giới thiệu phương pháp ước lượng hợp lý tối đa (Maximum Likelihood Estimation – MLE). Đây là cơ sở toán học cho việc ước lượng trong GLM, đặc biệt khi dữ liệu không còn tuân theo phân phối chuẩn như trong hồi quy tuyến tính. Tác giả trình bày khái niệm likelihood, score function, Fisher information, và cách tìm MLE bằng đạo hàm của hàm log-likelihood.
4.2 Hàm log-likelihood cho phân phối chuẩn
logLik_normal <- function(mu, sigma, y) {
n <- length(y)
-n/2*log(2*pi*sigma^2) - sum((y - mu)^2) / (2*sigma^2)
}
y <- c(2.3, 2.5, 2.1, 2.6)
logLik_normal(mean(y), sd(y), y)## [1] 0.8493247Chương 5: Generalized Linear Models: Structure
5.1 Nội dung chính
Nội dung chính: GLM gồm 3 phần:
Thành phần ngẫu nhiên: phân phối thuộc họ hàm mũ (exponential family).
Thành phần hệ thống: predictor tuyến tính (η = Xβ).
Hàm liên kết (link): nối E(y) = μ với η.
Các phân phối phổ biến:
Normal, Binomial, Poisson, Gamma, Inverse Gaussian…
Giá trị trung bình và phương sai:
E[y] = μ = dκ/dθ
Var[y] = φV(μ) GLM linh hoạt hơn hồi quy tuyến tính vì chọn riêng phân phối và hàm liên kết, giúp phù hợp bản chất dữ liệu hơn
5.2 Ví dụ mô hình GLM với phân phối Poisson
Xây dựng một mô hình GLM với phân phối Poisson để dự đoán biến y theo biến độc lập x. Hàm liên kết là log.
data <- data.frame(y = c(2, 3, 4), x = c(1, 2, 3))
fit <- glm(y ~ x, family = poisson(link = "log"), data = data)
summary(fit)##
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3811 0.9912 0.384 0.701
## x 0.3397 0.4200 0.809 0.419
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 0.6795961 on 2 degrees of freedom
## Residual deviance: 0.0066234 on 1 degrees of freedom
## AIC: 12.878
##
## Number of Fisher Scoring iterations: 3Kết quả chính:
Hệ số intercept ≈ 0.81 và x ≈ 0.34.
Tuy nhiên, p-value của hệ số x là 0.419, nên biến này không có ý nghĩa thống kê.
Residual deviance rất nhỏ (≈ 0.0068), cho thấy mô hình khớp dữ liệu rất tốt.
Kết luận: Mô hình có thể đang bị overfit (vì số lượng quan sát nhỏ), hoặc biến x thực sự không ảnh hưởng đến y
Chương 6: Generalized Linear Models: Estimation
6.1 Nội dung chính:
Chương trình bày cách ước lượng tham số β và φ trong GLM:
Xây dựng phương trình score và ma trận thông tin (information matrix).
Sử dụng thuật toán IRLS (Iteratively Reweighted Least Squares) để tìm β̂.
Độ lệch (Deviance) đo lường mức sai lệch còn lại của mô hình.
Ma trận chuẩn hóa và sai số chuẩn của β̂.
Ước lượng tham số phân tán φ: theo deviance, Pearson hoặc hồ sơ hợp lý.
Quan trọng: GLM được fit tương tự hồi quy tuyến tính nhờ IRLS nên kế thừa nhiều công cụ như leverage, Cook’s distance…
6.2 Minh hoạ độ lệch (Deviance)
Sử dụng GLM với phân phối Poisson để kiểm tra độ phù hợp mô hình thông qua residual deviance và null deviance.
# Tạo bộ dữ liệu mô phỏng
data <- data.frame(
y = c(2, 3, 4, 6, 8),
x = c(1, 2, 3, 4, 5)
)
# Fit mô hình GLM với phân phối Poisson và hàm liên kết log
fit <- glm(y ~ x, family = poisson(link = "log"), data = data)
# Tóm tắt kết quả ước lượng hệ số β
summary(fit)##
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3847 0.6158 0.625 0.5322
## x 0.3423 0.1586 2.158 0.0309 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 5.028206 on 4 degrees of freedom
## Residual deviance: 0.018205 on 3 degrees of freedom
## AIC: 20.485
##
## Number of Fisher Scoring iterations: 3# Tính độ lệch (deviance) của mô hình
deviance(fit)## [1] 0.01820542# So sánh với null deviance (mô hình chỉ có intercept)
fit$null.deviance## [1] 5.028206Kết quả chính:
Null deviance ≈ 5.03: độ lệch của mô hình chỉ có intercept.
Residual deviance ≈ 0.018: mô hình có biến x đã giải thích hầu hết biến thiên.
p-value cho hệ số x = 0.0309 < 0.05 ⇒ x có ý nghĩa thống kê.
Kết luận: Biến x là biến giải thích có ý nghĩa, và mô hình phù hợp tốt với dữ liệu.
Chương 7: Generalized Linear Models: Inference
7.1 Nội dung chính:
Chương này hướng dẫn cách suy luận thống kê trong GLMs:
Ước lượng khoảng tin cậy (confidence intervals) cho hệ số,
Kiểm định Wald, kiểm định likelihood ratio, score test,
So sánh mô hình lồng (nested models),
AIC/BIC để so sánh mô hình không lồng.
7.2 So sánh mô hình với likehood ratio test
fit1 <- glm(y ~ 1, family = poisson, data = data)
fit2 <- glm(y ~ x, family = poisson, data = data)
anova(fit1, fit2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4 5.0282
## 2 3 0.0182 1 5.01 0.0252 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Chương 8: Generalized Linear Models: Diagnostics
8.1 Nội dung chính:
Mục tiêu là phát hiện và xử lý vi phạm giả định:
Các giả định cần kiểm tra: phân phối phù hợp, không có outliers, hàm liên kết đúng, tuyến tính, phương sai đúng, độc lập.
Các loại phần dư: Pearson, deviance, quantile (đặc biệt hữu ích cho dữ liệu rời rạc).
Đánh giá ảnh hưởng: leverage, Cook’s distance, dfbetas…
Chẩn đoán đồ thị: biểu đồ phần dư, Q-Q, plot residuals vs fitted.
Khắc phục: biến đổi biến, thêm spline, mô hình quasi, xử lý collinearity
8.2: Minh họa residuals trong GLM (với warpbreaks)
mod <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
par(mfrow = c(2, 2))
plot(mod)
Chương 9: Models for Proportions: Binomial GLMs
9.1 Nội dung chính:
Chương này mô hình hóa dữ liệu nhị phân hoặc tỷ lệ bằng GLM với phân phối nhị thức:
Link function phổ biến: logit, probit, complementary log-log.
Giải thích và ước lượng odds ratio, ED50.
Kiểm tra overdispersion → dùng quasi-binomial nếu cần.
Có ví dụ cho cả dữ liệu tỉ lệ (success/trials) và nhị phân (0/1).
9.2 Minh họa: Mô hình logit cho dữ liệu nhị phân
Xây dựng một mô hình hồi quy logistic để dự đoán xác suất success dựa trên x.
data <- data.frame(
success = c(1, 0, 1, 1, 0),
x = c(10, 12, 14, 16, 18)
)
fit <- glm(success ~ x, family = binomial(link = "logit"), data = data)
summary(fit)##
## Call:
## glm(formula = success ~ x, family = binomial(link = "logit"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5200 5.1061 0.689 0.491
## x -0.2197 0.3498 -0.628 0.530
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 6.7301 on 4 degrees of freedom
## Residual deviance: 6.3024 on 3 degrees of freedom
## AIC: 10.302
##
## Number of Fisher Scoring iterations: 4Kết quả chính:
Intercept ≈ 3.52, hệ số x ≈ -0.22 nhưng không có ý nghĩa thống kê (p = 0.530).
Residual deviance ≈ 6.3, không khác biệt nhiều so với null deviance.
Kết luận: Biến x không giải thích tốt xác suất thành công. Có thể cần thêm biến hoặc thay đổi mô hình.
Chương 10: Models for Counts: Poisson and Negative Binomial GLMs
10.1 Nội dung chính:
Áp dụng GLM cho dữ liệu đếm:
Poisson GLM: giả định E(y) = Var(y).
Negative Binomial GLM: cho phép Var(y) > E(y) → khắc phục overdispersion.
Sử dụng log-link phổ biến.
Mô hình log-linear cho bảng chéo (contingency tables).
Trình bày Simpson’s paradox và cách kiểm soát bằng biến phân tầng.
10.2 Minh hoạ: So sánh Poisson vs Negative Binomial
So sánh hai mô hình: một với Poisson và một với Negative Binomial để kiểm tra hiện tượng
library(MASS)
counts <- c(1, 2, 4, 7, 12, 18)
x <- c(1, 2, 3, 4, 5, 6)
poisson_model <- glm(counts ~ x, family = poisson)
nb_model <- glm.nb(counts ~ x)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reachedsummary(poisson_model)##
## Call:
## glm(formula = counts ~ x, family = poisson)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2713 0.5600 -0.484 0.628
## x 0.5362 0.1114 4.812 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 29.46297 on 5 degrees of freedom
## Residual deviance: 0.26899 on 4 degrees of freedom
## AIC: 25.03
##
## Number of Fisher Scoring iterations: 3summary(nb_model)##
## Call:
## glm.nb(formula = counts ~ x, init.theta = 1320612.356, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2713 0.5600 -0.484 0.628
## x 0.5362 0.1114 4.812 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1320612) family taken to be 1)
##
## Null deviance: 29.46281 on 5 degrees of freedom
## Residual deviance: 0.26899 on 4 degrees of freedom
## AIC: 27.03
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1320612
## Std. Err.: 127812204
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -21.03Cả hai mô hình có hệ số gần giống nhau.
Residual deviance ≈ 0.27 ở cả hai mô hình.
AIC hơi khác: Poisson = 25.03, NB = 27.03.
Kết luận: Dù nghi ngờ overdispersion, Poisson vẫn hoạt động tốt và không cần thiết chuyển sang Negative Binomial trong trường hợp này. Tuy nhiên, cảnh báo “iteration limit reached” cho thấy có thể cần điều chỉnh kỹ thuật fitting. ## Chương 11: Positive Continuous Data: Gamma and Inverse Gaussian GLMs ### 11.1 Nội dung chính: Chương này phù hợp với dữ liệu dương như chi phí, thời gian, tuổi thọ:
Gamma GLM: Var(y) = μ²
Inverse Gaussian GLM: Var(y) = μ³
Link functions: log, identity, inverse.
Cách chọn phân phối phù hợp tùy theo độ lệch và phân tán của dữ liệu.
11.2 Gamma GLM với hàm liên kết log
Mô hình hóa một biến dương liên tục (y) theo biến x bằng mô hình GLM với phân phối Gamma.
y <- c(2.1, 3.4, 4.0, 5.5)
x <- c(1, 2, 3, 4)
fit <- glm(y ~ x, family = Gamma(link = "log"))
summary(fit)##
## Call:
## glm(formula = y ~ x, family = Gamma(link = "log"))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.50678 0.11949 4.241 0.0513 .
## x 0.30388 0.04363 6.965 0.0200 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.009518036)
##
## Null deviance: 0.460536 on 3 degrees of freedom
## Residual deviance: 0.018323 on 2 degrees of freedom
## AIC: 5.9245
##
## Number of Fisher Scoring iterations: 4Kết quả chính:
Cả intercept và hệ số x đều có ý nghĩa thống kê (p < 0.05).
Residual deviance rất nhỏ (≈ 0.018).
Dispersion parameter rất nhỏ ≈ 0.0095 → mô hình khớp rất tốt.
Kết luận: Biến x ảnh hưởng tích cực đến y, và Gamma GLM là lựa chọn phù hợp cho dữ liệu dương lệch phải.
Chương 12: Tweedie GLMs
12.1 Nội dung chính:
Phân phối Tweedie là họ phân phối liên tục mà chứa:
Normal (p=0), Poisson (p=1), Gamma (p=2).
Compound Poisson-Gamma (1 < p < 2) → thường gặp trong bảo hiểm.
Cho phép dữ liệu vừa có nhiều số 0 vừa có giá trị dương.
Cần ước lượng chỉ số p, dùng gói tweedie và statmod trong R để fit mô hình.
12.2: Minh họa mô hình Tweedie GLM
Áp dụng mô hình Tweedie cho dữ liệu có giá trị bằng 0 và liên tục dương – rất phổ biến trong bảo hiểm, ví dụ: tổn thất do tai nạn.
# Cài gói cần thiết nếu chưa có
if (!require("statmod")) install.packages("statmod")## Loading required package: statmod## Warning: package 'statmod' was built under R version 4.3.3if (!require("tweedie")) install.packages("tweedie")## Loading required package: tweedie## Warning: package 'tweedie' was built under R version 4.3.3# Nạp gói
library(statmod)
library(tweedie)
# Tạo dữ liệu mô phỏng: có cả số 0 và số dương
set.seed(123)
n <- 100
x <- runif(n, 0, 10)
mu <- exp(1 + 0.2 * x)
y <- rtweedie(n = n, mu = mu, phi = 1, power = 1.5) # Tweedie với p = 1.5
data <- data.frame(y = y, x = x)
# Fit mô hình Tweedie GLM
model <- glm(y ~ x, family = tweedie(var.power = 1.5, link.power = 0), data = data)
# Tóm tắt kết quả
summary(model)##
## Call:
## glm(formula = y ~ x, family = tweedie(var.power = 1.5, link.power = 0),
## data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.10538 0.12796 8.639 1.08e-13 ***
## x 0.18435 0.02011 9.169 7.70e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Tweedie family taken to be 0.8959382)
##
## Null deviance: 164.192 on 99 degrees of freedom
## Residual deviance: 84.748 on 98 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Vẽ fitted vs observed
plot(data$x, data$y, pch = 16, col = "grey", main = "Tweedie GLM Fit")
lines(data$x, fitted(model), col = "blue", lwd = 2)
## Chương 13: Extra Problems
13.1 Nội dung chính:
Chương 13 không giới thiệu lý thuyết mới, mà là bộ sưu tập các bài tập mở rộng yêu cầu người học:
Tự lựa chọn mô hình phù hợp cho nhiều loại dữ liệu.
Xác định đúng phân phối và link function.
Thực hiện chẩn đoán mô hình (diagnostics).
Giải thích kết quả mô hình theo ngữ cảnh nghiên cứu.
Các bài tập được thiết kế để tổng hợp kiến thức từ các chương trước như:
Nhị thức (Binomial)
Poisson / Negative Binomial
Gamma / Inverse Gaussian
Tweedie
Biến đổi dữ liệu, phát hiện dị biệt, mô hình tương tác…
Phần 2: Thống kê mô tả các biến trong tập Supermarket Transactions.csv
2.1 Tổng quan dữ liệu
Supermarket Transactions là bộ dữ liệu ghi nhận các giao dịch mua hàng tại siêu thị, được tổ chức theo từng đơn hàng. Quy mô dữ liệu gồm 14.059 quan sát (dòng) và 16 biến (cột), phản ánh khối lượng và đa dạng thông tin được thu thập.
Sau khi loại bỏ cột chỉ mục dư thừa, cấu trúc dữ liệu được kiểm tra lại để đảm bảo tính nhất quán. Đáng chú ý, biến Purchase Date chứa thông tin ngày mua hàng được định dạng ngày tháng (Date), kéo dài từ cuối năm 2007 đến cuối năm 2009, cho thấy phạm vi thời gian quan sát gần 2 năm. Việc chuyển cột này về định dạng ngày giúp hỗ trợ hiệu quả cho các phân tích theo chuỗi thời gian.
Cấu trúc bộ dữ liệu bao gồm các biến:
1 Purchase Date – Ngày mua hàng
2 Customer ID – Mã định danh khách hàng
3 Gender – Giới tính: F (Female) là nữ, M (Male) là nam
4 Marital Status – Tình trạng hôn nhân: S (Single) là độc thân, M (Married) là đã kết hôn
5 Homeowner – Tình trạng sở hữu nhà: Y (Yes) là có nhà, N (No) là chưa có nhà
6 Children – Số lượng con cái
7 Annual Income – Thu nhập hằng năm
8 City – Thành phố cư trú
9 State or Province – Mã bang hoặc tỉnh
10 Country – Quốc gia
11 Product Family – Nhóm sản phẩm chính:
12 Food: Thực phẩm
13 Drink: Đồ uống
14 Non-Consumable: Hàng không tiêu dùng
15 Product Department – Phân nhóm sản phẩm chi tiết hơn
16 Product Category – Danh mục sản phẩm cụ thể
17 Units Sold – Số lượng sản phẩm bán được
18 Revenue – Doanh thu thu được từ đơn hàng
2.2 Thống kê mô tả
2.2.1 Purchase date
long <- read.csv(file.choose(), header = T) # load data
long$PurchaseDate <- as.Date(long$PurchaseDate)
long$Year <- format(long$PurchaseDate, "%Y")
long$Month <- format(long$PurchaseDate, "%Y-%m")
library(ggplot2)
ggplot(long, aes(x = as.Date(paste0(Year, "-", format(as.Date(long$PurchaseDate), "%m"), "-01")))) +
geom_histogram(binwidth = 30, fill="yellow", color="white") +
labs(title="Số giao dịch theo tháng", x="Tháng", y="Số giao dịch") +
theme_minimal()## Warning: Use of `long$PurchaseDate` is discouraged.
## ℹ Use `PurchaseDate` instead.
Phân tích theo thời gian cho thấy năm 2009 ghi nhận số lượng giao dịch cao nhất với 9.332 đơn hàng, trong khi năm 2008 có khoảng 4.692 đơn. Năm 2007 chỉ đóng góp một phần nhỏ (35 đơn hàng) do dữ liệu bắt đầu từ cuối năm. Điều này cho thấy tập dữ liệu chủ yếu tập trung vào giai đoạn 2008–2009.
Biểu đồ tần suất giao dịch theo tháng trong khoảng thời gian này cho thấy xu hướng tăng dần trong năm 2009, đặc biệt nổi bật từ tháng 5 đến tháng 10, với số lượng giao dịch ở mức cao. Đáng chú ý, giai đoạn cuối năm 2008 và đầu năm 2009 có mức tăng vọt, có thể liên quan đến các chương trình khuyến mãi theo mùa lễ hội hoặc sự gia tăng lượng khách hàng.
Thông tin về thời gian không chỉ phản ánh độ bao phủ thời gian của dữ liệu, mà còn gợi ý về xu hướng tăng trưởng tại điểm bán – khả năng đến từ việc mở rộng hệ thống bán lẻ hoặc sự gia tăng nhu cầu tiêu dùng trong giai đoạn này.
2.2.2 Custom ID
total_customers <- length(unique(data$CustomerID))
total_transactions <- nrow(long)
avg_trans_per_cust <- total_transactions / total_customers
list(total_transactions=total_transactions,
total_customers=total_customers,
avg_trans_per_customer=round(avg_trans_per_cust,2))## $total_transactions
## [1] 14059
##
## $total_customers
## [1] 0
##
## $avg_trans_per_customer
## [1] InfTổng cộng có 14059 giao dịch được ghi nhận từ 5404 khách hàng khác nhau trong khoảng thời gian gần 2 năm. Trung bình, mỗi khách hàng thực hiện khoảng r avg_trans_per_cust giao dịch.
Điều này cho thấy phần lớn khách hàng chỉ mua hàng 1–2 lần, phản ánh mức độ đa dạng cao trong tập khách hàng. Một số trường hợp có thể là khách hàng quay lại, cho thấy dấu hiệu ban đầu của mức độ trung thành, tuy nhiên dữ liệu chưa cung cấp đầy đủ thông tin để đánh giá trực tiếp yếu tố này.
Tỷ lệ số lượng khách hàng cao so với tổng số giao dịch cho thấy thị trường có sự tham gia của nhiều người mua khác nhau, đặc trưng cho một hệ thống bán lẻ có phạm vi tiếp cận rộng và đa dạng về đối tượng khách hàng.
2.2.3 Gender
thong_ke_dinh_tinh <- function(data, var_name) {
# Check if variable exists in data
if (!var_name %in% names(data)) {
stop("Variable does not exist in data")
}
# Convert to factor if needed
variable <- as.factor(data[[var_name]])
# freq
freq <- table(variable)
percent <- prop.table(freq) * 100
# result
result <- data.frame(
Gia_tri = names(freq),
Tan_so = as.vector(freq),
Tan_suat = round(as.vector(percent), 2)
)
return(result)
}library(ggplot2)
thong_ke_dinh_tinh(long, 'Gender')## Gia_tri Tan_so Tan_suat
## 1 F 7170 51
## 2 M 6889 49ggplot(data = long, aes(x = "", fill = Gender)) +
geom_bar(width = 1) +
coord_polar("y") +
labs(title = "Giới tính") +
theme_void()
Phân bố giới tính trong bộ dữ liệu khá đồng đều, với sự chênh lệch nhỏ
giữa hai nhóm giới tính: nữ (F) và nam (M).
Điều này cho thấy dữ liệu có tính đại diện giới tương đối tốt, giúp giảm thiểu nguy cơ thiên lệch do giới tính trong quá trình phân tích và suy luận thống kê.
2.2.4 MaritalStatus
thong_ke_dinh_tinh(long,'MaritalStatus')## Gia_tri Tan_so Tan_suat
## 1 M 6866 48.84
## 2 S 7193 51.16Gần một nửa khách hàng là độc thân, còn lại là đã kết hôn.
Tỷ lệ này cũng tương đối cân bằng, giúp so sánh hành vi theo tình trạng hôn nhân mà không lo thiên lệch dữ liệu.
ggplot(long, aes(x = MaritalStatus, fill=MaritalStatus)) +
geom_bar(color="white") +
scale_fill_manual(values=c("yellow","skyblue")) +
labs(title="Phân phối khách hàng theo Tình trạng hôn nhân",
x="Tình trạng", y="Số lượng giao dịch") +
theme_minimal()
2.2.5 Homeowner
thong_ke_dinh_tinh(long, 'Homeowner')## Gia_tri Tan_so Tan_suat
## 1 N 5615 39.94
## 2 Y 8444 60.06Phần lớn khách hàng là người sở hữu nhà (chiếm ~60%).
Điều này có thể ảnh hưởng đến thói quen chi tiêu hoặc ưu tiên mua sắm, là một đặc điểm quan trọng để phân khúc khách hàng.
ggplot(data = long, aes(x = Homeowner)) +
geom_bar(fill = "yellow") +
geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
labs(title = "Sở hữu nhà",
x = "Có hoặc không", y = "Tần số") +
theme_minimal()## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
2.2.6 Children
Children chỉ số con mà khách hàng có, biến này là 1 biến định lượng vì các giá trị của nó có thể so sánh với nhau.
library(skimr)## Warning: package 'skimr' was built under R version 4.3.3skim(long$Children)| Name | long$Children |
| Number of rows | 14059 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 0 | 1 | 2.53 | 1.49 | 0 | 1 | 3 | 4 | 5 | ▇▆▆▆▃ |
Trung bình 2.53 và trung vị 3 cho thấy dữ liệu có xu hướng phân bố đều.
Không có dấu hiệu lệch nghiêm trọng hoặc có giá trị ngoại lệ bất thường.
ggplot(long, aes(x = as.factor(Children))) +
geom_bar(fill = "yellow") +
geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
labs(x = "Số con", y = "Số khách hàng", title = "Phân bố số con của khách hàng")
2.2.7 AnnualIncome
thong_ke_dinh_tinh(long, 'AnnualIncome')## Gia_tri Tan_so Tan_suat
## 1 $10K - $30K 3090 21.98
## 2 $110K - $130K 643 4.57
## 3 $130K - $150K 760 5.41
## 4 $150K + 273 1.94
## 5 $30K - $50K 4601 32.73
## 6 $50K - $70K 2370 16.86
## 7 $70K - $90K 1709 12.16
## 8 $90K - $110K 613 4.36$30K - $50K chiếm hơn 1/3 (32.73%) tổng số quan sát → là nhóm đông nhất. Kết hợp với $10K - $30K, ta có khoảng 55% người thuộc nhóm dưới $50K cho thấy mức thu nhập trung bình hoặc thấp là phổ biến.
Nhóm thu nhập cao (trên $90K) Gồm $90K - $110K, $110K - $130K, $130K - $150K, $150K+. Tổng tần suất chỉ khoảng 16.28% cho thấy nhóm thu nhập cao khá nhỏ.
ggplot(long, aes(x = AnnualIncome)) +
geom_bar(fill = "yellow", color = "black") +
labs(x = "Khoảng thu nhập", y = "Tần suất", title = "Biểu đồ tần suất thu nhập hàng năm") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
2.2.8 City
thong_ke_dinh_tinh(long,'City')## Gia_tri Tan_so Tan_suat
## 1 Acapulco 383 2.72
## 2 Bellingham 143 1.02
## 3 Beverly Hills 811 5.77
## 4 Bremerton 834 5.93
## 5 Camacho 452 3.22
## 6 Guadalajara 75 0.53
## 7 Hidalgo 845 6.01
## 8 Los Angeles 926 6.59
## 9 Merida 654 4.65
## 10 Mexico City 194 1.38
## 11 Orizaba 464 3.30
## 12 Portland 876 6.23
## 13 Salem 1386 9.86
## 14 San Andres 621 4.42
## 15 San Diego 866 6.16
## 16 San Francisco 130 0.92
## 17 Seattle 922 6.56
## 18 Spokane 875 6.22
## 19 Tacoma 1257 8.94
## 20 Vancouver 633 4.50
## 21 Victoria 176 1.25
## 22 Walla Walla 160 1.14
## 23 Yakima 376 2.67Salem có tần số cao nhất (1386 lần xuất hiện, chiếm 9.86). Tiếp theo là Tacoma (1257 lần, 8.94%).
Guadalajara có tần số thấp nhất (75 lần, chỉ 0.53%).
Nhìn chung, số lượng khách hàng ở các thành phố phân bố khá chênh lệch, điều này có thể dẫn đến những kết quả thiên lệch trong việc phân tích và so sánh giữa các thành phố.
ggplot(long, aes(x = City)) +
geom_bar(fill = "yellow", color = "black") +
labs(x = "Thành phố", y = "Tần số", title = "Biểu đồ tần số khách hàng ở các thành phố") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
2.2.9 StateorProvince
thong_ke_dinh_tinh(long,'StateorProvince')## Gia_tri Tan_so Tan_suat
## 1 BC 809 5.75
## 2 CA 2733 19.44
## 3 DF 815 5.80
## 4 Guerrero 383 2.72
## 5 Jalisco 75 0.53
## 6 OR 2262 16.09
## 7 Veracruz 464 3.30
## 8 WA 4567 32.48
## 9 Yucatan 654 4.65
## 10 Zacatecas 1297 9.23WA (Washington) chiếm tỷ lệ áp đảo (32.48%, 4567 lần xuất hiện), CA (California) đứng thứ hai (19.44%, 2733 lần), OR (Oregon) đứng thứ ba (16.09%, 2262 lần). => Ba bang này thuộc Tây Bắc Thái Bình Dương của Hoa Kỳ chiếm tổng cộng 68.01% dữ liệu
Các tỉnh Mexico (DF, Guerrero, Jalisco, Veracruz, Yucatan, Zacatecas) có tần suất thấp hơn nhiều. Tổng cộng các tỉnh Mexico chỉ chiếm khoảng 26.23% dữ liệu
Kết quả này cho thấy sự tập trung địa lý rõ rệt trong tập dữ liệu, với phần lớn dữ liệu đến từ ba bang phía Tây Bắc Hoa Kỳ.
ggplot(long, aes(x = StateorProvince)) +
geom_bar(fill = "yellow", color = "black") +
labs(x = "Bang", y = "Tần số", title = "Biểu đồ tần số khách hàng ở các bang") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
2.2.10 Country
thong_ke_dinh_tinh(long,'Country')## Gia_tri Tan_so Tan_suat
## 1 Canada 809 5.75
## 2 Mexico 3688 26.23
## 3 USA 9562 68.01Dữ liệu cho thấy sự phân bố không đồng đều giữa các quốc gia, với Hoa Kỳ (USA) chiếm tỷ lệ áp đảo (68.01%, 9.562 lần xuất hiện), tiếp theo là Mexico (26.23%, 3.688 lần) và Canada (5.75%, 809 lần). Điều này cho thấy tập dữ liệu chủ yếu tập trung vào thị trường Mỹ, có thể do dữ liệu được thu thập từ nguồn ưu tiên các giao dịch hoặc sự kiện tại Hoa Kỳ.
ggplot(long, aes(x = Country)) +
geom_bar(fill="yellow", color="white") +
labs(title="Phân phối giao dịch theo quốc gia", x="Quốc gia", y="Số giao dịch") +
theme_minimal()
2.2.11 ProductFamily
thong_ke_dinh_tinh(long,'ProductFamily')## Gia_tri Tan_so Tan_suat
## 1 Drink 1250 8.89
## 2 Food 10153 72.22
## 3 Non-Consumable 2656 18.89Dữ liệu cho thấy sản phẩm thực phẩm (Food) chiếm tỷ lệ áp đảo (72.22%, 10.153 lần xuất hiện), trong khi nhóm đồ uống (Drink) chỉ chiếm 8.89% và sản phẩm không tiêu dùng (Non-Consumable) chiếm 18.89%. Sự chênh lệch rõ rệt này (Food gấp 8 lần Drink và gần 4 lần Non-Consumable) phản ánh trọng tâm của tập dữ liệu chủ yếu tập trung vào các sản phẩm thực phẩm.
ggplot(long, aes(x = ProductFamily, fill=ProductFamily)) +
geom_bar(color="white") +
labs(title="Số giao dịch theo Product Family",
x="Nhóm sản phẩm", y="Số giao dịch") +
theme_minimal()
2.2.12 ProductDepartment
thong_ke_dinh_tinh(long,'ProductDepartment')## Gia_tri Tan_so Tan_suat
## 1 Alcoholic Beverages 356 2.53
## 2 Baked Goods 425 3.02
## 3 Baking Goods 1072 7.63
## 4 Beverages 680 4.84
## 5 Breakfast Foods 188 1.34
## 6 Canned Foods 977 6.95
## 7 Canned Products 109 0.78
## 8 Carousel 59 0.42
## 9 Checkout 82 0.58
## 10 Dairy 903 6.42
## 11 Deli 699 4.97
## 12 Eggs 198 1.41
## 13 Frozen Foods 1382 9.83
## 14 Health and Hygiene 893 6.35
## 15 Household 1420 10.10
## 16 Meat 89 0.63
## 17 Periodicals 202 1.44
## 18 Produce 1994 14.18
## 19 Seafood 102 0.73
## 20 Snack Foods 1600 11.38
## 21 Snacks 352 2.50
## 22 Starchy Foods 277 1.97Tập trung vào thực phẩm tươi và đồ gia dụng: Hai nhóm Produce và Household chiếm tổng cộng gần 25% dữ liệu, cho thấy đây có thể là các mặt hàng chủ lực.
Sự thiếu cân đối rõ rệt: Trong khi một số nhóm như Produce chiếm tới 14.18% thì nhiều nhóm khác như Meat, Seafood lại chiếm chưa tới 1%.
ggplot(long, aes(x = ProductDepartment)) +
geom_bar(fill = "yellow", color = "black") +
labs(x = "Sản phẩm", y = "Tần số", title = "Bộ phận sản phẩm khách hàng mua") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
2.2.13 ProductCategory
thong_ke_dinh_tinh(long,'ProductCategory')## Gia_tri Tan_so Tan_suat
## 1 Baking Goods 484 3.44
## 2 Bathroom Products 365 2.60
## 3 Beer and Wine 356 2.53
## 4 Bread 425 3.02
## 5 Breakfast Foods 417 2.97
## 6 Candles 45 0.32
## 7 Candy 352 2.50
## 8 Canned Anchovies 44 0.31
## 9 Canned Clams 53 0.38
## 10 Canned Oysters 35 0.25
## 11 Canned Sardines 40 0.28
## 12 Canned Shrimp 38 0.27
## 13 Canned Soup 404 2.87
## 14 Canned Tuna 87 0.62
## 15 Carbonated Beverages 154 1.10
## 16 Cleaning Supplies 189 1.34
## 17 Cold Remedies 93 0.66
## 18 Dairy 903 6.42
## 19 Decongestants 85 0.60
## 20 Drinks 135 0.96
## 21 Eggs 198 1.41
## 22 Electrical 355 2.53
## 23 Frozen Desserts 323 2.30
## 24 Frozen Entrees 118 0.84
## 25 Fruit 765 5.44
## 26 Hardware 129 0.92
## 27 Hot Beverages 226 1.61
## 28 Hygiene 197 1.40
## 29 Jams and Jellies 588 4.18
## 30 Kitchen Products 217 1.54
## 31 Magazines 202 1.44
## 32 Meat 761 5.41
## 33 Miscellaneous 42 0.30
## 34 Packaged Vegetables 48 0.34
## 35 Pain Relievers 192 1.37
## 36 Paper Products 345 2.45
## 37 Pizza 194 1.38
## 38 Plastic Products 141 1.00
## 39 Pure Juice Beverages 165 1.17
## 40 Seafood 102 0.73
## 41 Side Dishes 153 1.09
## 42 Snack Foods 1600 11.38
## 43 Specialty 289 2.06
## 44 Starchy Foods 277 1.97
## 45 Vegetables 1728 12.29Dữ liệu cho thấy Snack Foods (11.38%) và Vegetables (12.29%) là hai danh mục phổ biến nhất, chiếm tổng cộng gần 1/4 dữ liệu. Các danh mục đáng chú ý khác bao gồm Dairy (6.42%), Fruit (5.44%), và Meat (5.41%), trong khi nhiều danh mục như đồ hải sản đóng hộp (Canned Oysters 0.25%, Canned Shrimp 0.27%) hoặc các sản phẩm ít phổ biến (Candles 0.32%, Miscellaneous 0.30%) có tần suất rất thấp. Sự phân bố này phản ánh tập trung chính vào các mặt hàng tiêu dùng nhanh và thực phẩm tươi sống, trong khi các sản phẩm đặc biệt hoặc ít phổ biến chỉ chiếm tỷ lệ nhỏ.
ggplot(long, aes(x = ProductCategory)) +
geom_bar(fill = "yellow", color = "black") +
labs(x = "Sản phẩm", y = "Tần số", title = "Danh mục sản phẩm khách hàng mua") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
2.2.14 UnitsSold
skim(long$UnitsSold)| Name | long$UnitsSold |
| Number of rows | 14059 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 0 | 1 | 4.08 | 1.17 | 1 | 3 | 4 | 5 | 8 | ▁▃▇▂▁ |
Kết quả cho thấy trung bình 1 đơn mua sẽ mua khoảng 4.08 đơn vị sản phẩm, giá trị độ lệch chuẩn bằng 1.17 cho thấy mức độ phân tán tương đối thấp xung quanh giá trị trung bình.
Có thể thấy rằng có 50% giao dịch bán từ 3-5 đơn vị (khoảng tứ phân vị) và 75% giao dịch bán ≤5 đơn vị.
boxplot(long$UnitsSold, horizontal = TRUE, col = "yellow",
main = "Phân phối UnitsSold", xlab = "Số đơn vị bán ra")
2.2.15 Revenue
skim(long$Revenue)| Name | long$Revenue |
| Number of rows | 14059 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 0 | 1 | 13 | 8.22 | 0.53 | 6.84 | 11.25 | 17.37 | 56.7 | ▇▅▂▁▁ |
Biến Revenue (doanh thu) cho thấy một phân phối lệch phải rõ rệt với giá trị trung bình 13.00 và trung vị 11.25, phản ánh sự hiện diện của nhiều giao dịch có doanh thu cao kéo trung bình lên. Độ lệch chuẩn lớn (8.22) cùng khoảng biến thiên rộng (từ 0.53 đến 56.7) cho thấy mức độ chênh lệch đáng kể giữa các giao dịch.
boxplot(long$Revenue, horizontal = TRUE, col = "yellow",
main = "Phân phối Revenue", xlab = "Doanh thu")
Boxplot này 1 lần nữa xác nhận rõ ràng phân phối lệch phải mạnh của biến Revenue. Các điểm riêng lẻ xuất hiện ở vùng 30-50, trùng khớp với giá trị tối đa 56.7 trong thống kê. Khoảng cách xa từ p75 (~17) đến các điểm này cho thấy chúng thực sự là các giá trị cực trị