Chương 13 Hồi quy biến phụ thuộc là nhị phân

Trong chương này, ta sẽ xem xét một loại hồi quy đặc biệt nhằm giải thích biến phụ thuộc có giá trị rời rạc hoặc bị giới hạn. Cụ thể ta sẽ xem xét các mô hình sau:

  • Mô hình xác suất tuyến tính
  • Mô hình logit
  • Mô hình probit
  • phương pháp MLE (maximum likelihood estimation) cho mô hình hồi quy phi tuyến

Trong chương này, ta sẽ đặt vấn đề: liệu có sự phân biệt chủng tộc trong thị trường thế chấp Mỹ hay không.

Dữ liệu được sử dụng từ gói AER (Kleiber & Zeileis, 2018)12

13.1 Mô hình xác suất tuyến tính

Mô hình xác suất tuyến tính có dạng:

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \dots + \beta_k X_{ki} + u_{i}\]

với \(Y_i\) là biến nhị phân được gọi là mô hình xác suất tuyến tính. Cụm từ xác suất được lấy từ công thức tính kỳ vọng của \(Y_i\):

\[E(Y|X_i) = P(Y=1 | X_i)\]

trong đó:

\[P(Y=1 | X_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \dots + \beta_k X_{ki}\]

Theo đó, \(\beta_i\) có thể được giải thích là sự thay đổi xác suất mà \(Y_i=1\), giữ cố định các biến \(X\) khác. Phương pháp ước lượng giống như hồi quy đa biến. Cần chú ý là ta luôn phải sử dụng sai số chuẩn vững (robust standard errors) bởi vì \(u_i\) luôn luôn có phương sai thay đổi.

13.1.1 Dữ liệu thế chấp

Ta khảo sát dữ liệu hồ sơ đi vay thế chấp ở Boston trong năm 1990.

data(HMDA)
# inspect the data
head(HMDA)
##   deny pirat hirat     lvrat chist mhist phist unemp selfemp insurance
## 1   no 0.221 0.221 0.8000000     5     2    no   3.9      no        no
## 2   no 0.265 0.265 0.9218750     2     2    no   3.2      no        no
## 3   no 0.372 0.248 0.9203980     1     2    no   3.2      no        no
## 4   no 0.320 0.250 0.8604651     1     2    no   4.3      no        no
## 5   no 0.360 0.350 0.6000000     1     1    no   3.2      no        no
## 6   no 0.240 0.170 0.5105263     1     1    no   3.9      no        no
##   condomin afam single hschool
## 1       no   no     no     yes
## 2       no   no    yes     yes
## 3       no   no     no     yes
## 4       no   no     no     yes
## 5       no   no     no     yes
## 6       no   no     no     yes
summary(HMDA)
##   deny          pirat            hirat            lvrat        chist   
##  no :2095   Min.   :0.0000   Min.   :0.0000   Min.   :0.0200   1:1353  
##  yes: 285   1st Qu.:0.2800   1st Qu.:0.2140   1st Qu.:0.6527   2: 441  
##             Median :0.3300   Median :0.2600   Median :0.7795   3: 126  
##             Mean   :0.3308   Mean   :0.2553   Mean   :0.7378   4:  77  
##             3rd Qu.:0.3700   3rd Qu.:0.2988   3rd Qu.:0.8685   5: 182  
##             Max.   :3.0000   Max.   :3.0000   Max.   :1.9500   6: 201  
##  mhist    phist          unemp        selfemp    insurance  condomin  
##  1: 747   no :2205   Min.   : 1.800   no :2103   no :2332   no :1694  
##  2:1571   yes: 175   1st Qu.: 3.100   yes: 277   yes:  48   yes: 686  
##  3:  41              Median : 3.200                                   
##  4:  21              Mean   : 3.774                                   
##                      3rd Qu.: 3.900                                   
##                      Max.   :10.600                                   
##   afam      single     hschool   
##  no :2041   no :1444   no :  39  
##  yes: 339   yes: 936   yes:2341  
##                                  
##                                  
##                                  
## 

Biến ta quan tâm là deny, nếu hồ sơ bị từ chối cho vay thì deny = yes. Một biến giải thích có năng lực giải thích tốt đó là pirat, được định nghĩa tổng khoản trả hằng tháng so với thu nhập khách hàng. Mô hình được đề xuất như sau:

\[deny = \beta_0 + \beta_1 \times P/I ratio + u\]

Ta sử dụng câu lệnh lm() như bình thường, nhưng phải đảm bảo biến deny được chuyển đổi thành dạng numeric bằng cách sử dụng as.numeric() bởi vì lm() không chấp nhận dạng biến factor. Chú ý thêm rằng, hàm as.numeric() sẽ biến giá trị deny=no thành giá trị 1deny=yes thành giá trị 2, do đó ta dùng as.numeric()-1 để trả biến về giá trị 01.

# convert 'deny' to numeric
HMDA$deny <- as.numeric(HMDA$deny) - 1

# estimate a simple linear probabilty model
denymod1 <- lm(deny ~ pirat, data = HMDA)
denymod1
## 
## Call:
## lm(formula = deny ~ pirat, data = HMDA)
## 
## Coefficients:
## (Intercept)        pirat  
##    -0.07991      0.60353

Tiếp theo, ta vẽ đồ thị hồi quy.

# plot the data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Scatterplot Mortgage Application Denial and the Payment-to-Income Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.8)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add the estimated regression line
abline(denymod1, 
       lwd = 1.8, 
       col = "steelblue")

Theo mô hình, payment-to-income ratio có giá trị là \(1\) sẽ tương ứng xác suất từ chối phê duyệt kỳ vọng khoảng hơn \(50\%\). Mô hình cho thấy có mối liên hệ dương giữa payment-to-income ratio và xác suất hồ sơ bị từ chối. Tuy nhiên ta cần sử dụng ước lượng vững và do đó đoạn code trong R thực hiện như sau.

# print robust coefficient summary
coeftest(denymod1, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.079910   0.031967 -2.4998   0.01249 *  
## pirat        0.603535   0.098483  6.1283 1.036e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mô hình được giải thích như sau: khi tăng P/I ratio lên \(1\%\) thì xác suất từ chối tăng lên \(0.604 \times 0.01 = 0.6 \%\).

Tiếp theo ta xem xét biến black bằng \(1\) nếu là người Mỹ Phi. Nếu khảo sát thống kê biến này có ý nghĩa, thì đây là bằng chứng cho thấy có tồn tại sự phân biệt chủng tộc trong quan hệ tín dụng.

# rename the variable 'afam' for consistency
colnames(HMDA)[colnames(HMDA) == "afam"] <- "black"

# estimate the model
denymod2 <- lm(deny ~ pirat + black, data = HMDA)
coeftest(denymod2, vcov. = vcovHC)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.090514   0.033430 -2.7076  0.006826 ** 
## pirat        0.559195   0.103671  5.3939 7.575e-08 ***
## blackyes     0.177428   0.025055  7.0815 1.871e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rõ ràng, có sự phân biệt chủng tộc ở đây. Tuy nhiên, kết quả có thể bị thiên lệch do thiếu biến giải thích. Do đó kết luận có phân biệt chủng tộc là còn quá sớm!

13.2 Hồi quy Probit và Logit

Hồi quy xác suất tuyến tính có một nhược điểm quan trọng: nó giả định hàm xác suất có điều kiện là tuyến tính. Điều này không đảm bảo xác suất sẽ nằm trong vùng giá trị từ \(0\) đến \(1\). Như đồ thị hồi quy, nếu \(P/I ratio \ge 1.75\) thì xác suất từ chối hồ sơ sẽ lớn hơn \(1\). Chính vì vậy, ta cần dùng định dạng phi tuyến cho mô hình, và hai kỹ thuật được sử dụng nhiều nhất là logitprobit.

13.2.1 Hồi quy Probit

Trong hồi quy Probit, hàm phân phối Normal chuẩn tích luỹ \(\Phi()\) được sử dụng để mô hình hàm hồi quy.

\[E(Y|X) = P(Y=1|X) = \Phi(\beta_0 + \beta_1 X)\]

Khi đó \(\beta_0 + \beta_1 X\) đóng vai trò như một quantile \(z\).

Mô hình Probit khá là khó giải thích ý nghĩa hệ số hồi quy của \(X\). Ta cần giải thích thông quan nhiều bước tính như sau:

  • Tính xác suất dự đoán \(Y=1\) đối với một giá trị \(X\).
  • Tính xác suất dự đoán \(Y=1\) đối với một giá trị \(X + \Delta X\).
  • Tính chênh lệch kết quả hai xác suất vừa tính.

Trong R ta cần thực hiện như sau.

# estimate the simple probit model
denyprobit <- glm(deny ~ pirat, 
                  family = binomial(link = "probit"), 
                  data = HMDA)

coeftest(denyprobit, vcov. = vcovHC, type = "HC1")
## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -2.19415    0.18901 -11.6087 < 2.2e-16 ***
## pirat        2.96787    0.53698   5.5269 3.259e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta vẽ đồ thị hồi quy như sau.

# plot data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Probit Model of the Probability of Denial, Given P/I Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.85)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add estimated regression line
x <- seq(0, 3, 0.01)
y <- predict(denyprobit, list(pirat = x), type = "response")

lines(x, y, lwd = 1.5, col = "steelblue")

Ta tiếp tục bằng mô hình augmented Probit để ước lượng ảnh hưởng của chủng tộc lên xác suất từ chối hồ sơ.

denyprobit2 <- glm(deny ~ pirat + black, 
                   family = binomial(link = "probit"), 
                   data = HMDA)

coeftest(denyprobit2, vcov. = vcovHC, type = "HC1")
## 
## z test of coefficients:
## 
##              Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -2.258787   0.176608 -12.7898 < 2.2e-16 ***
## pirat        2.741779   0.497673   5.5092 3.605e-08 ***
## blackyes     0.708155   0.083091   8.5227 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.2.2 Hồi quy Logit

Mô hình hồi quy Logit được mô tả như sau.

\[P(Y=1|X_i) = F(\beta_0 + \sum \beta_i X_i) = \frac{1}{1+e^{-(\beta_0 + \sum \beta_i X_i)}}\]

Ý tưởng vẫn giống như mô hình Probit nhưng hàm số thay vì \(\Phi()\) thì ta dùng một hàm \(F\):

\[F(x) = \frac{1}{1+e^{-x}}\]

denylogit <- glm(deny ~ pirat, 
                 family = binomial(link = "logit"), 
                 data = HMDA)

coeftest(denylogit, vcov. = vcovHC, type = "HC1")
## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -4.02843    0.35898 -11.2218 < 2.2e-16 ***
## pirat        5.88450    1.00015   5.8836 4.014e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Đồ thị như sau.

# plot data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Probit and Logit Models Model of the Probability of Denial, Given P/I Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.9)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add estimated regression line of Probit and Logit models
x <- seq(0, 3, 0.01)
y_probit <- predict(denyprobit, list(pirat = x), type = "response")
y_logit <- predict(denylogit, list(pirat = x), type = "response")

lines(x, y_probit, lwd = 1.5, col = "steelblue")
lines(x, y_logit, lwd = 1.5, col = "black", lty = 2)

# add a legend
legend("topleft",
       horiz = TRUE,
       legend = c("Probit", "Logit"),
       col = c("steelblue", "black"), 
       lty = c(1, 2))

Bây giờ ta thêm biến black vào mô hình Logit.

# estimate a Logit regression with multiple regressors
denylogit2 <- glm(deny ~ pirat + black, 
                  family = binomial(link = "logit"), 
                  data = HMDA)

coeftest(denylogit2, vcov. = vcovHC, type = "HC1")
## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -4.12556    0.34597 -11.9245 < 2.2e-16 ***
## pirat        5.37036    0.96376   5.5723 2.514e-08 ***
## blackyes     1.27278    0.14616   8.7081 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Trong thực tiễn không có sự thiên vị nào giữa hai mô hình Probit và Logit, vì cơ bản ta cũng không thể nào biết được hàm hồi quy tổng thể của \(E(Y|X)\), ở đây ta đang cố gắng đưa vào yếu tố phi tuyến để khắc phục nhược điểm của hồi quy xác suất tuyến tính.

13.3 Ước lượng và suy luận trong hồi quy Logit và Probit

Vì cả hai mô hình đều phi tuyến đối với tham số do đó ta không thể sử dụng OLS để ước lượng. Thay vào đó một số người sử dụng ước lượng likelihood cực đại (MLE), một số khác sử dụng bình phương nhỏ nhất phi tuyến ( NLS - nonlinear least squares).

13.3.1 NLS

Tương tự như OLS, NLS ước lượng các tham số thông qua cực tiểu hoá tổng sai số:

\[\sum [Y_i - \Phi(b_0 + b_1 X_1 + \dots + b_k X_k)]^2\]

Ước lượng NLS là consistent vì cho ra ước lượng tuân theo phân phối Normal khi mẫu đủ lớn. Trong R hàm nls() từ package stats cung cấp thuật toán để giải quyết bài toán NLS. Tuy nhiên NLS lại không efficient nghĩa là có một kỹ thuật khác cho ra ước lượng có phương sai nhỏ hơn.

13.3.2 MLE

Trong MLE, ta tìm kiếm các ước lượng chưa biết bằng cách chọn những giá trị mà likelihood của việc lấy mẫu quan sát được là lớn nhất. Xác suất được đo lường bằng hàm likelihood. Nói cách khác, ước lượng MLE của những tham số chưa biết là các giá trị mà cho ra một mô hình có khả năng cao nhất sẽ cho ra dữ liệu đang quan sát. Do đó MLE sẽ efficient hơn NLS.

Bởi vì cùng có tính consistency nên suy luận dựa trên NLS hay MLE đều giống như OLS: ta có thể dùng thống kê \(t\) để tính toán khoảng tin cậy.

Nhiều phần mềm thống kê sử dụng các thuật toán MLE khác nhau, trong đó hàm glm() của R sử dụng thuật toán gọi là iteratively reweighted least squares.

13.3.3 Đo lường hiệu quả mô hình

Ta cần nhớ rằng \(R^2\)\(R^2\) hiệu chỉnh thông thường không áp dụng được đối với những mô hình hồi quy phi tuyến. Bởi cả hai thước đo đều giả định ngầm rằng quan hệ giữa các biến là tuyến tính.

Có rất nhiều thước đo đo lường hiệu quả mô hình thay vì dùng \(R^2\)\(\bar{R}^2\) nhưng không hề có thước đo nào vượt trội hơn cả, kể cả về ý nghĩa và cách giải thích. Phương pháp sơ khai nhất là dự đoán \(P(Y|X) > 0.5\) thì xem như \(Y_i = 1\)\(P(Y|X) <0.5\) thì xem như \(Y_i = 0\), còn lại thì không xác định. Tuy nhiên phương pháp này lại không đánh giá hiệu quả mô hình ở mức độ so sánh xác suất, nghĩa là chất lượng đầu ra một xác suất bằng \(0.51\) cũng hệt như chất lượng của một đầu ra xác suất \(0.99\).

Một phương pháp khác cũng khá nổi tiếng đó là sử dụng pseudo-\(R^2\), thay vì đo lường hiệu quả mô hình, ta so sánh giá trị logarithm của likelihood mô hình có đủ biến độc lập với mô hình không có biến độc lập nào. Công thức như sau.

\[\text{pseudo-}R^2= 1-\frac{\ln (f_{full}^{max})}{\ln(f_{null}^{max})}\]

Nguyên nhân là vì khi thêm một biến giải thích vào mô hình thì likelihood cực đại sẽ tăng lên, tương tự như việc giảm \(SSR\) trong trường hợp hồi quy tuyến tính. Nếu một mô hình \(full\) có likelihood cực đại như một mô hình \(null\) thì \(\text{pseudo-}R^2\approx 0\).

Nếu ta gọi:

\[\begin{aligned} \text{deviance} &= -2 \times [\ln (f_{saturated}^{max}) - \ln (f_{full}^{max})] \\ \text{null deviance} &= -2\times [\ln (f_{saturated}^{max}) - \ln (f_{null}^{max})] \end{aligned}\]

thì:

\[\text{pseudo-}R^2= 1 - \frac{\text{deviance}}{\text{null deviance}}\]

Trong R ta tính \(\text{pseudo-}R^2\) như sau.

# compute pseudo-R2 for the probit model of mortgage denial
pseudoR2 <- 1 - (denyprobit2$deviance) / (denyprobit2$null.deviance)
pseudoR2
## [1] 0.08594259

13.4 Câu hỏi lập trình

13.4.1 Câu 1: Dữ liệu sống sót sau sự kiện Titanic

Trong bộ phim Titanic, có một chính sách là ưu tiên phụ nữ và trẻ em trước khi tàu chuẩn bị chìm, tuy nhiên thực trạng cho thấy có rất nhiều phụ nữ và trẻ em ở tầng lớp thứ ba không được cứu sống.

  1. Phân bổ dữ liệu Titanic vào biến Titanic_1 và nhìn tổng quan về dữ liệu.
  2. Hình hoạ hoá tỷ lệ sống sót có điều kiện đối với tầng lớp Class, giới tính Sex, và tuổi tác Age bằng cách sử dụng mosaicplot().
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCB0aGUgZGF0YXNldCBgVGl0YW5pY2BcblxuIyBnZXQgYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YVxuXG4jIHBsb3QgYSBtb3NhaWMiLCJzb2x1dGlvbiI6IiMgYXR0YWNoIHRoZSBkYXRhc2V0IGBUaXRhbmljYFxuVGl0YW5pY18xIDwtIFRpdGFuaWNcblxuIyBnZXQgYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YVxuc3VtbWFyeShUaXRhbmljXzEpXG4jIG9yXG5zdHIoVGl0YW5pY18xKVxuIyBvclxuaGVhZChUaXRhbmljXzEpXG5cbiMgcHJpbnQgdGhlIGRhdGEgdGhlIGNvbnNvbGVcblRpdGFuaWNfMVxuXG4jIHBsb3QgYSBtb3NhaWNcbm1vc2FpY3Bsb3QoVGl0YW5pY18xLCBtYWluID0gXCJTdXJ2aXZhbCBvbiB0aGUgVGl0YW5pY1wiKSJ9

13.4.2 Câu 2: Dữ liệu sống sót sau sự kiện Titanic - Ctd

Dữ liệu Titanic ở trên chưa đủ đẹp để có thể chạy hồi quy. Ta sẽ sử dụng dữ liệu từ đường link sau: https://stanford.io/2O9RUCF .

File titanic.csv có chứa những biến như sau:

  • Survived - chỉ báo sống sót.
  • Pclass - hạng hành khách.
  • Name - tên hành khách.
  • Sex - giới tính hành khách.
  • Age - tuổi hành khách.
  • Siblings - số lượng anh em ở nước ngoài.
  • Parents.Children.Aboard - số lượng ba mẹ và con cái ở nước ngoài.
  • fare - tiền phí tàu.
  1. Nhập dữ liệu từ file csv bằng hàm read.csv2(). Lưu vào biến Titanic_2.
  2. Phân bổ tên cột vào Titanic_2: Survived, Class, Name, Sex, Age, Siblings, Parents and Fare.
  3. Khảo sát dữ liệu. Loại bỏ cột Name.
  4. Gọi package corrplotdplyr. Khảo sát đa cộng tuyến bằng hàm corrplot().
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCB0aGUgcGFja2FnZXMgYGRwbHlyYCBhbmQgYHJlYWRyYFxuICAgXG4gIFxuIyB1c2UgYHJlYWRfY3N2KClgIHRvIGltcG9ydCB0aGUgZGF0YXNldCBhbmQgYXNzaWduIHRoZSBkYXRhIHRvIGBUaXRhbmljXzJgXG5cblxuIyBnZXQgYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YSBhbmQgZHJvcCBgTmFtZWBcblxuICAgICAgICAgIFxuIyBjaGFuZ2UgdGhlIGNvbHVtbiBuYW1lc1xuXG5cbiMgYXR0YWNoIHRoZSBwYWNrYWdlIGBjb3JycGxvdGBcblxuXG4jIGNoZWNrIGNvcnJlbGF0aW9ucyB1c2luZyBgY29ycnBsb3QoKWAiLCJzb2x1dGlvbiI6IiMgYXR0YWNoIHRoZSBwYWNrYWdlcyBgZHBseXJgIGFuZCBgcmVhZHJgXG5saWJyYXJ5KHJlYWRyKVxubGlicmFyeShkcGx5cilcbiAgXG4jIHVzZSBgcmVhZF9jc3YoKWAgdG8gaW1wb3J0IHRoZSBkYXRhc2V0IGFuZCBhc3NpZ24gdGhlIGRhdGEgdG8gYFRpdGFuaWNfMmBcblRpdGFuaWNfMiA8LSByZWFkX2NzdihcImh0dHBzOi8vc3RhbmZvcmQuaW8vMk85UlVDRlwiKVxuXG4jIGdldCBhbiBvdmVydmlldyBvdmVyIHRoZSBkYXRhIGFuZCBkcm9wIGBOYW1lYFxuc3VtbWFyeShUaXRhbmljXzIpXG4jIG9yXG5zdHIoVGl0YW5pY18yKVxuIyBvclxuaGVhZChUaXRhbmljXzIpXG5cblRpdGFuaWNfMiA8LSBUaXRhbmljXzJbLCAtM11cblxuIyBjaGFuZ2UgdGhlIGNvbHVtbiBuYW1lc1xuY29sbmFtZXMoVGl0YW5pY18yKSA8LSBjKFwiU3Vydml2ZWRcIiwgXCJDbGFzc1wiLCBcIlNleFwiLCBcIkFnZVwiLCBcIlNpYmxpbmdzXCIsIFwiUGFyZW50c1wiLCBcIkZhcmVcIilcblxuIyBhdHRhY2ggdGhlIHBhY2thZ2UgYGNvcnJwbG90YFxubGlicmFyeShjb3JycGxvdClcbiAgXG4jIGNoZWNrIGNvcnJlbGF0aW9ucyB1c2luZyBgY29ycnBsb3QoKWBcbmNvcnJwbG90KGNvcihzZWxlY3RfaWYoVGl0YW5pY18yLCBpcy5udW1lcmljKSkpXG4jICh0aGUgaGlnaGVzdCBjb3JyZWxhdGlvbiBpcyBiZXR3ZWVuIGZhcmUgYW5kIHBhc3NlbmdlciBjbGFzcykifQ==

13.4.3 Câu 3: Dữ liệu sống sót sau sự kiện Titanic - Tỷ lệ sống sót

Bảng thống kê rất phù hợp để cho ta thấy phân phối có điều kiện của sống sót với các biến giải thích. Trong R bảng thống kê được thực hiện qua hàm table().

  1. Lập bảng thống kê giữa SurvivedClass. Lưu vào t_abs.
  2. Bảng t_abs báo cáo tần suất tuyệt đối. Chuyển đổi bảng thành dữ liệu tương đối, lưu vào biến t_rel.
  3. Hình hoạ hoá t_rel bằng barplot().
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGdlbmVyYXRlIGB0X2Fic2AsIGEgY29udGluZ2VuY3kgdGFibGUgb2YgYFN1cnZpdmVkYCBhbmQgYENsYXNzYFxuXG4jIGdlbmVyYXRlIGB0X3JlbGAsIHRoZSB0YWJsZSBvZiB0aGUgcmVsYXRpdmUgZnJlcXVlbmNpZXNcblxuIyBjcmVhdGUgdGhlIGJhcnBsb3QiLCJzb2x1dGlvbiI6IiMgZ2VuZXJhdGUgYHRfYWJzYCwgYSBjb250aW5nZW5jeSB0YWJsZSBvZiBgU3Vydml2ZWRgIGFuZCBgQ2xhc3NgXG50X2FicyA8LXRhYmxlKFRpdGFuaWNfMiRTdXJ2aXZlZCwgVGl0YW5pY18yJENsYXNzKVxuXG4jIGdlbmVyYXRlIGB0X3JlbGAsIHRoZSB0YWJsZSBvZiB0aGUgcmVsYXRpdmUgZnJlcXVlbmNpZXNcbnRfcmVsIDwtIHRfYWJzL25yb3coVGl0YW5pY18yKVxuXG4jIGNyZWF0ZSB0aGUgYmFycGxvdFxuYmFycGxvdCh0YWJsZShUaXRhbmljXzIkU3Vydml2ZWQsIFRpdGFuaWNfMiRDbGFzcykvbnJvdyhUaXRhbmljXzIpLCBcbmNvbCA9IGMoXCJkYXJrcmVkXCIsXCJkYXJrZ3JlZW5cIiksXG55bGltID0gYygwLDAuNiksIFxubWFpbiA9IFwiUmVsYXRpdmUgRnJlcXVlbmNpZXMgb2YgU3Vydml2YWxcIixcbnhsYWIgPSBcIkNsYXNzXCIpIn0=

13.4.4 Câu 4: Dữ liệu sống sót sau sự kiện Titanic - Phân phối có điều kiện của Age

Trong bài này ta sẽ hình hoạ phân phối của Age dựa vào điều kiện của Survived.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiJUaXRhbmljXzIgPC0gcmVhZF9jc3YoXCJodHRwczovL3N0YW5mb3JkLmlvLzJPOVJVQ0ZcIilcbiMgZXN0aW1hdGUgYm90aCBkZW5zaXRpZXMgYW5kIHNhdmUgdGhlIG91dGNvbWVzIHRvIGBkZW5zX2FnZV9zdXJ2YCBhbmQgYGRlbnNfYWdlX2RpZWRgXG5cbiMgcGxvdCBib3RoIGRlbnNpdGVzIHRvZ2V0aGVyIiwic29sdXRpb24iOiJUaXRhbmljXzIgPC0gcmVhZF9jc3YoXCJodHRwczovL3N0YW5mb3JkLmlvLzJPOVJVQ0ZcIilcbiMgZXN0aW1hdGUgYm90aCBkZW5zaXRpZXMgYW5kIHNhdmUgdGhlIG91dGNvbWVzIHRvIGBkZW5zX2FnZV9zdXJ2YCBhbmQgYGRlbnNfYWdlX2RpZWRgXG5kZW5zX2FnZV9zdXJ2IDwtIGRlbnNpdHkoVGl0YW5pY18yJEFnZVtUaXRhbmljXzIkU3Vydml2ZWQ9PTFdKVxuZGVuc19hZ2VfZGllZCA8LSBkZW5zaXR5KFRpdGFuaWNfMiRBZ2VbVGl0YW5pY18yJFN1cnZpdmVkPT0wXSlcblxuIyBwbG90IGJvdGggZGVuc2l0ZXMgdG9nZXRoZXJcbnBsb3QoZGVuc19hZ2Vfc3VydiwgXG55bGltID0gYygwLCAwLjAzNSksIFxuY29sID0gXCJkYXJrZ3JlZW5cIiwgXG5sd2QgPSAxLjUsIFxubWFpbiA9IFwiRGVuc2l0eSBFc3RpbWF0ZXMgb2YgXFwnQWdlXFwnXCIpXG5cbmxpbmVzKGRlbnNfYWdlX2RpZWQsIGNvbCA9IFwiZGFya3JlZFwiLCBsd2QgPSAxLjUpXG5cbmxlZ2VuZChsZWdlbmQgPSBjKFwiU3Vydml2ZWRcIiwgXCJEaWVkXCIpLCBcbnggPVwidG9wcmlnaHRcIiwgXG5jb2wgPSBjKFwiZGFya2dyZWVuXCIsIFwiZGFya3JlZFwiKSwgXG5sd2QgPSAxLjUpIn0=

13.4.5 Câu 5: Dữ liệu sống sót sau sự kiện Titanic - Mô hình xác suất tuyến tính 1

Ta cần khảo sát địa vị xã hội ảnh hưởng thế nào đến tỷ lệ sống sót.

  1. Gọi package AER.
  2. Đổi Class thành biến nhân tố.
  3. Ước lượng mô hình xác suất tuyến tính và lưu kết quả surv_mod.
  4. Ước lượng vững hệ số mô hình.
  5. Sử dụng mô hình để dự đoán xác suất sống sót cho ba cấp bậc xã hội.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCB0aGUgYEFFUmAgcGFja2FnZVxuXG5cbiMgZW5jb2RlIGBDbGFzc2AgYXMgYSBmYWN0b3JcblxuICBcbiMgZml0IHRoZSBsaW5lYXIgcHJvYmFiaWxpdHkgbW9kZWwsIGFzc2lnbiBpdCB0byBgc3Vydl9tb2RgXG5cblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgbW9kZWwgY29lZmZpY2llbnRzXG5cblxuIyBwcmVkaWN0IHRoZSBwcm9iYWJpbGl0eSBvZiBzdXJ2aXZhbCBmb3IgYWxsIHBhc3NlbmdlciBjbGFzc2VzIiwic29sdXRpb24iOiIjIGF0dGFjaCB0aGUgYEFFUmAgcGFja2FnZVxubGlicmFyeShBRVIpXG5cbiMgZW5jb2RlIGBDbGFzc2AgYXMgYSBmYWN0b3JcblRpdGFuaWNfMiRDbGFzcyA8LSBhcy5mYWN0b3IoVGl0YW5pY18yJENsYXNzKVxuXG4jIGZpdCB0aGUgbGluZWFyIHByb2JhYmlsaXR5IG1vZGVsLCBhc3NpZ24gaXQgdG8gYHN1cnZfbW9kYFxuc3Vydl9tb2QgPC0gbG0oU3Vydml2ZWQgfiBDbGFzcywgZGF0YSA9IFRpdGFuaWNfMilcblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgbW9kZWwgY29lZmZpY2llbnRzXG5jb2VmdGVzdChzdXJ2X21vZCwgdmNvdkhDKVxuXG4jIHByZWRpY3QgdGhlIHByb2JhYmlsaXR5IG9mIHN1cnZpdmFsIGZvciBhbGwgcGFzc2VuZ2VyIGNsYXNzZXNcbnByZWRpY3Qoc3Vydl9tb2QsIG5ld2RhdGEgPSBkYXRhLmZyYW1lKFwiQ2xhc3NcIiA9IGFzLmZhY3RvcigxOjMpKSkifQ==

13.4.6 Câu 6: Dữ liệu sống sót sau sự kiện Titanic - Mô hình xác suất tuyến tính 2

  1. Sử dụng surv_mod để có được ước lượng cụ thể cho xác suất sống sót.
  2. Ước lượng augmented LMP và lưu vào LMP_mod.
  3. Lấy hệ số vững.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbXB1dGUgdGhlIGNsYXNzIHNwZWNpZmljIGVzdGltYXRlcyAgIFxuICBcbiAgXG4jIGZpdCB0aGUgYXVnbWVudGVkIExNUCBhbmQgYXNzaWduIGl0IHRvIGBMUE1fbW9kYFxuXG5cbiMgb2J0YWluIGEgcm9idXN0IHN1bW1hcnkgb2YgdGhlIG1vZGVsIGNvZWZmaWNpZW50cyIsInNvbHV0aW9uIjoiIyBjb21wdXRlIHRoZSBjbGFzcyBzcGVjaWZpYyBlc3RpbWF0ZXMgICBcbnN1cnZfcHJvYl9jMSA8LSBzdXJ2X21vZCRjb2VmZmljaWVudHNbMV0gIFxuc3Vydl9wcm9iX2MyIDwtIHN1cnZfcHJvYl9jMSArIHN1cnZfbW9kJGNvZWZmaWNpZW50c1syXSAgXG5zdXJ2X3Byb2JfYzMgPC0gc3Vydl9wcm9iX2MxICsgc3Vydl9tb2QkY29lZmZpY2llbnRzWzNdICBcbiAgXG4jIGZpdCB0aGUgbGluZWFyIHByb2JhYmlsaXR5IG1vZGVsLCBhc3NpZ24gaXQgdG8gYHN1cnZfbW9kYFxuTFBNX21vZCA8LSBsbShTdXJ2aXZlZCB+IC4sIGRhdGEgPSBUaXRhbmljXzIpXG5cbiMgb2J0YWluIGEgcm9idXN0IHN1bW1hcnkgb2YgdGhlIG1vZGVsIGNvZWZmaWNpZW50c1xuY29lZnRlc3QoTFBNX21vZCwgdmNvdkhDKSJ9

13.4.7 Câu 7: Dữ liệu sống sót sau sự kiện Titanic - Mô hình Logit

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGZpdCB0aGUgTG9naXQgbW9kZWwgYW5kIGFzc2lnbiBpdCB0byBgTG9naXRfbW9kYFxuXG5cbiMgb2J0YWluIGEgcm9idXN0IHN1bW1hcnkgb2YgdGhlIG1vZGVsIGNvZWZmaWNpZW50c1xuXG5cbiMgcHJlZGljdCB0aGUgcHJvYmFiaWxpdHkgb2Ygc3Vydml2YWwgZm9yIHRoZSB0aHJlZSBwYXNzZW5nZXJzIiwic29sdXRpb24iOiIjIGZpdCB0aGUgTG9naXQgbW9kZWwsIGFzc2lnbiBpdCB0byBgTG9naXRfbW9kYFxuTG9naXRfbW9kIDwtIGdsbShTdXJ2aXZlZCB+IC4sIFxuZmFtaWx5ID0gYmlub21pYWwobGluayA9IFwibG9naXRcIiksIFxuZGF0YSA9IFRpdGFuaWNfMilcblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgbW9kZWwgY29lZmZpY2llbnRzXG5jb2VmdGVzdChMb2dpdF9tb2QsIHZjb3ZIQylcblxuIyBwcmVkaWN0IHRoZSBwcm9iYWJpbGl0eSBvZiBzdXJ2aXZhbCBmb3IgdGhlIGh5cG90aGVjaWFsIGluZGl2aWR1YWxzXG5wcmVkaWN0KExvZ2l0X21vZCwgbmV3ZGF0YSA9IHBhc3NlbmdlcnMsIHR5cGUgPSBcInJlc3BvbnNlXCIpIn0=

13.4.8 Câu 8: Dữ liệu sống sót sau sự kiện Titanic - Mô hình Probit

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGZpdCB0aGUgUHJvYml0IG1vZGVsIGFuZCBhc3NpZ24gaXQgdG8gYFByb2JpdF9tb2RgXG5cblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgbW9kZWwgY29lZmZpY2llbnRzXG5cblxuIyBwcmVkaWN0IHRoZSBwcm9iYWJpbGl0eSBvZiBzdXJ2aXZhbCBmb3IgdGhlIHRocmVlIHBhc3NlbmdlcnMiLCJzb2x1dGlvbiI6IiMgZml0IHRoZSBQcm9iaXQgbW9kZWwsIGFzc2lnbiBpdCB0byBgUHJvYml0X21vZGBcblByb2JpdF9tb2QgPC0gZ2xtKFN1cnZpdmVkIH4gLiwgZmFtaWx5ID0gYmlub21pYWwobGluayA9IFwicHJvYml0XCIpLCBkYXRhID0gVGl0YW5pY18yKVxuXG4jIG9idGFpbiBhIHJvYnVzdCBzdW1tYXJ5IG9mIHRoZSBtb2RlbCBjb2VmZmljaWVudHNcbmNvZWZ0ZXN0KFByb2JpdF9tb2QsIHZjb3ZIQylcblxuIyBwcmVkaWN0IHRoZSBwcm9iYWJpbGl0eSBvZiBzdXJ2aXZhbCBmb3IgdGhlIGh5cG90aGVjaWFsIGluZGl2aWR1YWxzXG5wcmVkaWN0KFByb2JpdF9tb2QsIG5ld2RhdGEgPSBwYXNzZW5nZXJzLCB0eXBlID0gXCJyZXNwb25zZVwiKSJ9

  1. Kleiber, C., & Zeileis, A. (2018). AER: Applied Econometrics with R