knitr::opts_chunk$set(
echo = TRUE,
warning = FALSE,
message = FALSE
)Bài thực hành — Horseshoe Crab Count Data & Poisson Regression
1 1. Mục tiêu bài thực hành
Bài thực hành này sử dụng bộ dữ liệu Horseshoe Crab để phân tích dữ liệu biến đếm bằng hồi quy Poisson. Biến phụ thuộc chính là sat, thể hiện số lượng cua đực vệ tinh đi cùng một cua cái. Đây là biến đếm không âm nên phù hợp để minh họa hồi quy Poisson.
Các mục tiêu chính:
- Đọc dữ liệu online hoặc offline một cách ổn định.
- Kiểm tra cấu trúc dữ liệu và nhận diện biến phụ thuộc, biến độc lập.
- Thực hiện thống kê mô tả rõ ràng cho biến đếm và các biến giải thích.
- Giải thích ý nghĩa thống kê mô tả.
- Xây dựng mô hình hồi quy Poisson.
- Diễn giải hệ số bằng
exp(beta), tức Incident Rate Ratio — IRR. - Kiểm tra overdispersion và so sánh với Quasi-Poisson nếu cần.
- Dự báo số cua đực vệ tinh kỳ vọng cho một trường hợp cụ thể.
2 2. Thiết lập môi trường làm việc
needed_packages <- c("dplyr", "ggplot2", "knitr")
for (pkg in needed_packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
library(dplyr)
library(ggplot2)
library(knitr)3 3. Đọc dữ liệu Horseshoe Crab
Dữ liệu online được lấy từ file Crabs.dat. Vì file raw có thể được lưu dưới dạng nhiều khoảng trắng thay vì bảng CSV chuẩn, ta viết một hàm đọc dữ liệu an toàn bằng cách tách toàn bộ nội dung thành token.
online_url <- "https://raw.githubusercontent.com/alanagresti/categorical-data/master/Crabs.dat"
offline_path <- "../../Data/Crabs.dat"
read_crabs_data <- function(path_or_url) {
raw_lines <- readLines(path_or_url, warn = FALSE)
tokens <- unlist(strsplit(paste(raw_lines, collapse = " "), "\\s+"))
tokens <- tokens[tokens != ""]
header <- tokens[1:7]
values <- tokens[-(1:7)]
if (length(values) %% length(header) != 0) {
stop("Số lượng giá trị không chia hết cho số cột. Cần kiểm tra lại định dạng dữ liệu.")
}
mat <- matrix(values, ncol = length(header), byrow = TRUE)
df <- as.data.frame(mat, stringsAsFactors = FALSE)
names(df) <- header
df[] <- lapply(df, type.convert, as.is = TRUE)
return(df)
}
if (file.exists(offline_path)) {
df <- read_crabs_data(offline_path)
} else {
df <- read_crabs_data(online_url)
}
dim(df)[1] 173 7head(df) crab sat y weight width color spine
1 1 8 1 3.05 28.3 2 3
2 2 0 0 1.55 22.5 3 3
3 3 9 1 2.30 26.0 1 1
4 4 0 0 2.10 24.8 3 3
5 5 4 1 2.60 26.0 3 3
6 6 0 0 2.10 23.8 2 3Giải thích:
Đoạn code trên ưu tiên đọc file offline tại ../../Data/Crabs.dat. Nếu file chưa có trong máy, R sẽ đọc trực tiếp từ đường dẫn online. Cách đọc này giúp hạn chế lỗi khi dữ liệu không phải dạng CSV thông thường.
4 4. Kiểm tra cấu trúc dữ liệu
str(df)'data.frame': 173 obs. of 7 variables:
$ crab : int 1 2 3 4 5 6 7 8 9 10 ...
$ sat : int 8 0 9 0 4 0 0 0 0 0 ...
$ y : int 1 0 1 0 1 0 0 0 0 0 ...
$ weight: num 3.05 1.55 2.3 2.1 2.6 2.1 2.35 1.9 1.95 2.15 ...
$ width : num 28.3 22.5 26 24.8 26 23.8 26.5 24.7 23.7 25.6 ...
$ color : int 2 3 1 3 3 2 1 3 2 3 ...
$ spine : int 3 3 1 3 3 3 1 2 1 3 ...summary(df) crab sat y weight width
Min. : 1 Min. : 0.000 Min. :0.0000 Min. :1.200 Min. :21.0
1st Qu.: 44 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:24.9
Median : 87 Median : 2.000 Median :1.0000 Median :2.350 Median :26.1
Mean : 87 Mean : 2.919 Mean :0.6416 Mean :2.437 Mean :26.3
3rd Qu.:130 3rd Qu.: 5.000 3rd Qu.:1.0000 3rd Qu.:2.850 3rd Qu.:27.7
Max. :173 Max. :15.000 Max. :1.0000 Max. :5.200 Max. :33.5
color spine
Min. :1.000 Min. :1.000
1st Qu.:2.000 1st Qu.:2.000
Median :2.000 Median :3.000
Mean :2.439 Mean :2.486
3rd Qu.:3.000 3rd Qu.:3.000
Max. :4.000 Max. :3.000 missing_table <- data.frame(
variable = names(df),
missing_count = colSums(is.na(df)),
missing_percent = round(colMeans(is.na(df)) * 100, 2)
)
kable(missing_table, caption = "Bảng kiểm tra giá trị khuyết thiếu")| variable | missing_count | missing_percent | |
|---|---|---|---|
| crab | crab | 0 | 0 |
| sat | sat | 0 | 0 |
| y | y | 0 | 0 |
| weight | weight | 0 | 0 |
| width | width | 0 | 0 |
| color | color | 0 | 0 |
| spine | spine | 0 | 0 |
Giải thích:
Bước này giúp kiểm tra dữ liệu có được đọc đúng hay không. Ta cần bảo đảm các biến như sat, weight, width, color, spine được nhận diện đúng kiểu số. Nếu có giá trị khuyết thiếu, cần xử lý trước khi xây dựng mô hình. Với hồi quy Poisson, biến phụ thuộc sat phải là biến đếm không âm.
5 5. Nhận diện biến nghiên cứu
Bộ dữ liệu có các biến chính sau:
| Biến | Ý nghĩa | Vai trò trong bài |
|---|---|---|
crab |
Mã số quan sát | Biến định danh, không đưa vào mô hình |
sat |
Số cua đực vệ tinh | Biến phụ thuộc dạng biến đếm |
y |
Có ít nhất một cua vệ tinh hay không | Biến nhị phân, không dùng làm biến phụ thuộc trong Poisson |
weight |
Cân nặng cua | Biến độc lập định lượng |
width |
Chiều rộng mai cua | Biến độc lập định lượng |
color |
Nhóm màu sắc | Biến độc lập phân loại |
spine |
Nhóm tình trạng gai/spine | Biến độc lập phân loại |
Lưu ý quan trọng:
Trong mô hình Poisson bên dưới, ta không đưa y vào mô hình vì y gần như được suy ra từ sat: nếu sat > 0 thì y = 1, còn nếu sat = 0 thì y = 0. Nếu dùng y để giải thích sat, ta đang đưa một biến chứa trực tiếp thông tin của biến phụ thuộc vào mô hình, làm sai ý nghĩa phân tích.
6 6. Tiền xử lý biến phân loại
df <- df %>%
mutate(
color_f = factor(color),
spine_f = factor(spine),
y_f = factor(y, levels = c(0, 1), labels = c("No satellite", "Has satellite"))
)
str(df)'data.frame': 173 obs. of 10 variables:
$ crab : int 1 2 3 4 5 6 7 8 9 10 ...
$ sat : int 8 0 9 0 4 0 0 0 0 0 ...
$ y : int 1 0 1 0 1 0 0 0 0 0 ...
$ weight : num 3.05 1.55 2.3 2.1 2.6 2.1 2.35 1.9 1.95 2.15 ...
$ width : num 28.3 22.5 26 24.8 26 23.8 26.5 24.7 23.7 25.6 ...
$ color : int 2 3 1 3 3 2 1 3 2 3 ...
$ spine : int 3 3 1 3 3 3 1 2 1 3 ...
$ color_f: Factor w/ 4 levels "1","2","3","4": 2 3 1 3 3 2 1 3 2 3 ...
$ spine_f: Factor w/ 3 levels "1","2","3": 3 3 1 3 3 3 1 2 1 3 ...
$ y_f : Factor w/ 2 levels "No satellite",..: 2 1 2 1 2 1 1 1 1 1 ...Giải thích:
Các biến color và spine là biến phân loại được mã hóa bằng số. Nếu để nguyên dạng số, mô hình có thể hiểu nhầm rằng khoảng cách giữa các mức là tuyến tính. Vì vậy, ta chuyển chúng thành factor để mô hình ước lượng tác động riêng cho từng nhóm so với nhóm tham chiếu.
7 7. Thống kê mô tả cho biến phụ thuộc sat
sat_desc <- df %>%
summarise(
n = n(),
min = min(sat),
q1 = quantile(sat, 0.25),
mean = mean(sat),
median = median(sat),
q3 = quantile(sat, 0.75),
max = max(sat),
variance = var(sat),
sd = sd(sat),
zero_count = sum(sat == 0),
zero_percent = mean(sat == 0) * 100
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
kable(sat_desc, caption = "Thống kê mô tả biến đếm sat")| n | min | q1 | mean | median | q3 | max | variance | sd | zero_count | zero_percent |
|---|---|---|---|---|---|---|---|---|---|---|
| 173 | 0 | 0 | 2.919 | 2 | 5 | 15 | 9.912 | 3.148 | 62 | 35.838 |
Giải thích thống kê mô tả:
Biến sat là số lượng cua đực vệ tinh nên chỉ nhận giá trị nguyên không âm. Trung bình cho biết số cua vệ tinh kỳ vọng thô trên mỗi quan sát. Trung vị giúp đánh giá xu hướng trung tâm khi phân phối bị lệch. Nếu trung bình lớn hơn trung vị và giá trị lớn nhất khá cao, dữ liệu có dấu hiệu lệch phải — đặc trưng phổ biến của dữ liệu đếm. Tỷ lệ quan sát bằng 0 cũng rất quan trọng vì nếu quá nhiều giá trị 0, ta có thể cần cân nhắc thêm mô hình Zero-Inflated bên cạnh Poisson thông thường.
8 8. So sánh mean và variance của biến đếm
mean_sat <- mean(df$sat)
var_sat <- var(df$sat)
mean_var_table <- data.frame(
statistic = c("Mean of sat", "Variance of sat", "Variance / Mean"),
value = c(mean_sat, var_sat, var_sat / mean_sat)
)
mean_var_table$value <- round(mean_var_table$value, 4)
kable(mean_var_table, caption = "So sánh trung bình và phương sai của biến sat")| statistic | value |
|---|---|
| Mean of sat | 2.9191 |
| Variance of sat | 9.9120 |
| Variance / Mean | 3.3956 |
Giải thích:
Phân phối Poisson có giả định quan trọng là trung bình xấp xỉ phương sai, tức Mean = Variance. Nếu phương sai lớn hơn trung bình rõ rệt, dữ liệu có thể có overdispersion. Tuy nhiên, so sánh mean và variance ở đây mới chỉ là kiểm tra mô tả thô của biến sat, chưa xét đến các biến giải thích như width, weight, color và spine. Vì vậy, kết luận chính thức về overdispersion nên dựa thêm vào kiểm tra sau khi đã ước lượng mô hình.
9 9. Bảng tần số biến đếm sat
sat_freq <- df %>%
count(sat, name = "frequency") %>%
mutate(percent = round(frequency / sum(frequency) * 100, 2))
kable(sat_freq, caption = "Bảng tần số số cua đực vệ tinh")| sat | frequency | percent |
|---|---|---|
| 0 | 62 | 35.84 |
| 1 | 16 | 9.25 |
| 2 | 9 | 5.20 |
| 3 | 19 | 10.98 |
| 4 | 19 | 10.98 |
| 5 | 15 | 8.67 |
| 6 | 13 | 7.51 |
| 7 | 4 | 2.31 |
| 8 | 6 | 3.47 |
| 9 | 3 | 1.73 |
| 10 | 3 | 1.73 |
| 11 | 1 | 0.58 |
| 12 | 1 | 0.58 |
| 14 | 1 | 0.58 |
| 15 | 1 | 0.58 |
Giải thích:
Bảng tần số cho biết có bao nhiêu con cua cái có 0, 1, 2, 3,… cua đực vệ tinh. Với dữ liệu biến đếm, phần lớn quan sát thường tập trung ở các giá trị nhỏ, đặc biệt là giá trị 0. Nếu bảng cho thấy số lượng giá trị 0 chiếm tỷ lệ lớn, đây là tín hiệu cần chú ý khi đánh giá độ phù hợp của mô hình Poisson.
10 10. Trực quan hóa phân phối biến đếm
ggplot(df, aes(x = sat)) +
geom_histogram(binwidth = 1, boundary = -0.5) +
labs(
title = "Phân phối số cua đực vệ tinh",
x = "Số cua đực vệ tinh sat",
y = "Tần số"
) +
theme_minimal()
Nhận xét:
Biểu đồ histogram giúp quan sát trực tiếp hình dạng phân phối của sat. Nếu phân phối lệch phải, có nhiều quan sát tại 0 và một số ít quan sát có giá trị rất lớn, điều này phù hợp với bản chất của dữ liệu đếm trong sinh học. Đây cũng là lý do hồi quy Poisson thường phù hợp hơn OLS, vì OLS không được thiết kế riêng cho biến phụ thuộc rời rạc và không âm.
11 11. Thống kê mô tả cho các biến định lượng width và weight
numeric_desc <- df %>%
summarise(
width_mean = mean(width),
width_sd = sd(width),
width_min = min(width),
width_median = median(width),
width_max = max(width),
weight_mean = mean(weight),
weight_sd = sd(weight),
weight_min = min(weight),
weight_median = median(weight),
weight_max = max(weight)
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
kable(numeric_desc, caption = "Thống kê mô tả width và weight")| width_mean | width_sd | width_min | width_median | width_max | weight_mean | weight_sd | weight_min | weight_median | weight_max |
|---|---|---|---|---|---|---|---|---|---|
| 26.299 | 2.109 | 21 | 26.1 | 33.5 | 2.437 | 0.577 | 1.2 | 2.35 | 5.2 |
cor_width_weight <- cor(df$width, df$weight)
round(cor_width_weight, 4)[1] 0.8869Giải thích:
width và weight đều phản ánh kích thước cơ thể của cua. Nếu hai biến này có tương quan cao, khi đưa cả hai vào cùng một mô hình, cần diễn giải hệ số cẩn thận vì chúng có thể cùng đại diện cho một đặc điểm sinh học chung là kích thước. Trong mô hình hồi quy, hệ số của width được hiểu là tác động của chiều rộng khi giữ cân nặng và các biến khác không đổi; ngược lại, hệ số của weight được hiểu là tác động của cân nặng khi giữ chiều rộng và các biến khác không đổi.
12 12. Thống kê mô tả theo nhóm màu sắc và spine
color_desc <- df %>%
group_by(color_f) %>%
summarise(
n = n(),
mean_sat = mean(sat),
variance_sat = var(sat),
mean_width = mean(width),
mean_weight = mean(weight),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
kable(color_desc, caption = "Thống kê mô tả theo nhóm color")| color_f | n | mean_sat | variance_sat | mean_width | mean_weight |
|---|---|---|---|---|---|
| 1 | 12 | 4.083 | 9.720 | 26.967 | 2.629 |
| 2 | 95 | 3.295 | 10.274 | 26.704 | 2.538 |
| 3 | 44 | 2.227 | 6.738 | 25.750 | 2.299 |
| 4 | 22 | 2.045 | 13.093 | 25.282 | 2.174 |
spine_desc <- df %>%
group_by(spine_f) %>%
summarise(
n = n(),
mean_sat = mean(sat),
variance_sat = var(sat),
mean_width = mean(width),
mean_weight = mean(weight),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
kable(spine_desc, caption = "Thống kê mô tả theo nhóm spine")| spine_f | n | mean_sat | variance_sat | mean_width | mean_weight |
|---|---|---|---|---|---|
| 1 | 37 | 3.649 | 11.512 | 27.111 | 2.679 |
| 2 | 15 | 2.000 | 5.571 | 24.727 | 2.153 |
| 3 | 121 | 2.810 | 9.822 | 26.245 | 2.398 |
Giải thích:
Các bảng theo nhóm giúp xem số lượng cua đực vệ tinh trung bình có khác nhau giữa các nhóm màu sắc hoặc nhóm spine hay không. Nếu một nhóm có mean_sat cao hơn rõ rệt, điều đó gợi ý nhóm này có thể liên hệ với số cua vệ tinh nhiều hơn. Tuy nhiên, đây mới là so sánh mô tả, chưa kiểm soát đồng thời các biến khác như width và weight. Vì vậy, cần hồi quy Poisson để phân tích mối liên hệ sau khi giữ các biến khác không đổi.
13 13. Trực quan hóa mối liên hệ giữa sat và width
ggplot(df, aes(x = width, y = sat)) +
geom_point(alpha = 0.65) +
geom_smooth(method = "glm", method.args = list(family = poisson), se = TRUE) +
labs(
title = "Mối liên hệ giữa chiều rộng mai cua và số cua đực vệ tinh",
x = "Width",
y = "Số cua đực vệ tinh sat"
) +
theme_minimal()
Nhận xét:
Đồ thị giúp quan sát trực quan liệu cua cái có kích thước lớn hơn có xu hướng có nhiều cua đực vệ tinh hơn hay không. Vì sat là biến đếm, đường xu hướng được vẽ theo mô hình Poisson thay vì đường hồi quy tuyến tính OLS. Nếu đường xu hướng đi lên, điều này gợi ý mối liên hệ dương giữa kích thước cua và số cua vệ tinh.
14 14. Xây dựng mô hình hồi quy Poisson
Mô hình Poisson được viết dưới dạng:
[ (_i) = _0 + _1 width_i + _2 weight_i + _3 color_i + _4 spine_i ]
Trong đó:
[ _i = E(sat_i X_i) ]
Tức là lambda là số cua đực vệ tinh kỳ vọng của quan sát thứ i.
model_p <- glm(
sat ~ width + weight + color_f + spine_f,
data = df,
family = poisson(link = "log")
)
summary(model_p)
Call:
glm(formula = sat ~ width + weight + color_f + spine_f, family = poisson(link = "log"),
data = df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.36180 0.96655 -0.374 0.70817
width 0.01675 0.04892 0.342 0.73207
weight 0.49647 0.16626 2.986 0.00283 **
color_f2 -0.26485 0.16811 -1.575 0.11515
color_f3 -0.51371 0.19536 -2.629 0.00855 **
color_f4 -0.53086 0.22692 -2.339 0.01931 *
spine_f2 -0.15037 0.21358 -0.704 0.48139
spine_f3 0.08728 0.11993 0.728 0.46674
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 632.79 on 172 degrees of freedom
Residual deviance: 549.59 on 165 degrees of freedom
AIC: 920.88
Number of Fisher Scoring iterations: 6Giải thích mô hình:
Mô hình trên sử dụng sat làm biến phụ thuộc. Các biến width, weight, color_f và spine_f được dùng làm biến độc lập. Vì Poisson dùng log link, hệ số trong bảng summary() nằm trên thang log(lambda), không phải trên thang số đếm gốc. Do đó, không nên đọc hệ số Poisson theo kiểu “tăng thêm beta đơn vị” như OLS. Muốn diễn giải trên thang dễ hiểu hơn, cần lấy exp(beta) để chuyển sang IRR.
15 15. Chuyển hệ số sang IRR
coef_mat <- summary(model_p)$coefficients
irr_table <- data.frame(
variable = rownames(coef_mat),
beta = coef_mat[, 1],
std_error = coef_mat[, 2],
z_value = coef_mat[, 3],
p_value = coef_mat[, 4],
IRR = exp(coef_mat[, 1]),
CI_95_low = exp(coef_mat[, 1] - 1.96 * coef_mat[, 2]),
CI_95_high = exp(coef_mat[, 1] + 1.96 * coef_mat[, 2])
)
irr_table <- irr_table %>%
mutate(across(where(is.numeric), ~ round(.x, 4)))
kable(irr_table, caption = "Bảng hệ số Poisson và IRR")| variable | beta | std_error | z_value | p_value | IRR | CI_95_low | CI_95_high | |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | (Intercept) | -0.3618 | 0.9666 | -0.3743 | 0.7082 | 0.6964 | 0.1047 | 4.6304 |
| width | width | 0.0167 | 0.0489 | 0.3424 | 0.7321 | 1.0169 | 0.9239 | 1.1192 |
| weight | weight | 0.4965 | 0.1663 | 2.9861 | 0.0028 | 1.6429 | 1.1860 | 2.2758 |
| color_f2 | color_f2 | -0.2649 | 0.1681 | -1.5755 | 0.1152 | 0.7673 | 0.5519 | 1.0668 |
| color_f3 | color_f3 | -0.5137 | 0.1954 | -2.6295 | 0.0086 | 0.5983 | 0.4079 | 0.8774 |
| color_f4 | color_f4 | -0.5309 | 0.2269 | -2.3395 | 0.0193 | 0.5881 | 0.3770 | 0.9175 |
| spine_f2 | spine_f2 | -0.1504 | 0.2136 | -0.7041 | 0.4814 | 0.8604 | 0.5661 | 1.3077 |
| spine_f3 | spine_f3 | 0.0873 | 0.1199 | 0.7278 | 0.4667 | 1.0912 | 0.8626 | 1.3804 |
Cách diễn giải IRR:
IRR = exp(beta) cho biết tỷ lệ thay đổi của số đếm kỳ vọng khi biến độc lập tăng thêm một đơn vị hoặc khi chuyển từ nhóm tham chiếu sang nhóm đang xét, giữ các biến khác không đổi. Nếu IRR > 1, số cua vệ tinh kỳ vọng tăng. Nếu IRR < 1, số cua vệ tinh kỳ vọng giảm. Nếu p-value nhỏ hơn 0.05, ta có bằng chứng thống kê rõ hơn về mối liên hệ giữa biến đó và số cua vệ tinh kỳ vọng.
Ví dụ, nếu hệ số của width có IRR lớn hơn 1 và có ý nghĩa thống kê, có thể diễn giải rằng cua cái có chiều rộng mai lớn hơn sẽ có số cua đực vệ tinh kỳ vọng cao hơn, sau khi đã kiểm soát cân nặng, màu sắc và spine.
16 16. Kiểm tra overdispersion
residual_deviance_ratio <- deviance(model_p) / df.residual(model_p)
pearson_chisq_ratio <- sum(residuals(model_p, type = "pearson")^2) / df.residual(model_p)
od_table <- data.frame(
metric = c(
"Residual deviance / df",
"Pearson chi-square / df",
"Mean of sat",
"Variance of sat",
"Variance / Mean"
),
value = c(
residual_deviance_ratio,
pearson_chisq_ratio,
mean(df$sat),
var(df$sat),
var(df$sat) / mean(df$sat)
)
)
od_table$value <- round(od_table$value, 4)
kable(od_table, caption = "Kiểm tra overdispersion")| metric | value |
|---|---|
| Residual deviance / df | 3.3308 |
| Pearson chi-square / df | 3.2353 |
| Mean of sat | 2.9191 |
| Variance of sat | 9.9120 |
| Variance / Mean | 3.3956 |
if (!requireNamespace("AER", quietly = TRUE)) {
cat("Package AER chưa được cài. Nếu muốn chạy dispersiontest(), hãy chạy install.packages('AER').")
} else {
AER::dispersiontest(model_p)
}
Overdispersion test
data: model_p
z = 5.1494, p-value = 1.307e-07
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
3.104872 Giải thích:
Trong mô hình Poisson, một giả định quan trọng là phương sai xấp xỉ trung bình. Nếu Residual deviance / df hoặc Pearson chi-square / df lớn hơn 1 rõ rệt, đặc biệt là lớn hơn khoảng 1.5 hoặc 2, ta cần nghi ngờ overdispersion. Overdispersion làm sai số chuẩn của mô hình Poisson bị đánh giá thấp, khiến p-value có thể trở nên quá lạc quan. Khi đó, Quasi-Poisson hoặc Negative Binomial có thể phù hợp hơn.
17 17. So sánh với mô hình Quasi-Poisson
model_qp <- glm(
sat ~ width + weight + color_f + spine_f,
data = df,
family = quasipoisson(link = "log")
)
summary(model_qp)
Call:
glm(formula = sat ~ width + weight + color_f + spine_f, family = quasipoisson(link = "log"),
data = df)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.36180 1.73852 -0.208 0.8354
width 0.01675 0.08799 0.190 0.8493
weight 0.49647 0.29905 1.660 0.0988 .
color_f2 -0.26485 0.30238 -0.876 0.3824
color_f3 -0.51371 0.35139 -1.462 0.1457
color_f4 -0.53086 0.40815 -1.301 0.1952
spine_f2 -0.15037 0.38415 -0.391 0.6960
spine_f3 0.08728 0.21571 0.405 0.6863
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 3.235255)
Null deviance: 632.79 on 172 degrees of freedom
Residual deviance: 549.59 on 165 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 6se_compare <- data.frame(
variable = names(coef(model_p)),
SE_Poisson = summary(model_p)$coefficients[, 2],
SE_QuasiPoisson = summary(model_qp)$coefficients[, 2]
)
se_compare <- se_compare %>%
mutate(across(where(is.numeric), ~ round(.x, 4)))
kable(se_compare, caption = "So sánh sai số chuẩn Poisson và Quasi-Poisson")| variable | SE_Poisson | SE_QuasiPoisson | |
|---|---|---|---|
| (Intercept) | (Intercept) | 0.9666 | 1.7385 |
| width | width | 0.0489 | 0.0880 |
| weight | weight | 0.1663 | 0.2990 |
| color_f2 | color_f2 | 0.1681 | 0.3024 |
| color_f3 | color_f3 | 0.1954 | 0.3514 |
| color_f4 | color_f4 | 0.2269 | 0.4081 |
| spine_f2 | spine_f2 | 0.2136 | 0.3842 |
| spine_f3 | spine_f3 | 0.1199 | 0.2157 |
Giải thích:
Quasi-Poisson thường không làm thay đổi hệ số ước lượng chính, nhưng điều chỉnh sai số chuẩn theo mức độ phân tán của dữ liệu. Nếu dữ liệu có overdispersion, sai số chuẩn của Quasi-Poisson thường lớn hơn Poisson, làm kiểm định thống kê thận trọng hơn. Nếu kết luận giữa Poisson và Quasi-Poisson không thay đổi nhiều, kết quả có tính ổn định cao hơn.
18 18. Dự báo số cua đực vệ tinh kỳ vọng
Ta dự báo cho một cua cái điển hình có width và weight bằng giá trị trung bình, color_f và spine_f bằng nhóm xuất hiện nhiều nhất trong dữ liệu.
most_common_color <- names(which.max(table(df$color_f)))
most_common_spine <- names(which.max(table(df$spine_f)))
new_crab <- data.frame(
width = mean(df$width),
weight = mean(df$weight),
color_f = factor(most_common_color, levels = levels(df$color_f)),
spine_f = factor(most_common_spine, levels = levels(df$spine_f))
)
new_crab width weight color_f spine_f
1 26.29884 2.437191 2 3pred_link <- predict(model_p, newdata = new_crab, type = "link")
pred_response <- predict(model_p, newdata = new_crab, type = "response")
prediction_table <- data.frame(
predicted_log_lambda = as.numeric(pred_link),
predicted_lambda = as.numeric(pred_response)
)
prediction_table <- prediction_table %>%
mutate(across(where(is.numeric), ~ round(.x, 4)))
kable(prediction_table, caption = "Dự báo trên thang log(lambda) và lambda")| predicted_log_lambda | predicted_lambda |
|---|---|
| 1.1111 | 3.0377 |
Giải thích:
Với mô hình Poisson, predict(..., type = "link") trả về giá trị trên thang tuyến tính, tức log(lambda). Ngược lại, predict(..., type = "response") trả về lambda, tức số cua đực vệ tinh kỳ vọng trên thang gốc của biến đếm. Khi báo cáo kết quả dự báo cho người đọc, nên dùng type = "response" vì kết quả này có ý nghĩa trực tiếp hơn.
19 19. Vẽ đường dự báo theo width và color
pred_grid <- expand.grid(
width = seq(min(df$width), max(df$width), length.out = 100),
weight = mean(df$weight),
color_f = levels(df$color_f),
spine_f = factor(most_common_spine, levels = levels(df$spine_f))
)
pred_grid$predicted_sat <- predict(model_p, newdata = pred_grid, type = "response")
ggplot(pred_grid, aes(x = width, y = predicted_sat, color = color_f)) +
geom_line(linewidth = 1.1) +
labs(
title = "Số cua đực vệ tinh dự báo theo width và color",
x = "Width",
y = "Số cua đực vệ tinh kỳ vọng",
color = "Color group"
) +
theme_minimal()
Nhận xét:
Đồ thị trên biểu diễn số cua đực vệ tinh kỳ vọng theo chiều rộng mai cua, trong khi giữ weight ở mức trung bình và spine_f ở nhóm phổ biến nhất. Nếu các đường dự báo đi lên khi width tăng, điều đó cho thấy cua cái có chiều rộng lớn hơn có xu hướng có nhiều cua đực vệ tinh hơn. Nếu các đường theo màu sắc cách xa nhau, màu sắc có thể liên hệ đến số cua vệ tinh kỳ vọng sau khi kiểm soát các biến còn lại.
20 20. Báo cáo kết quả phân tích
Dữ liệu Horseshoe Crab là một ví dụ phù hợp cho hồi quy Poisson vì biến phụ thuộc sat là biến đếm không âm, thể hiện số cua đực vệ tinh đi cùng một cua cái. Thống kê mô tả ban đầu giúp nhận diện phân phối của sat, mức độ tập trung ở giá trị 0, cũng như sự khác biệt sơ bộ giữa các nhóm color và spine. Việc so sánh mean và variance của sat là bước kiểm tra ban đầu để xem dữ liệu có dấu hiệu overdispersion hay không, nhưng đây chưa phải kết luận cuối cùng vì chưa kiểm soát các biến giải thích.
Mô hình hồi quy Poisson được xây dựng với các biến độc lập width, weight, color_f và spine_f. Do mô hình dùng log link, hệ số hồi quy ban đầu nằm trên thang log(lambda). Vì vậy, khi diễn giải kết quả cần chuyển hệ số sang exp(beta), tức IRR. IRR cho biết số cua đực vệ tinh kỳ vọng thay đổi theo tỷ lệ bao nhiêu lần khi biến độc lập tăng thêm một đơn vị hoặc khi chuyển từ nhóm tham chiếu sang nhóm khác, giữ các biến còn lại không đổi.
Khi kiểm tra overdispersion, cần xem cả Residual deviance / df và Pearson chi-square / df. Nếu các tỷ số này lớn hơn 1 rõ rệt, mô hình Poisson có thể đánh giá thấp sai số chuẩn, từ đó làm p-value quá nhỏ. Trong trường hợp đó, Quasi-Poisson là lựa chọn hợp lý để kiểm tra độ nhạy của kết luận. Cuối cùng, kết quả hồi quy chỉ nên được hiểu là mối liên hệ thống kê giữa kích thước, màu sắc, spine và số cua vệ tinh kỳ vọng; không nên diễn giải trực tiếp thành quan hệ nhân quả tuyệt đối nếu không có thiết kế nghiên cứu phù hợp.