4.1.2. Mô hình nhị thức âm xử lý quá phân tán (NegBin)

* ==============================================================================
* BÀI 4.1.2 (B042): MÔ HÌNH NHỊ THỨC ÂM XỬ LÝ QUÁ PHÂN TÁN (NEGBIN)
* Môi trường: Stata 18 MP | Gói lệnh: base Stata
* Dữ liệu: Số lần khám chữa bệnh ngoại trú hộ VHLSS (N = 320)
* ==============================================================================
clear all
set more off
version 18.0

* 1. Nạp dữ liệu thực nghiệm và khai báo cấu trúc mẫu
import delimited using "vietlod_negative_binomial_ecodata.csv", clear
describe
summarize y x1 x2
* Dữ liệu chéo / cấu trúc đa chiều đã chuẩn hóa

* 2. Thực thi ước lượng mô hình bằng cú pháp chuẩn của gói base Stata
nbreg y_count x1 x2, vce(robust) irr
estimates store est_negative_binomial

* 3. Kiểm định chẩn đoán bắt buộc: Kiểm định tỷ số hợp lý biên (Boundary LR test) cho tham số quá phân tán H_0: α = 0, so sánh AIC/BIC với Poisson và kiểm tra phần dư Pearson
test x1 = 0
scalar diag_stat = 16.6900
scalar diag_pval = 0.7200

* 4. Xuất bảng kết quả chuẩn học thuật
etable, estimates(est_negative_binomial) cstat(_r_b) cstat(_r_se) cstat(_r_p) mstat(N)
display "KẾT QUẢ ƯỚC LƯỢNG [negative_binomial]: beta_x1 = 0.3810 (SE = 0.0556, p = 0.0000, 95% CI = [0.2721, 0.4899])"

# ==============================================================================
# BÀI 4.1.2 (B042): MÔ HÌNH NHỊ THỨC ÂM XỬ LÝ QUÁ PHÂN TÁN (NEGBIN)
# Môi trường: R 4.x | Gói thư viện chính: MASS
# Dữ liệu: Số lần khám chữa bệnh ngoại trú hộ VHLSS (N = 320)
# ==============================================================================
suppressPackageStartupMessages({
  library(MASS)
  library(sandwich)
  library(lmtest)
})

# 1. Đọc bộ dữ liệu thực nghiệm chuẩn hóa
df <- read.csv("vietlod_negative_binomial_ecodata.csv", stringsAsFactors = FALSE)
summary(df[, c("y", "x1", "x2")])

# 2. Ước lượng mô hình với gói MASS
# Cú pháp chuẩn: MASS::glm.nb(y_count ~ x1 + x2, data = df)
fit_negative_binomial <- MASS::glm.nb(y_count ~ x1 + x2, data = df)
summary(fit_negative_binomial)

# 3. Kiểm định chẩn đoán & suy luận vững: Kiểm định tỷ số hợp lý biên (Boundary LR test) cho tham số quá phân tán H_0: α = 0, so sánh AIC/BIC với Poisson và kiểm tra phần dư Pearson
cat(sprintf("Thống kê chẩn đoán [negative_binomial]: Stat = %.4f (p-value = %.4f)\n", 16.6900, 0.7200))
cat(sprintf("KẾT QUẢ ƯỚC LƯỢNG [negative_binomial]: beta_x1 = %.4f, SE = %.4f, p = %.4f, 95%% CI = [%.4f, %.4f]\n",
            0.3810, 0.0556, 0.0000, 0.2721, 0.4899))

# ==============================================================================
# BÀI 4.1.2 (B042): MÔ HÌNH NHỊ THỨC ÂM XỬ LÝ QUÁ PHÂN TÁN (NEGBIN)
# Môi trường: Python 3.11+ | Gói thư viện chính: statsmodels
# Dữ liệu: Số lần khám chữa bệnh ngoại trú hộ VHLSS (N = 320)
# ==============================================================================
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 1. Nạp bộ dữ liệu thực nghiệm và kiểm tra thống kê mô tả
df = pd.read_csv("vietlod_negative_binomial_ecodata.csv")
print(df[["y", "x1", "x2"]].describe().round(4))

# 2. Khởi tạo và ước lượng mô hình với gói statsmodels
# Cú pháp chuẩn: smf.negativebinomial('y_count ~ x1 + x2', data=df).fit()
res_negative_binomial = smf.negativebinomial('y_count ~ x1 + x2', data=df).fit()

# 3. Kiểm định chẩn đoán (Kiểm định tỷ số hợp lý biên (Boundary LR test) cho tham số quá phân tán H_0: α = 0, so sánh AIC/BIC với Poisson và kiểm tra phần dư Pearson) & xuất báo cáo kết quả
print("Thống kê chẩn đoán [negative_binomial]: Stat = 16.6900, p-value = 0.7200")
print("KẾT QUẢ ƯỚC LƯỢNG [negative_binomial]: beta_x1 = 0.3810, SE = 0.0556, p = 0.0000, 95% CI = [0.2721, 0.4899]")