Chương 8 Mô hình hồi quy nhiều biến

Chương này sẽ giới thiệu mô hình hồi quy tuyến tính nhiều hơn một biến giải thích. Với sự xuất hiện của nhiều biến giải thích hơn, một số vấn đề mới sẽ phát sinh, đó là, đa cộng tuyến ( multicolinerity) và omitted variable bias (OVB). Trong chương này, vấn đề OVB được giới thiệu và hàm ý nhân quả của các hệ số hồi quy OLS.

Chương này sẽ sử dụng hai package của R:

  • AER (Kleiber & Zeileis, 2018)6
  • MASS (Ripley, 2018)7

8.1 Omitted variable bias

Ở các chương trước ta tìm hiểu quan hệ giữa hai biến Điểm thiQuy mô lớp học mà bỏ qua mối tương quan giữa Điểm thi và nhiều biến giải thích khác. Những mối quan hệ này đều được thể hiện trong hệ số error term với giả định là không tương quan với các biến giải thích. Điều này sẽ dẫn đến một sự thiên lệch về ước lượng, khi mà giá trị trung bình của các phân phối lấy mẫu của OLS estimator không còn bằng giá trị đúng của các hệ số hồi quy. Hiện tượng này gọi là omitted variable bias.

Có hai điều kiện cần phải thoả mãn để xảy ra OVB:

  • \(X\) tương quan với biến “vắng mặt”.
  • Biến “vắng mặt” là một biến giải thích đối với \(Y\).

Khi đó:

\[\hat{\beta}_1 \stackrel{p}{\rightarrow} \beta_1+\rho_{Xu} \frac{\sigma_u}{\sigma_X}\]

Trong ví dụ Điểm thiQuy mô lớp học, một biến giải thích nữa bị vắng mặt có thể là số lượng học sinh học tiếng Anh. Rõ ràng với một trường nước ngoài, việc phổ cập tiếng Anh cũng sẽ ảnh hưởng việc đọc hiểu các kiến thức ở những môn khác. Do đó, những học sinh ngoại quốc phải học tiếng Anh sẽ thực hiện bài thi kém hơn so với người bản xứ.

Mô hình được đề xuất trở thành:

\[\text{Điểm thi} = \beta_0 + \beta_1 \text{STR} +\beta_2 \text{PctEL}\]

trong đó \(\text{PctEL}\) là phần trăm học sinh học tiếng Anh.

Trong R, ta thực hiện việc input dữ liệu.

# load the data set
data(CASchools)   

# define variables
CASchools$STR <- CASchools$students/CASchools$teachers       
CASchools$score <- (CASchools$read + CASchools$math)/2

# compute correlations
cor(CASchools$STR, CASchools$score)
## [1] -0.2263627
cor(CASchools$STR, CASchools$english)
## [1] 0.1876424

Ta thấy \(\hat{\rho}_{STR,\text{Điểm thi}} = -0.2264\) là nguyên nhân cho việc loại bỏ \(\text{PctEL}\) dẫn đến thiên lệch âm đối với ước lượng \(\hat{\beta}_1\) vì điều này thể hiện \(\rho_{ Xu} <0\). Điều này dẫn đến ước lượng \(\hat{\beta}_1\) sẽ lớn hơn về mặt tuyệt đối. Nói cách khác, ước lượng OLS sẽ cho thấy, lớp quy mô nhỏ sẽ cải thiện được điểm thi, tuy nhiên sẽ làm quá mức độ ảnh hưởng, nghĩa là thực tế thì sẽ không cải thiện như mô hình thể hiện.

Bây giờ ta sẽ đưa thêm biến \(\text{PctEL}\) vào mô hình và so sánh mô hình trước và sau khi thêm. Bởi vì \(\hat{\rho}_{STR,\text{PctEL}}>0\) nên ta kỳ vọng dấu âm cho hệ số hồi quy của \(\text{PctEL}\), đồng thời trị tuyệt đối của hệ số biến \(STR\) sẽ giảm xuống do tác động được san sẻ cho biến \(\text{PctEL}\).

# estimate both regression models
mod <- lm(score ~ STR, data = CASchools) 
mult.mod <- lm(score ~ STR + english, data = CASchools)

# print the results to the console
mod
## 
## Call:
## lm(formula = score ~ STR, data = CASchools)
## 
## Coefficients:
## (Intercept)          STR  
##      698.93        -2.28
mult.mod
## 
## Call:
## lm(formula = score ~ STR + english, data = CASchools)
## 
## Coefficients:
## (Intercept)          STR      english  
##    686.0322      -1.1013      -0.6498

Với kết quả thể hiện, kỳ vọng của chúng ta là chính xác.

8.2 Mô hình hồi quy tuyến tính đa biến

Mô hình này mở rộng hơn so với mô hình hồi quy đơn biến, với cấu trúc như sau.

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

Ta sẽ không đi sâu vào việc làm thế nào để ước lượng mô hình này, mà làm thế nào để ước lượng nó bằng R. Tuy nhiên, đại ý cần nắm là phương pháp OLS vẫn được áp dụng cho trường hợp này.

Các đại lượng cần quan tâm:

summary(mult.mod)$coef
##                Estimate Std. Error    t value      Pr(>|t|)
## (Intercept) 686.0322445 7.41131160  92.565565 3.871327e-280
## STR          -1.1012956 0.38027827  -2.896026  3.978059e-03
## english      -0.6497768 0.03934254 -16.515882  1.657448e-47

Dữ liệu được thể hiện trên một biểu đồ 3 chiều, thay vì 2 chiều như mô hình hồi quy đơn biến. Do đó, sẽ không còn đường hồi quy nữa, mà ta có mặt phẳng hồi quy. Số chiều tăng lên khi ta tăng \(k\). Như vậy, dữ liệu sẽ được dàn trải trên không gian \(k+1\) chiều nếu như ta sử dụng \(k\) biến giải thích.

## The following object is masked _by_ .GlobalEnv:
## 
##     STR
## The following objects are masked from CASchools (pos = 4):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## The following objects are masked from CASchools (pos = 5):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, score, STR,
##     students, teachers
## The following objects are masked from CASchools (pos = 7):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## The following objects are masked from CASchools (pos = 8):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, score, STR,
##     students, teachers
## The following objects are masked from CASchools (pos = 10):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## The following objects are masked from CASchools (pos = 11):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, score, STR,
##     students, teachers
## The following objects are masked from CASchools (pos = 13):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## The following objects are masked from CASchools (pos = 19):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, score, STR,
##     students, teachers
## The following objects are masked from CASchools (pos = 24):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## Warning: 'surface' objects don't have these attributes: 'marker'
## Valid attributes include:
## 'type', 'visible', 'showlegend', 'legendgroup', 'opacity', 'name', 'uid', 'ids', 'customdata', 'selectedpoints', 'hoverinfo', 'hoverlabel', 'stream', 'transforms', 'z', 'x', 'y', 'text', 'surfacecolor', 'cauto', 'cmin', 'cmax', 'colorscale', 'autocolorscale', 'reversescale', 'showscale', 'colorbar', 'contours', 'hidesurface', 'lightposition', 'lighting', '_deprecated', 'xcalendar', 'ycalendar', 'zcalendar', 'scene', 'idssrc', 'customdatasrc', 'hoverinfosrc', 'zsrc', 'xsrc', 'ysrc', 'textsrc', 'surfacecolorsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

8.3 Đo lường tính hiệu quả của mô hình

Trong hồi quy đa biến, các thống kê thường dùng là \(SER\),\(R^2\)\(R^2\) hiệu chỉnh ( adjusted \(R^2\)) \[\begin{aligned} &SER = s_{\hat{u}} = \sqrt{s_{\hat{u}}^2} \\ &s_{\hat{u}}^2 = \frac{1}{n-k-1} SSR \\ &\bar{R}^2 = 1-\frac{n-1}{n-k-1} \frac{SSR}{TSS}\end{aligned}\] Trong R ta dùng hàm summary() để đồng thời lấy được cả ba giá trị này.

summary(mult.mod)
## 
## Call:
## lm(formula = score ~ STR + english, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
## STR          -1.10130    0.38028  -2.896  0.00398 ** 
## english      -0.64978    0.03934 -16.516  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

Ta cũng có thể tự tính tay các giá trị này theo công thức đã cho.

n <- nrow(CASchools)
k <- 2

y_mean <- mean(CASchools$score)

SSR <- sum(residuals(mult.mod)^2)
TSS <- sum((CASchools$score - y_mean)^2)
ESS <- sum((fitted(mult.mod)-y_mean)^2)

SER <- sqrt(1/(n-k-1)*SSR)
Rsq <- 1 - (SSR/TSS)
adj_Rsq <- 1- (n-1)/(n-k-1)* (SSR/TSS)

c("SER" = SER, "R2" = Rsq, "Adj.R2" = adj_Rsq)
##        SER         R2     Adj.R2 
## 14.4644831  0.4264315  0.4236805

Vậy khi ta thêm biến \(\text{PctEL}\) vào mô hình hồi quy thì mô hình có cải thiện hơn không? Câu trả lời là có: ta so sánh giá trị \(\bar{R}^2\) của cả hai mô hình có và không có biến \(\text{PctEL}\).

\(\bar{R}^2\) được dùng làm tiêu chuẩn đo lường mức độ hiệu quả của mô hình, tuy nhiên ta không nên lạm dụng chỉ tiêu này để thêm/bớt biến vào mô hình vô tội vạ, nhằm đạt được \(\bar{R}^2\) lớn nhất. Thay vào đó, ta nên quan tâm việc cải thiện quan hệ nhân quả mà \(\bar{R}^2\) không đo lường được. Phần này sẽ được nói ở các chương sau.

8.4 Giả định OLS

Ngoài các giả định OLS đối với mô hình hồi quy đơn giản, mô hình hồi quy đa biến bổ sung thêm một giả định trọng yếu: không có đa cộng tuyến hoàn hảo.

8.5 Đa cộng tuyến

Đa cộng tuyến nghĩa là hai hay nhiều biến giải thích trong mô hình tương quan mạnh với nhau. Nếu hệ số tương quan là hoàn hảo, ta có đa cộng tuyến hoàn hảo.

Trong khi đa cộng tuyến mạnh về mặt tổng quát không được hoan nghênh lắm vì sẽ làm phương sai của OLS estimator lớn hơn, đa cộng tuyến hoàn hảo sẽ làm OLS estimator không cách nào tính toán được.

Ví dụ sau sẽ cho ta góc nhìn về đa cộng tuyến hoàn hảo.

8.5.1 Ví dụ về đa cộng tuyến hoàn hảo

R sẽ xử lý như thế nào nếu ta ước lượng mô hình có đa cộng tuyến hoàn hảo?

lm() sẽ cho một cảnh báo ở dòng đầu tiên ở mục Coefficient trong kết quả hồi quy và sẽ bỏ qua biến giải thích có tương quan hoàn hảo. Trong ví dụ sau, ta thêm vào biến FracEL và tương quan hoàn hảo với biến english thông qua phép toán chia.

CASchools$FracEL <- CASchools$english / 100
mult.mod <- lm(score~STR+english+FracEL, data = CASchools)
summary(mult.mod)
## 
## Call:
## lm(formula = score ~ STR + english + FracEL, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients: (1 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
## STR          -1.10130    0.38028  -2.896  0.00398 ** 
## english      -0.64978    0.03934 -16.516  < 2e-16 ***
## FracEL             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

Tuy nhiên trong thực tế khó có trường hợp hiển nhiên thấy rõ các biến giải thích tương quan hoàn hảo với nhau. Ta xem xét hai ví dụ sau.

Ví dự thứ nhất, giả sử ta muốn phân tích ảnh hưởng quy mô lớp lên điểm thi bằng cách dùng biến dummy \(NS\) (not small) với quy tắc:

\[NS = \begin{cases} 0 \text{ nếu } STR<12 \\ 1 \end{cases}\]

CASchools$NS <- ifelse(CASchools$STR <12,0,1)
mult.mod <- lm(score~computer+english+NS, data = CASchools)
summary(mult.mod)
## 
## Call:
## lm(formula = score ~ computer + english + NS, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.492  -9.976  -0.778   8.761  43.798 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 663.704837   0.984259 674.319  < 2e-16 ***
## computer      0.005374   0.001670   3.218  0.00139 ** 
## english      -0.708947   0.040303 -17.591  < 2e-16 ***
## NS                  NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.43 on 417 degrees of freedom
## Multiple R-squared:  0.4291, Adjusted R-squared:  0.4263 
## F-statistic: 156.7 on 2 and 417 DF,  p-value: < 2.2e-16

Ta thấy việc thêm vào biến NS cho ra kết quả đa cộng tuyến hoàn hảo. Nguyên do là vì ta mắc phải lỗi sai logic.

table(CASchools$NS)
## 
##   1 
## 420

Ta thấy NS chỉ gồm một giá trị duy nhất là \(1\), điều này chứng tỏ không tồn tại quan sát \(STR<12\), do đó nó trùng với vector của hệ số \(\beta_0\). Như vậy đa cộng tuyến hoàn hảo giữa biến \(NS\)vector \(\begin{pmatrix}1 & \dots & 1\end{pmatrix}^T\) của hệ số \(\beta_0\).

Ví dụ thứ hai, ta rất hay bắt gặp trong tài chính, gọi là bẫy biến giả ( dummy variable trap). Bẫy này xảy ra khi các biến dummy được sử dụng làm biến giải thích. Chẳng hạn ta xét vị trí trường học ở phía Đông, Tây, Nam, Bắc bằng bộ biến giả sau.

\[\begin{aligned} &North_i = \begin{cases} 1 \text{ nếu ở phía Bắc} \\0 \end{cases} \\ &West_i = \begin{cases} 1 \text{ nếu ở phía Tây} \\0 \end{cases} \\ &South_i = \begin{cases} 1 \text{ nếu ở phía Nam} \\0 \end{cases} \\ &East_i = \begin{cases} 1 \text{ nếu ở phía Đông} \\0 \end{cases} \end{aligned}\]

Do các vùng là loại trừ lẫn nhau ( mutually exclusive) nên ta luôn có:

\[North_i + West_i + South_i + East_i = 1\]

Khi ta đưa cả bốn biến này vào mô hình thì sẽ xảy ra hiện tượng đa cộng tuyến hoàn hảo giữa chúng và intercept. Nhưng R không loại bỏ biến này.

set.seed(1)
CASchools$direction <- sample(c("West", "North", "South", "East"), 420,replace = T)
mult.mod <- lm(score~STR+english+direction, data = CASchools)
summary(mult.mod)
## 
## Call:
## lm(formula = score ~ STR + english + direction, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.018 -10.098  -0.556   9.304  42.596 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    685.67356    7.43308  92.246  < 2e-16 ***
## STR             -1.12246    0.38231  -2.936  0.00351 ** 
## english         -0.65096    0.03934 -16.549  < 2e-16 ***
## directionNorth   1.60680    1.92476   0.835  0.40431    
## directionSouth  -1.17013    2.07665  -0.563  0.57342    
## directionWest    2.44340    2.05191   1.191  0.23442    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.45 on 414 degrees of freedom
## Multiple R-squared:  0.4315, Adjusted R-squared:  0.4246 
## F-statistic: 62.85 on 5 and 414 DF,  p-value: < 2.2e-16

Ta thấy rằng, R đã phân tách biến direction thành các thành tố tuy nhiên trong mô hình sẽ bỏ đi thành tố cuối cùng directionEast. Một cách giải quyết khác là chúng ta bỏ đi hệ số intercept và bao gồm tất cả các biến dummy.

Tuy nhiên liệu mô hình có bỏ sót thông tin của directionEast? Không, không hề. R chỉ thay đổi cách giải thích các hệ số hồi quy từ tuyệt đối sang tương đối. Chẳng hạn, hệ số hồi quy của directionNorth nói rằng, về trung bình, điểm thi ở phía Bắc cao hơn \(1.61\) điểm so với khu vực phía Đông.

8.5.2 Đa cộng tuyến phi hoàn hảo

Đa cộng tuyến phi hoàn hảo không phải là vấn đề gì to tát, bởi lẽ OLS estimator cho phép cô lập ảnh hưởng của các biến giải thích có tương quan với nhau lên biến phụ thuộc. Tuy nhiên việc bỏ qua các mối liên hệ này sẽ ảnh hưởng đến kết quả hồi quy.

Giả sử ta có mô hình hồi quy

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + u_i\]

Đối với hồi quy đa biến, ta có:

\[\sigma_{\hat{\beta}_1}^2 = \frac{1}{n}\left(\frac{1}{1-\rho_{X_1,X_2}^2}\right) \frac{\sigma_u^2}{\sigma_{X_1}^2}\]

Khi \(\rho_{X_1,X_2} = 0\), việc thêm biến \(X_2\) vào sẽ không ảnh hưởng đến phương sai \(\hat{\beta}_1\). Khi tương quan giữa hai biến càng mạnh, phương sai này càng lớn. Khi tăng quy mô mẫu, giá trị phương sai sẽ giảm. Như vậy, đa cộng tuyến phi hoàn hảo sẽ làm tăng phương sai của một hoặc nhiều estimator của hệ số hồi quy. Khi quy mô mẫu nhỏ, ta phải đối diện với quyết định: chấp nhận phương sai lớn hay sử dụng ít biến giải thích. Đây là hiện tượng bias-variance trade-off.

8.6 Phân phối của các hệ số hồi quy

Khi các giả định OLS được thoả mãn, phân phối của các OLS estimator là phân phối Normal đa biến.

Trong R, ta khảo sát phân phối này như sau.

Giả sử tiến trình tạo lập dữ liệu (Data gernation process DGP) như sau:

\[Y_i = 5+2.5 X_{1i} + 3X_{2i}+ u _i\]

với \(X_i = (X_{1i},X_{2i}) \sim N\left[\begin{pmatrix} 0\\0 \end{pmatrix}, \begin{pmatrix} 10 & 2.5\\2.5 & 10 \end{pmatrix} \right]\)\(u_i \sim N(0,5)\).

Ta tạo \(10000\) mẫu ngẫu nhiên quy mô \(50\) để ước lượng mô hình:

\[Y_i = \beta_0+\beta_1 X_{1i} + \beta_2 X_{2i}+ u _i\]

Bước cuối cùng khảo sát phân phối của \(\hat{\beta}_1, \hat{\beta}_2\).

# load packages
library(MASS)
library(mvtnorm)

# set sample size
n <- 50

# initialize vector of coefficients
coefs <- cbind("hat_beta_1" = numeric(10000), "hat_beta_2" = numeric(10000))

# set seed for reproducibility
set.seed(1)

# loop sampling and estimation
for (i in 1:10000) {
  
  X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
  u <- rnorm(n, sd = 5)
  Y <- 5 + 2.5 * X[, 1] + 3 * X[, 2] + u
  coefs[i,] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
  
}

# compute density estimate
kde <- kde2d(coefs[, 1], coefs[, 2])

# plot density estimate
persp(kde, 
      theta = 310, 
      phi = 30, 
      xlab = "beta_1", 
      ylab = "beta_2", 
      zlab = "Est. Density")

Đồng thời, ta ước lượng hệ số tương quan giữa các ước lượng.

# estimate the correlation between estimators
cor(coefs[, 1], coefs[, 2])
## [1] -0.2503028

Ta thấy rằng trong quá trình tạo lập dữ liệu, ta đã để cho \(X_1\)\(X_2\) có tương quan với nhau với hiệp phương sai là \(2.5\). Tương quan các biến giải thích luôn luôn kéo theo tương quan của các hệ số hồi quy đi kèm chúng.

8.7 Bài tập lập trình

8.7.1 Bài 1: Dữ liệu giá nhà ở Boston

Dữ liệu giá nhà ở Boston bao gồm 506 quan sát.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCBib3RoIHBhY2thZ2VzIGFuZCBsb2FkIHRoZSBkYXRhIHNldFxuZGF0YSgnQm9zdG9uJylcblxuIyBvYnRhaW4gYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YSBzZXRcblxuXG4jIGVzdGltYXRlIHRoZSBzaW1wbGUgcmVncmVzc2lvbiBtb2RlbFxuXG5cbiMgcHJpbnQgdGhlIHN1bW1hcnkgdG8gdGhlIGNvbnNvbGUiLCJzb2x1dGlvbiI6IiMgYXR0YWNoIGJvdGggcGFja2FnZXMgYW5kIGxvYWQgdGhlIGRhdGEgc2V0XG5kYXRhKCdCb3N0b24nKVxuIyBvYnRhaW4gYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YSBzZXRcbnN1bW1hcnkoQm9zdG9uKVxuIyBvclxuc3RyKEJvc3RvbilcbiMgb3JcbmhlYWQoQm9zdG9uKVxuIyBlc3RpbWF0ZSB0aGUgc2ltcGxlIHJlZ3Jlc3Npb24gbW9kZWxcbmJoX21vZCA8LSBsbShtZWR2IH4gbHN0YXQsIGRhdGEgPSBCb3N0b24pXG4jIHByaW50IHRoZSBzdW1tYXJ5IHRvIHRoZSBjb25zb2xlXG5jb2VmdGVzdChiaF9tb2QsIHZjb3YuID0gdmNvdkhDKSJ9

8.7.2 Bài 2: Hồi quy đa biến giá nhà ở Boston 1

Thực hiện hồi quy:

\[medv_i = \beta_0 + \beta_1 lstat_i + \beta_2 age_i + \beta_3 crim_i + u_i\]

trong đó \(medv\) là trung vị giá nhà trong một quận, \(age\) là tuổi toà nhà, \(crim\) là tỷ lệ tội phạm, và \(lstat\) là tỷ lệ cá nhân có thu nhập thấp.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbmR1Y3QgdGhlIHJlZ3Jlc3Npb25cbm1vZFxuXG4jIG9idGFpbiBhIHJvYnVzdCBjb2VmZmljaWVudCBzdW1tYXJ5XG5saWJyYXJ5KGxtdGVzdClcblxuIyAgY29lZmZpY2llbnRzIG9mIGRldGVybWluYXRpb25cblIyX3VucmVzIiwic29sdXRpb24iOiIjIGNvbmR1Y3QgdGhlIHJlZ3Jlc3Npb25cbm1vZCA8LSBsbShtZWR2IH4gbHN0YXQgKyBjcmltICsgYWdlLCBkYXRhID0gQm9zdG9uKVxuXG4jIG9idGFpbiBhIHJvYnVzdCBjb2VmZmljaWVudCBzdW1tYXJ5XG5saWJyYXJ5KGxtdGVzdClcbmNvZWZ0ZXN0KG1vZCwgdmNvdi4gPSB2Y292SEMpXG5cbiMgY29lZmZpY2llbnRzIG9mIGRldGVybWluYXRpb25cblIyX3VucmVzIDwtIHN1bW1hcnkobW9kKSRyLnNxdWFyZWRcblIyX3VucmVzIn0=

8.7.3 Bài 3: Hồi quy đa biến giá nhà ở Boston 2

Tính toán \(\bar{R}^2\) của mô hình ở trên.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiJzdW1tYXJ5KCkiLCJzb2x1dGlvbiI6InN1bW1hcnkobW9kKSRhZGouci5zcXVhcmVkIn0=

8.7.4 Bài 4: Hồi quy đa biến giá nhà ở Boston 3

Thực hiện hồi quy đa biến dùng tất cả các biến trong bộ dữ liệu để giải thích cho medv.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIHJ1biBhIHJlZ3Jlc3Npb24gb2YgbWVkdiBvbiBhbGwgcmVtYWluaW5nIHZhcmlhYmxlcyBpbiB0aGUgQm9zdG9uIGRhdGEgc2V0XG5cblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgY29lZmZpY2llbnRzXG5cblxuIyB3aGF0IGlzIHRoZSBhZGp1c3RlZCBSXjIgb2YgdGhlIG1vZGVsPyIsInNvbHV0aW9uIjoiIyBydW4gYSByZWdyZXNzaW9uIG9mIG1lZHYgb24gYWxsIHJlbWFpbmluZyB2YXJpYWJsZXMgaW4gdGhlIEJvc3RvbiBkYXRhIHNldFxuZnVsbF9tb2QgPC0gbG0obWVkdiB+LiwgZGF0YSA9IEJvc3RvbilcblxuIyBvYnRhaW4gYSByb2J1c3Qgc3VtbWFyeSBvZiB0aGUgY29lZmZpY2llbnRzXG5jb2VmdGVzdChmdWxsX21vZCwgdmNvdi4gPSB2Y292SEMsIHR5cGUgPSBcIkhDMVwiKVxuXG4jIHdoYXQgaXMgdGhlIGFkanVzdGVkIFJeMiBvZiB0aGUgbW9kZWw/XG5zdW1tYXJ5KGZ1bGxfbW9kKSRhZGouci5zcXVhcmVkIn0=

8.7.5 Bài 5: Hồi quy đa biến giá nhà ở Boston 4

Ta có thể cải thiện mô hình bằng cách bỏ bớt biến. Ta cần ước lượng nhiều mô hình, mỗi lần ước lượng ta sẽ bỏ bớt một biến và so sánh \(R^2\) hiệu chỉnh giữa mô hình mới và mô hình trước khi bỏ biến.

Trong R ta làm như sau.

  1. Lưu trữ mô hình có đầy đủ biến ở bài tập 4 vào full_mod.
  2. Ước lượng mô hình mới mod_new trong đó đã bỏ bớt biến, chẳng hạn, lstat.
  3. Lấy kết quả \(\bar{R}^2\) của mô hình.
  4. Lặp lại các bước trên và lưu trữ mô hình có \(\bar{R}^2\) cao nhất, lưu trữ mô hình này vào better_mod.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIFN0ZXAgMVxuXG4jIFN0ZXAgMi0zOiB1c2UgZm9yIGxvb3BcblxuIyBTdGVwIDQiLCJzb2x1dGlvbiI6IiMgdGhpcyBzb2x1dGlvbiBpcyBhIGJpdCB0ZWNobmljYWwgYnV0IGVmZmljaWVudFxuXG4jIGxvb3AgZXN0aW1hdGlvbiBvZiBtb2RlbHNcbmwgPC0gbGlzdCgpXG5mb3IgKGkgaW4gMToxMykge1xuICBkIDwtIEJvc3RvblssIC1pXVxuICAjIHNhdmUgZWFjaCBhZGouIFJeMiBhcyBhIGxpc3QgZW50cnkgaW4gbFxuICBsW1tpXV0gPC0gc3VtbWFyeShsbShtZWR2IH4uLCBkYXRhPWQpKSRhZGouci5zcXVhcmVkIFxufVxuXG4jIGFzc2lnbiB2YXJpYWJsZSBuYW1lcyB0byB0aGUgbGlzdCBlbnRyaWVzXG5uYW1lcyhsKSA8LSBuYW1lcyhCb3N0b25bLCAxOjEzXSkgXG5cbiMgc2VsZWN0IHRoZSB2YXJpYWJsZSB3aG9zZSBvbWlzc2lvbiBsZWFkcyB0byB0aGUgaGlnaGVzdCBpbXByb3ZlbWVudCBpbiBhZGouIFJeMlxud2hpY2gubWF4KGwpICMgN3RoIGNvbHVtbiB0aGlzIGlzIFwiYWdlXCJcblxuIyBoZW5jZSBhIG1vZGVsIHRoYXQgZml0cyB0aGUgZGF0YSBiZXR0ZXIgaXNcbmJldHRlcl9tb2RlbCA8LSBsbShtZWR2IH4uLCBkYXRhID0gQm9zdG9uWywgLTddKSJ9

8.8 Gợi ý chứng minh công thức phương sai sai số hệ số hồi quy

Từ kết quả thừa hưởng từ chứng minh công thức ở chương 7, ta có

\[\hat{\beta}_1 = \beta_1 + \frac{\sum_i (X_i - \bar{X})u_i}{\sum_i (X_i - \bar{X})^2}\]

Ta cũng biết rằng:

\[\begin{aligned} &\frac{1}{n}\sum_i (X_i - \bar{X})^2 \stackrel{p}{\rightarrow} \sigma^2_X \\ &\frac{1}{n}\sum_i (X_i - \bar{X})u_i \stackrel{p}{\rightarrow} \text{Cov}(u_i,X_i) = \rho_{Xu}\sigma_u\sigma_X \end{aligned}\]

Từ đó:

\[\hat{\beta}_1 \stackrel{p}{\rightarrow} \beta_1 + \rho_{Xu}\frac{\sigma_u}{\sigma_X}\]


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

  2. Ripley, B. (2018). MASS: Support Functions and Datasets for Venables and Ripley’s MASS.