Chương 12 Hồi quy dữ liệu bảng

Hồi quy sử dụng dữ liệu bảng có thể giảm thiểu omitted variable bias khi không có thông tin các biến tương quan với các biến được giải thích và biến giải thích và nếu những biến này là hằng số về chiều thời gian và các chủ thể đối tượng.

Chương này mô tả các chủ đề sau:

  • ký hiệu liên quan đến dữ liệu bảng.
  • hồi quy tác động cố định ( fixed effect regression) bằng cách cố định chiều thời gian, hoặc chiều không gian, hoặc cả hai.
  • tính sai số chuẩn đối với hồi quy tác động cố định.

Ta ứng dụng bộ dữ liệu Fatalities từ package AER liên quan đến các quan sát tai nạn giao thông Mỹ từ 1982 đến 1988. Ta khảo sát tác động của thuế bia và luật lái xe khi say.

Trong chương này ta sử dụng các package sau.

12.1 Dữ liệu bảng

Khác với dữ liệu cross-section data khi ta quan sát \(n\) chủ thể, dữ liệu bảng có \(n\) chủ thể tại \(T \ge 2\) thời điểm khác nhau. Ta ký hiệu: \((X_{it}, Y_{it})\).

Ta khảo sát dữ liệu Fatalities trong R như sau.

# load the package and the dataset
library(AER)
data(Fatalities)
# obtain the dimension and inspect the structure
is.data.frame(Fatalities)
## [1] TRUE
dim(Fatalities)
## [1] 336  34
str(Fatalities)
## 'data.frame':    336 obs. of  34 variables:
##  $ state       : Factor w/ 48 levels "al","az","ar",..: 1 1 1 1 1 1 1 2 2 2 ...
##  $ year        : Factor w/ 7 levels "1982","1983",..: 1 2 3 4 5 6 7 1 2 3 ...
##  $ spirits     : num  1.37 1.36 1.32 1.28 1.23 ...
##  $ unemp       : num  14.4 13.7 11.1 8.9 9.8 ...
##  $ income      : num  10544 10733 11109 11333 11662 ...
##  $ emppop      : num  50.7 52.1 54.2 55.3 56.5 ...
##  $ beertax     : num  1.54 1.79 1.71 1.65 1.61 ...
##  $ baptist     : num  30.4 30.3 30.3 30.3 30.3 ...
##  $ mormon      : num  0.328 0.343 0.359 0.376 0.393 ...
##  $ drinkage    : num  19 19 19 19.7 21 ...
##  $ dry         : num  25 23 24 23.6 23.5 ...
##  $ youngdrivers: num  0.212 0.211 0.211 0.211 0.213 ...
##  $ miles       : num  7234 7836 8263 8727 8953 ...
##  $ breath      : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ jail        : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 2 2 2 ...
##  $ service     : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 2 2 2 ...
##  $ fatal       : int  839 930 932 882 1081 1110 1023 724 675 869 ...
##  $ nfatal      : int  146 154 165 146 172 181 139 131 112 149 ...
##  $ sfatal      : int  99 98 94 98 119 114 89 76 60 81 ...
##  $ fatal1517   : int  53 71 49 66 82 94 66 40 40 51 ...
##  $ nfatal1517  : int  9 8 7 9 10 11 8 7 7 8 ...
##  $ fatal1820   : int  99 108 103 100 120 127 105 81 83 118 ...
##  $ nfatal1820  : int  34 26 25 23 23 31 24 16 19 34 ...
##  $ fatal2124   : int  120 124 118 114 119 138 123 96 80 123 ...
##  $ nfatal2124  : int  32 35 34 45 29 30 25 36 17 33 ...
##  $ afatal      : num  309 342 305 277 361 ...
##  $ pop         : num  3942002 3960008 3988992 4021008 4049994 ...
##  $ pop1517     : num  209000 202000 197000 195000 204000 ...
##  $ pop1820     : num  221553 219125 216724 214349 212000 ...
##  $ pop2124     : num  290000 290000 288000 284000 263000 ...
##  $ milestot    : num  28516 31032 32961 35091 36259 ...
##  $ unempus     : num  9.7 9.6 7.5 7.2 7 ...
##  $ emppopus    : num  57.8 57.9 59.5 60.1 60.7 ...
##  $ gsp         : num  -0.0221 0.0466 0.0628 0.0275 0.0321 ...
# list the first few observations
head(Fatalities)
##   state year spirits unemp   income   emppop  beertax baptist  mormon
## 1    al 1982    1.37  14.4 10544.15 50.69204 1.539379 30.3557 0.32829
## 2    al 1983    1.36  13.7 10732.80 52.14703 1.788991 30.3336 0.34341
## 3    al 1984    1.32  11.1 11108.79 54.16809 1.714286 30.3115 0.35924
## 4    al 1985    1.28   8.9 11332.63 55.27114 1.652542 30.2895 0.37579
## 5    al 1986    1.23   9.8 11661.51 56.51450 1.609907 30.2674 0.39311
## 6    al 1987    1.18   7.8 11944.00 57.50988 1.560000 30.2453 0.41123
##   drinkage     dry youngdrivers    miles breath jail service fatal nfatal
## 1    19.00 25.0063     0.211572 7233.887     no   no      no   839    146
## 2    19.00 22.9942     0.210768 7836.348     no   no      no   930    154
## 3    19.00 24.0426     0.211484 8262.990     no   no      no   932    165
## 4    19.67 23.6339     0.211140 8726.917     no   no      no   882    146
## 5    21.00 23.4647     0.213400 8952.854     no   no      no  1081    172
## 6    21.00 23.7924     0.215527 9166.302     no   no      no  1110    181
##   sfatal fatal1517 nfatal1517 fatal1820 nfatal1820 fatal2124 nfatal2124
## 1     99        53          9        99         34       120         32
## 2     98        71          8       108         26       124         35
## 3     94        49          7       103         25       118         34
## 4     98        66          9       100         23       114         45
## 5    119        82         10       120         23       119         29
## 6    114        94         11       127         31       138         30
##    afatal     pop  pop1517  pop1820  pop2124 milestot unempus emppopus
## 1 309.438 3942002 208999.6 221553.4 290000.1    28516     9.7     57.8
## 2 341.834 3960008 202000.1 219125.5 290000.2    31032     9.6     57.9
## 3 304.872 3988992 197000.0 216724.1 288000.2    32961     7.5     59.5
## 4 276.742 4021008 194999.7 214349.0 284000.3    35091     7.2     60.1
## 5 360.716 4049994 203999.9 212000.0 263000.3    36259     7.0     60.7
## 6 368.421 4082999 204999.8 208998.5 258999.8    37426     6.2     61.5
##           gsp
## 1 -0.02212476
## 2  0.04655825
## 3  0.06279784
## 4  0.02748997
## 5  0.03214295
## 6  0.04897637
# summarize the variables 'state' and 'year'
summary(Fatalities[, c(1, 2)])
##      state       year   
##  al     :  7   1982:48  
##  az     :  7   1983:48  
##  ar     :  7   1984:48  
##  ca     :  7   1985:48  
##  co     :  7   1986:48  
##  ct     :  7   1987:48  
##  (Other):294   1988:48

Ta thấy bộ dữ liệu bao gồm \(336\) quan sát và \(34\) biến. Chú ý rằng state là biến nhân tố với \(48\) bậc, ứng với \(48\) tiểu bang kề nhau của Mỹ. Biến year cũng là biến nhân tố với \(7\) bậc xác định thời gian quan sát. Do đó ta có \(7\times 48 = 336\) quan sát. Bởi vì không có dữ liệu trống, nên bộ dữ liệu xác định cân bằng ( balanced), trong trường hợp có dữ liệu trống, bộ dữ liệu trở thành phi cân bằng ( unbalanced).

12.1.1 Ví dụ thuế bia và tai nạn giao thông

Ta ước lượng hồi quy đơn giản giữa thuế bia (được điều chỉnh so với đồng đô la 1988) và tỷ lệ tai nạn giao thông.

# define the fatality rate
Fatalities$fatal_rate <- Fatalities$fatal / Fatalities$pop * 10000

# subset the data
Fatalities1982 <- subset(Fatalities, year == "1982")
Fatalities1988 <- subset(Fatalities, year == "1988")

# estimate simple regression models using 1982 and 1988 data
fatal1982_mod <- lm(fatal_rate ~ beertax, data = Fatalities1982)
fatal1988_mod <- lm(fatal_rate ~ beertax, data = Fatalities1988)

coeftest(fatal1982_mod, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.01038    0.14957 13.4408   <2e-16 ***
## beertax      0.14846    0.13261  1.1196   0.2687    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fatal1988_mod, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.85907    0.11461 16.2205 < 2.2e-16 ***
## beertax      0.43875    0.12786  3.4314  0.001279 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mô hình hồi quy như sau

\[\begin{aligned} &\hat{FatalityRate} = 2.01 + 0.15 \times BeerTax \text{ (1982)} \\ &\hat{FatalityRate} = 1.86 + 0.44 \times BeerTax \text{ (1988)} \end{aligned}\]

# plot the observations and add the estimated regression line for 1982 data
plot(x = Fatalities1982$beertax, 
     y = Fatalities1982$fatal_rate, 
     xlab = "Beer tax (in 1988 dollars)",
     ylab = "Fatality rate (fatalities per 10000)",
     main = "Traffic Fatality Rates and Beer Taxes in 1982",
     ylim = c(0, 4.5),
     pch = 20, 
     col = "steelblue")

abline(fatal1982_mod, lwd = 1.5)

# plot observations and add estimated regression line for 1988 data
plot(x = Fatalities1988$beertax, 
     y = Fatalities1988$fatal_rate, 
     xlab = "Beer tax (in 1988 dollars)",
     ylab = "Fatality rate (fatalities per 10000)",
     main = "Traffic Fatality Rates and Beer Taxes in 1988",
     ylim = c(0, 4.5),
     pch = 20, 
     col = "steelblue")

abline(fatal1988_mod, lwd = 1.5)

Trong cả hai hình, mỗi điểm tượng trưng cho quan sát của thuế và tỷ lệ tai nạn cho một bang đối với một năm tương ứng. Hồi quy cho thấy có tương quan dương giữa thuế bia và tỷ lệ tai nạn cho cả hai năm quan sát. Tuy nhiên hệ số hồi quy cho năm 1988 gấp ba lần hệ số năm 1982. Điều này là ngược với kỳ vọng của chúng ta: thuế bia rượu sẽ làm giảm tỷ lệ tai nạn giao thông. Điều này có thể xuất phát từ omitted variable bias, bởi vì cả hai mô hình đều không có biến giải thích về điều kiện kinh tế. Điều này có thể được sửa chữa nhờ vào hồi quy đa biến.

Tuy nhiên, phương pháp hồi quy đa biến không tính được những nhân tố không quan sát được do sự khác biệt về địa lý. Và hồi quy dữ liệu bảng cho phép chúng ta cố định những nhân tố này, để loại bỏ sự khác biệt không gian hoặc/ và thời gian quan sát.

12.2 Dữ liệu bảng với hai chiều thời gian

Giả sử ta chỉ có hai khung thời gian quan sát, \(T=2\): \(t=1982,1988\). Điều này cho phép ta phân tích sự khác biệt giữa sự thay đổi của tỷ lệ tai nạn từ năm 1982 đến 1988. Mô hình hồi quy tổng thể như sau.

\[Fatality_{it} = \beta_0 + \beta_1 BeerTax_{it} + \beta_2 Z_i + u_{it}\]

trong đó \(Z_i\) là yếu tố xác định tiểu bang và là hằng số theo thời gian. Ta có thể loại bỏ biến này bằng cách:

\[Fatality_{i1988} - Fatality_{i1982} = \beta_1(BeerTax_{i1988}-BeerTax_{i1982}) +u_{i1988} - u_{i1982}\]

Mô hình này cho ra ước lượng vững \(\beta_1\) bị thiên lệch đôi chút do bỏ đi yếu tố địa lý \(Z_i\).

# compute the differences 
diff_fatal_rate <- Fatalities1988$fatal_rate - Fatalities1982$fatal_rate
diff_beertax <- Fatalities1988$beertax - Fatalities1982$beertax

# estimate a regression using differenced data
fatal_diff_mod <- lm(diff_fatal_rate ~ diff_beertax)

coeftest(fatal_diff_mod, vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.072037   0.065355 -1.1022 0.276091   
## diff_beertax -1.040973   0.355006 -2.9323 0.005229 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plot the differenced data
plot(x = diff_beertax, 
     y = diff_fatal_rate, 
     xlab = "Change in beer tax (in 1988 dollars)",
     ylab = "Change in fatality rate (fatalities per 10000)",
     main = "Changes in Traffic Fatality Rates and Beer Taxes in 1982-1988",
     xlim = c(-0.6, 0.6),
     ylim = c(-1.5, 1),
     pch = 20, 
     col = "steelblue")

# add the regression line to plot
abline(fatal_diff_mod, lwd = 1.5)

Kết quả cho thấy hệ số hồi quy đối với BeerTax đã âm và khác \(0\) thống kê tại \(5\%\). Ý nghĩa hệ số là nếu tăng thuế lên \(\$1\) thì tỷ lệ tai nạn giảm \(1.04\) mỗi \(10000\) người.

Một lần nữa, kết quả này có thể là hậu quả của việc bỏ qua các yếu tố trong hồi quy một năm mà có ảnh hưởng đến tỷ lệ tử vong và tương quan với thuế bia và thay đổi theo thời gian. Thông điệp ở đây là ta cần cẩn thận và kiểm soát các yếu tố như vậy trước khi đưa ra kết luận về tác động của việc tăng thuế bia.

12.3 Hồi quy tác động cố định

Xét mô hình hồi quy

\[Y_{it} = \alpha_{it} + \beta_1 X_{it} + \beta_2 Z_i + u_{it}\]

trong đó \(Z_i\) là biến số không đồng nhất phi biến thời gian không quan sát. Mục tiêu là ước lượng \(\beta_1\). Đặt \(\alpha_i = \beta_0 + \beta_2 Z_i\) ta sẽ có mô hình

\[Y_{it} = \alpha_i + \beta_1 X_{it} + u_{it}\]

Với mỗi một intercept nhất định, mỗi giá trị này có thể hiểu là tác động cố định của mỗi chủ thể \(i\). Khi đó mô hình này được gọi là mô hình hồi quy tác động cố định. Sự biến thiên của \(\alpha_i\) xuất phát từ sự biến thiên của \(Z_i\). Mô hình hồi quy tác động cố định có thể được viết lại công thức chứa \(n-1\) biến giả và hằng số:

\[Y_{it} = \beta_0 + \beta_1 X_{it} + \gamma_2 D2_i + \gamma_3 D3_i + \dots + \gamma_n Dn_i + u_{it}\]

Mô hình này có \(n\) hệ số intercept khác nhau, mỗi intercept đại diện cho một chủ thể.

12.3.1 Ước lượng và suy luận

Các phần mềm thống kê sử dụng thuật toán entity-demeaned OLS để ước lượng mô hình có \(k+n\) biến giải thích. Từ mô hình hồi quy tác động cố định rút gọn, ta lấy trung bình hai vế:

\[\begin{aligned} \frac{1}{n} \sum_i Y_{it} &= \beta_1 \frac{1}{n} \sum_i X_{it} + \frac{1}{n}\sum_i \alpha _i + \frac{1}{n} \sum_i u_{it} \\ \bar{Y}_i &= \beta_1 \bar{X}_i + \alpha_i +\bar{u}_{it} \end{aligned}\]

Suy ra:

\[\begin{aligned} Y_{it} - \bar{Y}_i &= \beta_1 (X_{it} -\bar{X}_i ) + (u_{it}-\bar{u}_i) \\ \tilde{Y}_{it} &= \beta_1 \tilde{X}_{it} + \tilde{u}_{it} \end{aligned}\]

Trong mô hình này, ước lượng OLS của \(\beta_1\) bằng với ước lượng của mô hình biến giả mà không cần ước lượng \(n-1\) biến giả và intercept.

Ta kết luận có hai cách ước lượng \(\beta_1\) trong hồi quy tác động cố định:

  • OLS của mô hình biến giả.
  • OLS của mô hình loại bỏ trung bình.

Miễn là giả định OLS được thoả mãn, phân phối lấy mẫu của OLS estimator sẽ tuân theo Normal với quy mô lớn.

12.3.2 Áp dụng vào dữ liệu tai nạn giao thông

Mô hình hồi quy tác động cố định đơn giản được đề xuất như sau

\[FatalityRate_{it} = \beta_1 BeerTax_{it} + StateFixedEffects + u_{it}\]

fatal_fe_lm_mod <- lm(fatal_rate ~ beertax + state - 1, data = Fatalities)
fatal_fe_lm_mod
## 
## Call:
## lm(formula = fatal_rate ~ beertax + state - 1, data = Fatalities)
## 
## Coefficients:
## beertax  stateal  stateaz  statear  stateca  stateco  statect  statede  
## -0.6559   3.4776   2.9099   2.8227   1.9682   1.9933   1.6154   2.1700  
## statefl  statega  stateid  stateil  statein  stateia  stateks  stateky  
##  3.2095   4.0022   2.8086   1.5160   2.0161   1.9337   2.2544   2.2601  
## statela  stateme  statemd  statema  statemi  statemn  statems  statemo  
##  2.6305   2.3697   1.7712   1.3679   1.9931   1.5804   3.4486   2.1814  
## statemt  statene  statenv  statenh  statenj  statenm  stateny  statenc  
##  3.1172   1.9555   2.8769   2.2232   1.3719   3.9040   1.2910   3.1872  
## statend  stateoh  stateok  stateor  statepa  stateri  statesc  statesd  
##  1.8542   1.8032   2.9326   2.3096   1.7102   1.2126   4.0348   2.4739  
## statetn  statetx  stateut  statevt  stateva  statewa  statewv  statewi  
##  2.6020   2.5602   2.3137   2.5116   2.1874   1.8181   2.5809   1.7184  
## statewy  
##  3.2491

Như đã trình bày, ta có thể ước lượng \(\beta_1\) bằng cách dùng dữ liệu loại bỏ trung bình, tức thực hiện hồi quy

\[\tilde{FatalityRate} = \beta_1 \tilde{BeerTax}_{it} + u_{it}\]

# obtain demeaned data
Fatalities_demeaned <- with(Fatalities,
            data.frame(fatal_rate = fatal_rate - ave(fatal_rate, state),
            beertax = beertax - ave(beertax, state)))

# estimate the regression
summary(lm(fatal_rate ~ beertax - 1, data = Fatalities_demeaned))
## 
## Call:
## lm(formula = fatal_rate ~ beertax - 1, data = Fatalities_demeaned)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58696 -0.08284 -0.00127  0.07955  0.89780 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## beertax  -0.6559     0.1739  -3.772 0.000191 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 335 degrees of freedom
## Multiple R-squared:  0.04074,    Adjusted R-squared:  0.03788 
## F-statistic: 14.23 on 1 and 335 DF,  p-value: 0.0001913

Một cách khác trong R là sử dụng hàm plm(). Như với hàm lm() ta cũng phải định hình quy cách công thức và dữ liệu được dùng cho plm(). Tuy nhiên ta sẽ phải chỉ ra biến chủ thể và biến thời gian. Trong Fatalities biến chủ thể là state trong khi biến thời gian là year. Bởi vì estimator cho tác động cố định cũng được gọi là within estimator, ta cài đặt mô hình model = "within".

# install and load the 'plm' package
## install.packages("plm")
library(plm)
# estimate the fixed effects regression with plm()
fatal_fe_mod <- plm(fatal_rate ~ beertax, 
                    data = Fatalities,
                    index = c("state", "year"), 
                    model = "within")

# print summary using robust standard errors
coeftest(fatal_fe_mod, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##         Estimate Std. Error t value Pr(>|t|)  
## beertax -0.65587    0.28880  -2.271  0.02388 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hệ số ước lượng theo R chính là \(-0.6559\). plm() dùng thuật toán loại trừ trung bình cho nên sẽ không khai báo về biến giả.

Ta thấy rằng, việc cố định không gian cũng chưa loại bỏ hết thiên lệch trong ước lượng, vì hệ số tác động cũng khá bé. Do đó ta nghi ngờ việc một tác động nào đó đã bị bỏ qua và nó thay đổi theo thời gian, dẫn đến việc vẫn còn thiên lệch trong ước lượng.

12.4 Hồi quy cố định thời gian

Nếu chỉ có tác động cố định thời gian, mô hình hồi quy như sau

\[Y_{it} = \beta_0 + \beta_1 X_{it} + \delta_2B2_t+\dots +\delta_TBT_t + u_{it}\]

trong đó có \(T-1\) biến giả.

Kết hợp cả hai tác động, ta được mô hình dưới đây.

\[Y_{it} = \beta_0 + \beta_1X_{it} + \gamma_2D2_i + \dots + \gamma_nDT_i + \delta_2B2_t + \dots + \delta_T BT_i + u_{it}\]

Mô hình kết hợp cho phép ta loại bỏ sự thiên lệch từ những thay đổi không quan sát được qua thời gian nhưng cố định qua chủ thể và từ những nhân tố khác nhau bởi chủ thể nhưng cố định qua thời gian. Mô hình này cũng được thực hiện nhờ vào thuật toán OLS.

Trong R, ta thực hiện hồi quy mô hình kết hợp như sau.

\[FatalityRate_{it} = \beta_0 BeerTax_{it} + StateEffects + TimeFixedEffects + u_{it}\]

# estimate a combined time and entity fixed effects regression model

# via lm()
fatal_tefe_lm_mod <- lm(fatal_rate ~ beertax + state + year - 1, data = Fatalities)
fatal_tefe_lm_mod
## 
## Call:
## lm(formula = fatal_rate ~ beertax + state + year - 1, data = Fatalities)
## 
## Coefficients:
##  beertax   stateal   stateaz   statear   stateca   stateco   statect  
## -0.63998   3.51137   2.96451   2.87284   2.02618   2.04984   1.67125  
##  statede   statefl   statega   stateid   stateil   statein   stateia  
##  2.22711   3.25132   4.02300   2.86242   1.57287   2.07123   1.98709  
##  stateks   stateky   statela   stateme   statemd   statema   statemi  
##  2.30707   2.31659   2.67772   2.41713   1.82731   1.42335   2.04488  
##  statemn   statems   statemo   statemt   statene   statenv   statenh  
##  1.63488   3.49146   2.23598   3.17160   2.00846   2.93322   2.27245  
##  statenj   statenm   stateny   statenc   statend   stateoh   stateok  
##  1.43016   3.95748   1.34849   3.22630   1.90762   1.85664   2.97776  
##  stateor   statepa   stateri   statesc   statesd   statetn   statetx  
##  2.36597   1.76563   1.26964   4.06496   2.52317   2.65670   2.61282  
##  stateut   statevt   stateva   statewa   statewv   statewi   statewy  
##  2.36165   2.56100   2.23618   1.87424   2.63364   1.77545   3.30791  
## year1983  year1984  year1985  year1986  year1987  year1988  
## -0.07990  -0.07242  -0.12398  -0.03786  -0.05090  -0.05180
# via plm()
fatal_tefe_mod <- plm(fatal_rate ~ beertax, 
                      data = Fatalities,
                      index = c("state", "year"), 
                      model = "within", 
                      effect = "twoways")

coeftest(fatal_tefe_mod, vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##         Estimate Std. Error t value Pr(>|t|)  
## beertax -0.63998    0.35015 -1.8277  0.06865 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# check the class of 'state' and 'year'
class(Fatalities$state)
## [1] "factor"
class(Fatalities$year)
## [1] "factor"

Ở đây ta dùng lm() vì đây là mô hình hồi quy mở rộng với các biến giả. Tuy nhiên ta cũng có thể dùng hàm plm() và đặt argument effect = 'twoways' cho việc cố định hai chiều thời gian và không gian. Lưu ý ta cần đảm bảo stateyear có dạng factor. Hàm lm() tự động chuyển các biến factor thành dummies. Ở đây ta không ước lượng intercept nên trong công thức lm() ta đặt -1 ở cuối.

Ta thấy so với cố định không gian, khi cố định thêm chiều thời gian thì hệ số ước lượng không khác nhau nhiều, nhưng sai số to hơn và may mắn là vẫn đạt ý nghĩa thống kê thường thấy. Do đó ta có thể khẳng định, quan hệ giữa hai biến khảo sát không bị tác động bởi omitted variable bias từ những nhân tố cố định theo thời gian.

12.5 Giả định hồi quy tác động cố định

Đối với hồi quy cố định không gian, nhằm OLS cho ra các ước lượng không thiên lệch và tiệm cận Normal khi quy mô mẫu lớn, thì cần đảm bảo các giả định giống như với hồi quy đa biến. Cụ thể:

  • Sai số \(u_{it}\) có kỳ vọng có điều kiện bằng \(0\).
  • Các biến giải thích và sai số tuân theo \(i.i.d\).
  • \((X_{it},u_{it})\)moment thứ tư xác định và khác \(0\).
  • Không có đa cộng tuyến hoàn hảo.

Giả định đầu tiên đảm bảo rằng không có omitted variable bias. Giả định thứ hai chỉ yêu cầu các biến phải \(i.i.d\) qua ( across) các chủ thể, chứ không yêu cầu không tương quan trong ( within) chủ thể. Nói cách khác, \(X_{it}\) cho phép tự tương quan ( autocorrelated) trong chủ thể. Đây là một đặc trưng thường thấy trong dữ liệu chuỗi thời gian. Giả định thứ hai được đảm bảo nếu như chủ thể được chọn từ lấy mẫu ngẫu nhiên đơn giản. Giả định thứ ba và thứ tư đồng nhất về ý nghĩa như hồi quy đa biến.

12.5.1 Sai số chuẩn

Tương tự như phương sai sai số thay đổi, tự tương quan cũng gây nhiều phiền phức cho mô hình hồi quy tác động cố định.

Khi có cả phương sai sai số thay đổi và tự tương quan, sai số chuẩn HAC ( heteroskedasticity and autocorrelation-consistent (HAC) standard errors) được sử dụng, chẳng hạn như clustered standard errors CSE. Nó cho phép sai số do phương sai thay đổi và tự tương quan trong cùng một chủ thể, chứ không cho phép trong trường hợp các chủ thể. Trong R ta mượn hàm vcovHC() trong package sandwich để đo lường CSE.

# check class of the model object
class(fatal_tefe_lm_mod)
## [1] "lm"
# obtain a summary based on heteroskedasticity-robust standard errors 
# (no adjustment for heteroskedasticity only)
coeftest(fatal_tefe_lm_mod, vcov = vcovHC, type = "HC1")[1, ]
##   Estimate Std. Error    t value   Pr(>|t|) 
## -0.6399800  0.2547149 -2.5125346  0.0125470
# check class of the (plm) model object
class(fatal_tefe_mod)
## [1] "plm"        "panelmodel"
# obtain a summary based on clusterd standard errors 
# (adjustment for autocorrelation + heteroskedasticity)
coeftest(fatal_tefe_mod, vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##         Estimate Std. Error t value Pr(>|t|)  
## beertax -0.63998    0.35015 -1.8277  0.06865 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta thấy kết quả so sánh khá khác biệt. Không có tự tương quan ta đạt SE là \(0.25\) với hệ số \(\hat{\beta}_1\) có ý nghĩa mức \(5\%\). Tuy nhiên sử dụng CSE ta thấy kết quả SE là \(0.35\) dẫn đến hệ số này không còn ý nghĩa ở \(5\%\) nữa.

12.6 Luật giao thông và Tai nạn giao thông

Có hai nguồn chủ yếu gây ra OVB đối với tất cả mô hình liên quan đến tai nạn giao thông và thuế bia rượu: điều kiện kinh tế và luật giao thông. Trong bộ dữ liệu của chúng ta, Fatalities có chứa dữ liệu về tuổi uống rượu hợp pháp của từng bang drinkage, mức phạt jail, service và các chỉ số kinh tế như tỷ lệ thất nghiệp unemp và thu nhập bình quân đầu người income. Ta sẽ sử dụng các biến sau để mở rộng phân tích.

  • unemp: tỷ lệ thất nghiệp từng bang.
  • log(income): logarithm thu nhập bình quân so với năm 1988.
  • miles: trung bình dặm mỗi tài xế.
  • drinkage: mức tuổi tối thiểu được phép uống.
  • drinkagc: phiên bản rời rạc của drinkage.
  • punish: biến dummy sẽ là yes nếu bị phạt tù hoặc lao động công ích.
# discretize the minimum legal drinking age
Fatalities$drinkagec <- cut(Fatalities$drinkage,
                            breaks = 18:22, 
                            include.lowest = TRUE, 
                            right = FALSE)

# set minimum drinking age [21, 22] to be the baseline level
Fatalities$drinkagec <- relevel(Fatalities$drinkagec, "[21,22]")

# mandadory jail or community service?
Fatalities$punish <- with(Fatalities, factor(jail == "yes" | service == "yes", 
                                             labels = c("no", "yes")))

# the set of observations on all variables for 1982 and 1988
Fatalities_1982_1988 <- Fatalities[with(Fatalities, year == 1982 | year == 1988), ]

Sau đó, ta ước lượng \(7\) mô hình sử dụng plm().

# estimate all seven models
fatalities_mod1 <- lm(fatal_rate ~ beertax, data = Fatalities)

fatalities_mod2 <- plm(fatal_rate ~ beertax + state, data = Fatalities)

fatalities_mod3 <- plm(fatal_rate ~ beertax + state + year,
                       index = c("state","year"),
                       model = "within",
                       effect = "twoways", 
                       data = Fatalities)

fatalities_mod4 <- plm(fatal_rate ~ beertax + state + year + drinkagec 
                       + punish + miles + unemp + log(income), 
                       index = c("state", "year"),
                       model = "within",
                       effect = "twoways",
                       data = Fatalities)

fatalities_mod5 <- plm(fatal_rate ~ beertax + state + year + drinkagec 
                       + punish + miles,
                       index = c("state", "year"),
                       model = "within",
                       effect = "twoways",
                       data = Fatalities)

fatalities_mod6 <- plm(fatal_rate ~ beertax + year + drinkage 
                       + punish + miles + unemp + log(income), 
                       index = c("state", "year"),
                       model = "within",
                       effect = "twoways",
                       data = Fatalities)

fatalities_mod7 <- plm(fatal_rate ~ beertax + state + year + drinkagec 
                       + punish + miles + unemp + log(income), 
                       index = c("state", "year"),
                       model = "within",
                       effect = "twoways",
                       data = Fatalities_1982_1988)

Kết quả như sau:

Linear Panel Regression Models of Traffic Fatalities due to Drunk Driving
fatal_rate
OLS panel
linear
(1) (2) (3) (4) (5) (6) (7)
beertax 0.365*** -0.656** -0.640* -0.445 -0.690** -0.456 -0.926***
(0.053) (0.289) (0.350) (0.291) (0.345) (0.301) (0.337)
drinkagec[18,19) 0.028 -0.010 0.037
(0.068) (0.081) (0.101)
drinkagec[19,20) -0.018 -0.076 -0.065
(0.049) (0.066) (0.097)
drinkagec[20,21) 0.032 -0.100* -0.113
(0.050) (0.055) (0.123)
drinkage -0.002
(0.021)
punishyes 0.038 0.085 0.039 0.089
(0.101) (0.109) (0.101) (0.161)
miles 0.000 0.000* 0.000 0.000***
(0.000) (0.000) (0.000) (0.000)
unemp -0.063*** -0.063*** -0.091***
(0.013) (0.013) (0.021)
log(income) 1.816*** 1.786*** 0.996
(0.624) (0.631) (0.666)
Constant 1.853***
(0.047)
Observations 336 336 336 335 335 335 95
R2 0.093 0.041 0.036 0.360 0.066 0.357 0.659
Adjusted R2 0.091 -0.120 -0.149 0.217 -0.134 0.219 0.157
Residual Std. Error 0.544 (df = 334)
F Statistic 34.394*** (df = 1; 334) 12.190*** (df = 1; 287) 10.513*** (df = 1; 281) 19.194*** (df = 8; 273) 3.252*** (df = 6; 275) 25.423*** (df = 6; 275) 9.194*** (df = 8; 38)
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.

Ta đã khảo sát kết quả cột (2) và (3), cột (1) thể hiện ước lượng OLS thuần tuý cho hồi quy tỷ lệ tử vong và thuế bia không có bất kỳ tác động cố định nào. Ta có thể thấy rõ, ước lượng dương của hệ số thuế bia có khả năng thiên lệch hướng lên. Hệ số \(\bar{R}^2\) cũng không cao lắm. Sang mô hình (2) và (3) ta thấy dấu của hệ số hồi quy thay đổi.

Mô hình (4) đến (7) là mô hình mở rộng thêm đã đề cập ở trên. Ở mô hình (4) ta quan sát nhiều kết quả thú vị:

  • việc thêm vào các biến không nhất thiết làm giảm ảnh hưởng của hệ số thuế bia. Hệ số này không có ý nghĩa thống kê cho thấy khả năng cao là ước lượng không chính xác.
  • mức tuổi cho phép không tác động đến tỷ lệ tử vong: không có biến giả nào có ý nghĩa thống kê cả. Để xác nhận trường hợp này, ta xem xét kiểm định F.
# test if legal drinking age has no explanatory power
linearHypothesis(fatalities_mod4,
                 test = "F",
                 c("drinkagec[18,19)=0", "drinkagec[19,20)=0", "drinkagec[20,21)"), 
                 vcov. = vcovHC, type = "HC1")
## Linear hypothesis test
## 
## Hypothesis:
## drinkagec[18,19) = 0
## drinkagec[19,20) = 0
## drinkagec[20,21) = 0
## 
## Model 1: restricted model
## Model 2: fatal_rate ~ beertax + state + year + drinkagec + punish + miles + 
##     unemp + log(income)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F Pr(>F)
## 1    276                 
## 2    273  3 0.3782 0.7688
  • việc đề ra các mức phạt có lẽ chưa đủ răn đe: hệ số không có ý nghĩa thống kê.
  • các biến kinh tế giải thích rất tốt tỷ lệ tai nạn: các hệ số đều có ý nghĩa thống kê cao.

Kết quả mô hình (5) loại bỏ các biến kinh tế cho thấy các chỉ báo kinh tế nên được thêm vào, vì ảnh hưởng đến trị tuyệt đối của hệ số \(\hat{\beta}_1\). Mô hình (6) cho thấy tuổi luật pháp có rất ít năng lực giải thích. Mô hình (7) giảm lượng quan sát bằng cách thu hẹp khung thời gian cho thấy làm tăng SE và có sự chênh lệch nhiều về ước lượng hệ số.

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

12.7.1 Bài 1: Bộ dữ liệu Guns

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCB0aGUgcGFja2FnZSBhbmQgbG9hZCB0aGUgZGF0YXNldFxubGlicmFyeShBRVIpXG5cblxuIyBvYnRhaW4gYW4gb3ZlcnZpZXcgb3ZlciB0aGUgZGF0YXNldFxuXG5cbiMgdmVyaWZ5IHRoYXQgR3VucyBpcyBhIGJhbGFuY2VkIHBhbmVsXG55ZWFycyAgPC0gXG5zdGF0ZXMgPC0iLCJzb2x1dGlvbiI6IiMgYXR0YWNoIHRoZSBgQUVSYCBwYWNrYWdlIGFuZCBsb2FkIHRoZSBgR3Vuc2AgZGF0YXNldFxubGlicmFyeShBRVIpICBcbmRhdGEoXCJHdW5zXCIpICBcbiAgXG4jIG9idGFpbiBhbiBvdmVydmlldyBvdmVyIHRoZSBkYXRhc2V0XG5zdW1tYXJ5KEd1bnMpXG5cbiMgdmVyaWZ5IHRoYXQgYEd1bnNgIGlzIGEgYmFsYW5jZWQgcGFuZWxcbnllYXJzICA8LSBsZW5ndGgobGV2ZWxzKEd1bnMkeWVhcikpXG5zdGF0ZXMgPC0gbGVuZ3RoKGxldmVscyhHdW5zJHN0YXRlKSlcbnllYXJzKnN0YXRlcyA9PSBucm93KEd1bnMpIn0=

12.7.2 Bài 2: Thắt chặt hay thả lỏng 1?

Có nhiều cuộc tranh luận về vấn đề cho phép hay không cho phép việc mang súng ảnh hưởng đến tỷ lệ tội phạm. Đạo luật CCW nói rằng hạn chế sử dụng súng sẽ ngăn cản tội phạm trong khi những ý kiến trái ngược cho rằng việc cho phép sử dụng súng sẽ làm quá trình truy bắt tội phạm dễ dàng hơn.

Xem xét mô hình sau:

\[\log(\hat{violent}_i) = 6.135 - 0.443 \times law_i\]

trong đó violent là tỷ lệ tội phạm trên 100000 cư dân, law là biến dummy sẽ là yes nếu thực hiện luật CCW.

Thực hiện yêu cầu dưới đây.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGVzdGltYXRlIGEgbW9kZWwgd2l0aCBzdGF0ZSBmaXhlZCBlZmZlY3RzIHVzaW5nIHBsbSgpXG5tb2RlbF9zZSA8LSBcblxuIyBwcmludCBhIHN1bW1hcnkgdXNpbmcgcm9idXN0IHN0YW5kYXJkIGVycm9yc1xuXG5cbiMgdGVzdCB3aGV0aGVyIHRoZSBzdGF0ZSBmaXhlZCBlZmZlY3RzIGFyZSBqb2ludGx5IHNpZ25pZmljYW50IGZyb20gemVybyIsInNvbHV0aW9uIjoiIyBlc3RpbWF0ZSBhIG1vZGVsIHdpdGggc3RhdGUgZml4ZWQgZWZmZWN0cyB1c2luZyBwbG0oKVxubW9kZWxfc2UgPC0gcGxtKGxvZyh2aW9sZW50KSB+IGxhdywgZGF0YSA9IEd1bnMsIGluZGV4ID0gYyhcInN0YXRlXCIsIFwieWVhclwiKSwgbW9kZWwgPSBcIndpdGhpblwiKVxuXG4jIHByaW50IGEgc3VtbWFyeSB1c2luZyByb2J1c3Qgc3RhbmRhcmQgZXJyb3JzXG5jb2VmdGVzdChtb2RlbF9zZSwgdmNvdi4gPSB2Y292SEMsIHR5cGUgPSBcIkhDMVwiKVxuXG4jIHRlc3Qgd2hldGhlciB0aGUgc3RhdGUgZml4ZWQgZWZmZWN0cyBhcmUgam9pbnRseSBzaWduaWZpY2FudCBmcm9tIHplcm9cbnBGdGVzdChtb2RlbF9zZSwgbW9kZWwpIn0=

12.7.3 Bài 3: Thắt chặt hay thả lỏng 2?

Kết quả từ Bài 2 cho thấy ta nên thêm vào mô hình tác động thời gian:

\[\log(violent_i) = \beta_1 \times law_i + \alpha_i + \lambda_t + u_i\]

Thực hiện yêu cầu dưới đây.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGVzdGltYXRlIGEgbW9kZWwgd2l0aCBzdGF0ZSBhbmQgdGltZSBmaXhlZCBlZmZlY3RzIHVzaW5nIHBsbSgpXG5tb2RlbF9zZXRlIDwtIFxuXG4jIHByaW50IGEgc3VtbWFyeSB1c2luZyByb2J1c3Qgc3RhbmRhcmQgZXJyb3JzXG5cblxuIyB0ZXN0IHdoZXRoZXIgc3RhdGUgYW5kIHRpbWUgZml4ZWQgZWZmZWN0cyBhcmUgam9pbnRseSBzaWduaWZpY2FudCBmcm9tIHplcm8iLCJzb2x1dGlvbiI6IiMgZXN0aW1hdGUgYSBtb2RlbCB3aXRoIHN0YXRlIGFuZCB0aW1lIGZpeGVkIGVmZmVjdHMgdXNpbmcgcGxtKClcbm1vZGVsX3NldGUgPC0gcGxtKGxvZyh2aW9sZW50KSB+IGxhdywgZGF0YSA9IEd1bnMsIGluZGV4ID0gYyhcInN0YXRlXCIsIFwieWVhclwiKSwgbW9kZWwgPSBcIndpdGhpblwiLCBlZmZlY3QgPSBcInR3b3dheXNcIilcblxuIyBwcmludCBhIHN1bW1hcnkgdXNpbmcgcm9idXN0IHN0YW5kYXJkIGVycm9yc1xuY29lZnRlc3QobW9kZWxfc2V0ZSwgdmNvdi4gPSB2Y292SEMsIHR5cGUgPSBcIkhDMVwiKVxuXG4jIHRlc3Qgd2hldGhlciBzdGF0ZSBhbmQgdGltZSBmaXhlZCBlZmZlY3RzIGFyZSBqb2ludGx5IHNpZ25pZmljYW50IGZyb20gemVyb1xucEZ0ZXN0KG1vZGVsX3NldGUsIG1vZGVsKSJ9

12.7.4 Bài 4: Thắt chặt hay thả lỏng 3?

Có thể có OVB đối với mô hình ở Bài 3, do đó ta xem xét thêm mô hình dưới đây.

\[\begin{aligned}\log(violent_i) = &\beta_1 \times law_i + \beta_2\times density_i + \beta_3 \times income_i + \beta_4 \times population_i \\&+ \beta_5 \times afarm_i + \beta_6 \times cauc_i + beta_7 \times male_i + \alpha_i + \lambda_t + u_i \end{aligned}\]

Để biết thông tin các biến, dùng ?Guns.

Thực hiện yêu cầu dưới đây.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGVzdGltYXRlIHRoZSBleHRlbmRlZCBtb2RlbFxubW9kZWxfc2V0ZV9leHQgPC0gXG5cbiMgcHJpbnQgYSBzdW1tYXJ5IHVzaW5nIHJvYnVzdCBzdGFuZGFyZCBlcnJvcnMiLCJzb2x1dGlvbiI6IiMgZXN0aW1hdGUgdGhlIGV4dGVuZGVkIG1vZGVsXG5tb2RlbF9zZXRlX2V4dCA8LSBwbG0obG9nKHZpb2xlbnQpIH4gbGF3ICsgcHJpc29uZXJzICsgZGVuc2l0eSArIGluY29tZSArIHBvcHVsYXRpb24gKyBhZmFtICsgY2F1YyArIG1hbGUsIGRhdGEgPSBHdW5zLCBpbmRleCA9IGMoXCJzdGF0ZVwiLCBcInllYXJcIiksIG1vZGVsID0gXCJ3aXRoaW5cIiwgZWZmZWN0ID0gXCJ0d293YXlzXCIpXG5cbiMgcHJpbnQgYSBzdW1tYXJ5IHVzaW5nIHJvYnVzdCBzdGFuZGFyZCBlcnJvcnNcbmNvZWZ0ZXN0KG1vZGVsX3NldGVfZXh0LCB2Y292LiA9IHZjb3ZIQywgdHlwZSA9IFwiSEMxXCIpIn0=

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

  2. Croissant, Y., Millo, G., & Tappe, K. (2017). plm: Linear Models for Panel Data.