GIÁO TRÌNH TRANG 66 - 70

library(kableExtra)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(gridExtra)

Khi đó các giá trị phù hợp cho bốn loại ngáy ngủ sẽ thay đổi một chút. Nhưng chúng sẽ không thay đổi nếu khoảng cách tương đối giữa các điểm số là như nhau, chẳng hạn như {0, 4, 8, 10} hoặc {1, 3, 5, 6}. Cho dù vậy, đối với những dữ liệu này, bất cứ điểm số nào tăng đều cho kết luận rằng nguy cơ mắc bệnh tim tăng lên khi mức ngáy tăng lên.

Nếu nhập dữ liệu là 2484 quan sát nhị phân của 0 hoặc 1 và ước lượng mô hình xác suất tuyến tính bằng OLS:

\[\begin{array}{lrllr} \hline \text{Dependent Variable}: \text{𝜋} \\ \text{Method: Least Squares} \\ \text{Included observations}: 2484 & & & & \\ \hline \text{Variable} & \text{Coefficient} & \text{Std. Error} & \text{z-Statistic} & \text{Prob.} \\ \hline\hline \text{C} & 0.016872 & 0.020038 & 3.271697 & 0.0011 \\ \text{X} & 0.020038 & 0.002317 & ~~~8.649795 & 0.0000 \\ \hline\hline \text{R-squared} & 0.029263 & \text{Mean dependent var} & & 0.044283 & \\ \text{S.E. of regression} & 0.202773 & \text{Akaike info criterion} & & -0.352651 & \\ \text{Log likelihood} & 439.9923 & \text{Hannan-Quinn criter} & & -0.350950 & \\ \text{Prob(F-statistic)} & 0.000000 & & & \\ \end{array}\]

Bảng 3.3. Kết quả hồi quy mô hình xác suất tuyến tính cho dữ liệu bảng 3.1 bằng OLS

Chúng ta nhận được mô hình xác suất tuyến tính ước lượng:

\[ \hat{\pi} = 0,016872+0,020038.x \] Trong thực tế, khi mẫu phù hợp, ước lượng OLS và ước lượng ML thường gần nhau.

Mô hình hồi quy logistic

Như trên chúng ta đã thấy sự hạn chế của mô hình xác suất tuyến tính. Sự thay đổi của xác suất thành công tuy phụ thuộc đơn điệu theo X, nhưng thường là phi tuyến. Các đường cong hình chữ S thể hiện trong hình 3.2 thường là hình dạng thực tế cho mối quan hệ. Để khắc phục các hạn chế trên của LPM, người ta muốn tìm một hàm liên kết mà nó đồng biến với xác suất 𝜋 và giá trị biến thiên của nó là tùy ý (từ \(-\infty\) đến \(+\infty\) ). Hàm liên kết được đề xuất là \(\log(\frac{\pi(x)}{1-\pi(x)})\), được gọi là hàm \(logit(\pi)\) và ký hiệu là và mô hình hồi quy tương ứng là:

\[\begin{equation*} \log(\frac{\pi(x)}{1-\pi(x)})= \alpha + \beta x \tag{3.2.2} \end{equation*}\]

Đây được gọi là mô hình logistic hay mô hình logit. Từ mô hình, ta có: \[\pi(x) = \frac{1}{1 + e^{-\alpha + \beta x}}\]

Mô hình hồi quy logistic là trường hợp đặc biệt của GLM. Thành phần ngẫu nhiên (thành công, thất bại) là nhị phân. Hàm liên kết logit của \(\pi\) là phép biến đổi \(\log[\frac{\pi}{1 - \pi}]\), ký hiệu bởi \(logit(\pi)\) . Các mô hình hồi quy logistic thường được gọi là mô hình logit. Logit là tham số tự nhiên của phân phối nhị thức. Trong đó \(\pi\) bị giới hạn trong phạm vi (0, 1), logit có thể là số thực bất kỳ. Tập số thực cũng là khoảng giá trị có thể có cho các dự báo tuyến tính \((\alpha + \beta x)\), do đó, mô hình này không có vấn đề cấu trúc mà mô hình xác suất tuyến tính mắc phải.

# Hàm logistic
logistic <- function(x, beta0, beta1) {
  1 / (1 + exp(-(beta0 + beta1 * x)))
}

# Tạo dữ liệu cho biến độc lập X
x <- seq(-10, 10, by = 0.1)

# Hệ số cho mô hình logistic
beta0 <- 0
beta1_positive <- 1
beta1_negative <- -1

# Tính xác suất P(Y=1) cho mỗi giá trị của X
y_positive <- logistic(x, beta0, beta1_positive)
y_negative <- logistic(x, beta0, beta1_negative)

# Tạo data frame để dùng ggplot2
data_positive <- data.frame(x = x, y = y_positive)
data_negative <- data.frame(x = x, y = y_negative)

# Vẽ đồ thị
plot1 <- ggplot(data_positive, aes(x = x, y = y)) +
  geom_line() +
  labs(title = expression(beta > 0), y = "Y", x = "X") +
  theme_minimal() +
  ylim(0, 1)

plot2 <- ggplot(data_negative, aes(x = x, y = y)) +
  geom_line() +
  labs(title = expression(beta < 0), y = "Y", x = "X") +
  theme_minimal() +
  ylim(0, 1)

grid.arrange(plot1, plot2, ncol = 2)

Hình 3.2 (a,b)

Tham số \(\beta\) trong (3.2.2) là tốc độ biến thiên của logit nhưng logit đồng biến theo \(\pi(x)\) nên khi nó cũng thể hiện xu hướng biến thiên của \(\pi(x)\) theo x. Khi \(\beta>0\) , \(\pi(x)\) tăng khi tăng, như trong hình 3.2a. Khi \(\beta<0\), \(\pi(x)\) giảm khi x tăng, như trong hình 3.2b. Độ lớn của xác định tốc độ tăng hoặc giảm đường cong. Khi \(|\beta|\) tăng, đường cong có tỷ lệ thay đổi dốc hơn. Khi \(\beta=0\) , đường cong phẳng gần như đường thẳng nằm ngang.

Phương pháp ML để tìm ước lượng \(\hat{\alpha},\hat{\beta}\) cho các hệ số \(\alpha,\beta\) là trên cơ sở mẫu quan sát: {\((x_1,y_1),(x_2,y_2),...,(x_n,y_n)\)}, tìm nghiệm của phương trình hợp lý: \[\begin{align}\begin{cases}\frac{\partial log(L(\alpha,\beta))}{\partial\alpha}=0\\\frac{\partial log(L(\alpha,\beta))}{\partial\beta}=0\end{cases}\end{align}\]

Với hàm hợp lý: \(L(\alpha,\beta) = \prod \limits_{i}\pi^{y_i}(x_i).[1- \pi(x_i)]^{1-y_i} (\pi(x)=\frac{1}{1+e^{-\alpha-\beta x}})\)

Ví dụ 2: Đối với dữ liệu bệnh ngáy và bệnh tim ở Bảng 3.1, phần mềm cho thấy sự phù hợp của ML đối với mô hình hồi quy logistic

\[\begin{array}{lrllr} \hline \text{Dependent Variable}: \text{Y} \\ \text{Method: ML - Binary Probit (Quadratic hill climbing)} \\ \text{Included observations}: 2484 & & & & \\ \hline \text{Variable} & \text{Coefficient} & \text{Std. Error} & \text{z-Statistic} & \text{Prob.} \\ \hline\hline \text{C} & -3.866248 & 0.397337 & 0.166214 & 0.050011 \\ \text{X} & -23.26061 & 7.945039 & ~~~40.0000 & 0.0000 \\ \hline\hline \text{McFadden R-squared} & 0.070042 & \text{Mean dependent var} & & 0.044283 & \\ \text{Akaike info criterion} & 0.338861 & \text{Sum squared resid} & & 102.1917 & \\ \text{Schwarz criterion} & 0.343545 & \text{Log likelihood} & & -418.8658 & \\ \text{Hannan-Quinn criter.} & 0.340562 & \text{Deviance} & & 837.7316 & \\ \text{Restr. deviance} & 900.8272 & \text{Restr. log likelihood} & & -450.4136 & \\ \text{LR statistic} & 63.09557 & \text{Avg. log likelihood} & & -0.168626 & \\ \text{Prob(LR statistic)} & 0.000000 & & & \\ \hline \text{Obs with Dep}=0 & 2374 & \text{Total obs} & & 2484 & \\ \text{Obs with Dep}=1 & 110 \\ \hline \end{array}\]

Bảng 3.4. Hồi quy logistic cho dữ liệu bảng 3.1

\[logit[\hat\pi(x)] = −3,866248+0,397337. x \] - Diễn giải kết quả hồi quy: - Giá trị \(\hat\beta=0,397337\) phản ánh sự gia tăng nguy cơ mắc bệnh tim ở mức độ ngáy cao hơn.

• Tính xác suất mắc bệnh tim ứng với các mức ngáy:

  • \(x = 0\), \(logit[\hat\pi(0)] = −3,866248+0,397337.0=−3,866248\) nên \(\hat\pi(0)=\frac{e^{-3,866248}}{1 + e^{-3,866248}} = 0,020507417\), nghĩa là: xác suất mắc bệnh tim đối với người không ngáy là khoảng 2,051%.

  • \(x = 2\), \(logit[\hat\pi(2)] = −3,866248+0,397337.2=−3,071574\) nên \(\hat\pi(2)=\frac{e^{-3,071574}}{1 + e^{-3,071574}} = 0,044295\), nghĩa là xác suất mắc bệnh tim đối với người thỉnh thoảng ngáy là khoảng 4,43%.

  • \(x = 4\), \(logit[\hat\pi(4)] = −3,866248+0,397337.4=−2,2769\) nên \(\hat\pi(2)=\frac{e^{−2,2769}}{1 + e^{−2,2769}} = 0,093054\), nghĩa là xác suất mắc bệnh tim đối với người ngáy gần như mỗi đêm là khoảng 9,31%.

  • \(x = 5\), \(logit[\hat\pi(4)] = −3,866248+0,397337.5=−1,879563\) nên \(\hat\pi(2)=\frac{e^{−1,879563}}{1 + e^{−1,879563}} = 0,132439\),nghĩa là xác suất mắc bệnh tim đối với người ngáy mỗi đêm là khoảng 13,244%.

Khoảng tin cậy cho các hệ số trong mô hình được tính bằng phần mềm và kết quả như sau:

data_frame5 <- data.frame(
  Variable = c("C", "N"),
  Coefficient = c(-3.866248, 0.397337),
  Low = c(-4.139748, 0.315046),
  High = c(-3.592748, 0.479628),
  Low = c(-4.192181, 0.299270),
  High = c(-3.540315, 0.495404),
  Low = c(-4.294717, 0.268419),
  High = c(-3.437779, 0.526255)
)

data_frame5 %>%
  kbl(caption = "Coefficient Confidence Intervals") %>%
  kable_classic_2(full_width = F) %>%
  add_header_above(c(" " = 1, " " = 1, "90% CI" = 2, "95% CI" = 2, "99% CI" = 2))

Coefficient Confidence Intervals

		90% CI		95% CI		99% CI
Variable	Coefficient	Low	High	Low.1	High.1	Low.2	High.2
C	-3.8662	-4.140	-3.5927	-4.1922	-3.5403	-4.2947	-3.4378
N	0.3973	0.315	0.4796	0.2993	0.4954	0.2684	0.5263

Các liên kết nhị phân thay thế

Đối với các đường cong hồi quy logit được minh họa trong Hình 3.2, xác suất thành công tăng liên tục hoặc giảm liên tục khi tăng x. Từ mô hình logit: \[\log(\frac{\pi(x)}{1-\pi(x)})=\alpha+\beta x \]

ta có: \(\pi(x)=\frac{1}{1+e6{-\alpha+\beta x}}= F(x)\)

thì F(x) biến thiên từ 0 đến 1 và là hàm đơn điệu (tăng nếu $ > 0 $, giảm nếu $ < 0 $). Nói cách khác \(\pi(x)\) (hoặc \(1-\pi(x)\)) là hàm phân phối xác suất của một biến ngẫu nhiên. Điều này gợi mở cho ý tưởng xem \(\pi(x)\) (hoặc \(1-\pi(x)\)) là hàm phân phối xác suất F(x) của một biến ngẫu nhiên nào đó:

\[\begin{equation*} \pi(x)=F(x) \tag{3.2.3} \end{equation*}\]

Khi thay đổi hàm phân phối F(x), (3.2.3) xác định một lớp các mô hình.

Mô hình từ (3.2.3) xảy ra một cách tự nhiên khi phân phối dung sai (hay mức chịu đựng) áp dụng cho các đáp ứng của đối tượng. Chẳng hạn, đối với một loại côn trùng, mỗi con có thể chịu đựng một lượng thuốc diệt côn trùng nhất định và nó sẽ chết nếu liều lượng thuốc vượt quá mức chịu đựng. Như vậy lượng X về thuốc diệt côn trùng mà một con côn trùng có thể chịu đựng là một biến ngẫu nhiên mà: \(\pi(x)=P(X \leq x)= F(x)\) là xác suất để một con côn trùng loại này chết ở mức độ 𝑥 về lượng thuốc, tức là mô hình xác suất này có dạng (3.2.3). Chẳng hạn, nếu mức chịu đựng khác nhau giữa các con muỗi theo phân phối logistic, thì mô hình hồi quy logit sẽ được áp dụng.

Các mô hình Probit

Hồi quy Probit, đôi khi còn được gọi là mô hình probit, được sử dụng để ước lượng mô hình có biến phụ thuộc dạng nhị phân. Trong hồi quy probit, xem hàm xác suất \(\pi(x)\) là hàm phân phối xác suất của một biến ngẫu nhiên có phân phối chuẩn, tức là có dạng:

\[\begin{equation*} \pi(x)= \Phi(\alpha+\beta x) \tag{3.2.4} \end{equation*}\]

trong đó \(\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^x \mathrm{e}^{-\frac{-1}{2}t^2}\,\mathrm{d}t\) là hàm chuẩn tắc.

Ký hiệu: \(Probit(\pi(x))=\Phi^{-1}(\pi(x))\) và mô hình Probit được viết dươi dạng: \[Probit(\pi(x))=\alpha+\beta x\] Khi biết giá trị của xác suất \(\pi\), ta tìm được probit tương ứng và ngược lại. Chẳng hạn \(\pi(x)=0,05\), thì \(prpbit(0,05) = \Phi^{-1} (0,05) = −1,645\) ; $ = 0,5 $, thì \(probit(0,5)=\Phi^{-1}(0,5)\) ; $ probit((x))=1,96 $ thì \(\pi(x)=\Phi(1,96)=0,975,...\)

Mô hình Probit là một GLM với thành phần ngẫu nhiên nhị phân và liên kết probit.

Việc ước lượng mô hình Probit chuyển về ước lượng các tham số \(\alpha,\beta\). Để ước lượng cho hai tham số này, người ta sử dụng phương pháp ML, trong đó ước lượng hợp lý cực đại cho \(\alpha,\beta\) là nghiệm của phương trình hợp lý:

\[\begin{equation*} \begin{cases} \frac{\partial \log(L(\alpha,\beta))}{\partial \alpha} = 0 \\ \frac{\partial \log(L(\alpha,\beta))}{\partial \beta} = 0 \end{cases} \tag{3.2.5} \end{equation*}\]

với hàm lý là: \[L(\alpha,\beta)= \prod\limits_{i}\pi^{y_i}[1-\pi(x_i)]^{1-y_i} \space\space\space\space\space(\pi(x))=\Phi(\alpha+\beta x)\]

Việc giải (3.2.5) và do đó việc chạy mô hình probit phải dựa vào các phần mềm thống kê.

Ví dụ 3: Chúng ta minh họa bằng cách sử dụng dữ liệu bệnh ngáy và bệnh tim. ML phù hợp với mô hình probit, sử dụng điểm số {0,2,4,5} cho mức ngáy.

Từ bảng 3.5. về kết quả hồi quy ta có mô hình Probit ước lượng: \[Probit[\hat\pi(x)]= −2,060552+0,187770 x\]

• Khi \(x = 0 \ \Rightarrow \hat\pi(2)=\Phi((−2,061)=0,0197\)
• Khi \(x = 2 \ \Rightarrow \hat\pi(2)=\Phi((−1,6865)=0,0455\)
• Khi \(x = 4 \ \Rightarrow \hat\pi(4)=\Phi(−1,309)= 0,0951\)
• Khi \(x = 5 \ \Rightarrow \hat\pi(5)=\Phi(−1,122)=0,1314\)

\[\begin{array}{lrllr} \hline \text{Dependent Variable}: \text{Y} \\ \text{Method: ML - Binary Probit (Quadratic hill climbing)} \\ \text{Included observations}: 2484 \\ \text{Convergence achieved after 5 iterations}\\ \text{Covariance matrix computed using second derivatives} \\ & & & & \\ \hline \text{Variable} & \text{Coefficient} & \text{Std. Error} & \text{z-Statistic} & \text{Prob.} \\ \hline\hline \text{C} & -3.866248 & 0.397337 & 0.166214 & 0.050011 \\ \text{X} & -23.26061 & 7.945039 & ~~~40.0000 & 0.0000 \\ \hline\hline \text{McFadden R-squared} & 0.071082 & \text{Mean dependent var} & & 0.044283 & \\ \text{Akaike info criterion} & .338484 & \text{Sum squared resid} & & 102.1504 & \\ \text{Schwarz criterion} & 0.343168 & \text{Log likelihood} & & -418.3971 & \\ \text{Hannan-Quinn criter.} & 0.340185 & \text{Deviance} & & 836.7943 & \\ \text{Restr. deviance} & 900.8272 & \text{Restr. log likelihood} & & -450.4136 & \\ \text{LR statistic} & 64.03292 & \text{Avg. log likelihood} & & -0.168437 & \\ \text{Prob(LR statistic)} & 0.000000 & & & \\ \hline \text{Obs with Dep}=0 & 2374 & \text{Total obs} & & 2484 & \\ \text{Obs with Dep}=1 & 110 \\ \hline \end{array}\]

Bảng 3.5. Hồi quy Probit cho dữ liệu bảng 3.1

Ước lượng khoảng tin cậy cho các hệ số trong mô hình:

data_frame6 <- data.frame(
  Variable = c("C", "X"),
  Coefficient = c(-2.060552 , 0.187770),
  Low = c(-2.176474, 0.148888),
  High = c(-1.944630, 0.226653),
  Low = c(-2.1998697, 0.141434),
  High = c(-1.922406, 0.234107),
  Low = c(-2.242156, 0.126857),
  High = c(-1.878947, 0.248684)
)

data_frame6 %>%
  kbl(caption = "Coefficient Confidence Intervals") %>%
  kable_classic_2(full_width = F) %>%
  add_header_above(c(" " = 1, " " = 1, "90% CI" = 2, "95% CI" = 2, "99% CI" = 2))

Coefficient Confidence Intervals

		90% CI		95% CI		99% CI
Variable	Coefficient	Low	High	Low.1	High.1	Low.2	High.2
C	-2.0606	-2.1765	-1.9446	-2.1999	-1.9224	-2.2422	-1.8789
X	0.1878	0.1489	0.2267	0.1414	0.2341	0.1269	0.2487

Xác suất bị bệnh tim ứng với các mức ngáy được ước lượng qua các mô hình: Xác suất tuyến tính, logit và probit được tổng kết trong bảng 3.2 sau:

Ngáy	Bệnh tim		Tỷ lệ Có	Ước lượng tuyến tính	Ước lượng logit	Ước lượng probit
	Có	Không
Không ngáy	24	1355	0.0174	0.0172	0.0204	0.0197
Thỉnh thoảng	35	603	0.0549	0.0568	0.0444	0.0455
Gần như mỗi đêm	21	192	0.0990	0.0964	0.0936	0.0951
Mỗi đêm	30	224	0.118	0.1161	0.1335	0.1314

Bảng 3.6 Ước lượng xác suất bị bệnh tim theo các mức ngáy

Bảng 3.6 cho thấy rằng các ươc lượng probit và logistic khá gần nhau và mô hình tuyến tính phù hợp với những dữ liệu này. (Phần sau thảo luận về các mô hình phù hợp cho việc phân tích các GLM với biến đáp ứng là biến nhị phân).

Các mô hình tuyến tính tổng quát cho dữ liệu đếm

Nhiều biến đáp ứng rời rạc đóng vai trò đếm số lần xuất hiện một biến cố ngẫu nhiên A nào đó trong một khoảng thời gian hay không gian hay một bối cảnh nhất định, hay nói cách khác là các biến ngẫu nhiên có phân phối Poisson. Những dữ liệu thu được khi quan sát các biến dạng này thường được gọi là dữ liệu đếm. Ví dụ, số vụ tai nạn giao thông xảy ra trên một cung đường trong một ngày, số vụ trộm xe ô tô trong năm 2017, hay số khách hàng đến một hệ thống dịch vụ trong một tuần làm việc,…. Phần này giới thiệu GLM cho dữ liệu đếm. Các GLM giả sử một phân phối Poisson cho thành phần ngẫu nhiên. Cũng giống như các số đếm, biến ngẫu nhiên Poisson có thể lấy bất kỳ giá trị số nguyên không âm.