Chương 6 Hồi quy tuyến tính một biến

Những packages sau đây sẽ được dùng trong chương này, nhằm rút ngắn các hàm lập trình phức tạp và tận dụng bộ dữ liệu:

  • AER - package đi kèm cùng cuốn sách Applied Econometrics with R (Kleiber và Zeileis, 20082).
  • MASS - package chứa nhiều hàm ứng dụng.

6.1 Hồi quy tuyến tính đơn giản

Để bắt đầu với một ví dụ đơn giản, ta xem xét tổ hợp điểm trung bình và tỷ lệ học sinh - giáo viên trung bình trong một số trường học giả tưởng.

STR <- c(15, 17, 19, 20, 22, 23.5, 25)
TestScore <- c(680, 640, 670, 660, 630, 660, 635)

Trong mô hình hồi quy đơn giản, ta mô hình mối quan hệ giữa hai biến này bằng một đường thẳng:

\[Y = a + bX\]

Giả sử hàm số cho ví dụ chúng ta là:

\[\text{TestScore} = 713 -3 \times \text{STR}\]

Trước tiên ta cần hình dung bộ dữ liệu của chúng ta được đồ thị hoá như thế nào. Trong R ta thực hiện như sau:

plot(TestScore~STR)
abline(a=713, b=-3)

Rõ ràng ta thấy đường thẳng không hề đi qua điểm nào mặc dù chúng ta xác nhận rằng nó đại diện cho quan hệ tổng thể (quan hệ hệ thống). Nguyên nhân chính là sự ngẫu nhiên. Trong thực tế ta sẽ không bao giờ tìm thấy một mối quan hệ tuyệt đối nào mà không có sự hiện diện của tính ngẫu nhiên. Đồng thời, hầu hết các trường hợp ta sẽ thấy ngoài tính ngẫu nhiên còn có một số tác động cho thấy không có quan hệ nhị biến giữa hai biến.

Khi tính đến sự khác biệt giữa dữ liệu và quan hệ hệ thống, ta thêm vào mô hình một sai số (error term) \(u\) thể hiện tính ngẫu nhiên. Về bản chất mà nói, \(u\) ngoài tính ngẫu nhiên còn thể hiện sự chênh lệch do sai số phép đo, hoặc do bỏ biến có liên quan đến biến được giải thích \(Y\).

Trong ví dụ có thể có những nhân tố khác tác động đến TestScore, chẳng hạn như chất lượng giáo viên và trình độ học sinh. Cũng có khả năng đề thi dễ hơn vào một số ngày. Để tổng quát hoá, ta dùng mô hình sau đây:

\[\text{TestScore} = 713 -3 \times \text{STR} + \text{Other factors}\]

Như vậy, mô hình hồi quy tuyến tính cơ bản nhất được hình thành theo dạng sau:

\[Y_i = \beta_0 + \beta_1 X_i + u_i\]

trong đó: \(Y\) là biến phụ thuộc, hay biến được giải thích, \(X\) là biến độc lập hay biến giải thích, \(\beta_0\) gọi là hệ số chặn (trong môn học sẽ gọi tên quốc tế là intercept) và \(\beta_1\) là độ nghiêng tuyến tính (trong môn học sẽ gọi tên quốc tế là slope).

6.2 Ước lượng hệ số hồi quy

Trong thực tế, chỉ khi ta khảo sát một tổng thể, ta mới biết chắc chắn giá trị các hệ số \(\beta\), nhưng điều này là bất khả thi. Bài toán trở thành, dựa vào một mẫu quan sát ta ước lượng các tham số hồi quy chưa biết. Để cụ thể hơn, ta dùng dữ liệu cacs trường học tại Canada làm ví dụ.

library(AER)
data(CASchools)
head(CASchools)
##   district                          school  county grades students
## 1    75119              Sunol Glen Unified Alameda  KK-08      195
## 2    61499            Manzanita Elementary   Butte  KK-08      240
## 3    61549     Thermalito Union Elementary   Butte  KK-08     1550
## 4    61457 Golden Feather Union Elementary   Butte  KK-08      243
## 5    61523        Palermo Union Elementary   Butte  KK-08     1335
## 6    62042         Burrel Union Elementary  Fresno  KK-08      137
##   teachers calworks   lunch computer expenditure    income   english  read
## 1    10.90   0.5102  2.0408       67    6384.911 22.690001  0.000000 691.6
## 2    11.15  15.4167 47.9167      101    5099.381  9.824000  4.583333 660.5
## 3    82.90  55.0323 76.3226      169    5501.955  8.978000 30.000002 636.3
## 4    14.00  36.4754 77.0492       85    7101.831  8.978000  0.000000 651.9
## 5    71.50  33.1086 78.4270      171    5235.988  9.080333 13.857677 641.8
## 6     6.40  12.3188 86.9565       25    5580.147 10.415000 12.408759 605.7
##    math
## 1 690.0
## 2 661.9
## 3 650.9
## 4 643.5
## 5 639.9
## 6 605.4

Ta có thể thấy dữ liệu bao gồm khá nhiều biến, và hầu hết chúng đều là ở dạng numeric. Ở đây cả hai biến STRTestScore đều không có sẵn, bắt buộc ta phải tự tính từ dữ liệu.

CASchools$STR <- CASchools$students/CASchools$teachers 
CASchools$Score <- (CASchools$read + CASchools$math)/2

Bước tiếp theo, ta có thể nhìn lướt qua thống kê mô tả hai biến này:

summary(data.frame(CASchools$STR,CASchools$Score))
##  CASchools.STR   CASchools.Score
##  Min.   :14.00   Min.   :605.5  
##  1st Qu.:18.58   1st Qu.:640.0  
##  Median :19.72   Median :654.5  
##  Mean   :19.64   Mean   :654.2  
##  3rd Qu.:20.87   3rd Qu.:666.7  
##  Max.   :25.80   Max.   :706.8

Như trên, ta tiếp tục nhìn đồ thị điểm rải của dữ liệu:

plot(Score~STR, data = CASchools, main = "Đồ thị điểm rải", xlab = "STR", ylab = "Score")

Từ kiến thức các chương trước, ta có thể thấy hai biến quan tâm có tương quan âm (nhưng khá yếu). Có nghĩa là điểm số sẽ thấp hơn nếu như quy mô lớp học lớn hơn. Để tính toán chính xác hệ số tương quan, ta dùng hàm cor trong R.

cor(CASchools$STR, CASchools$Score)
## [1] -0.2263627

Bây giờ ta sẽ đi ước lượng mô hình hồi quy bằng phương pháp OLS (Bình phương nhỏ nhất - Ordinary Least Squares). Phương pháp này sẽ chọn các hệ số hồi quy sao cho đường hồi quy gần với các giá trị quan sát nhất có thể. Việc đo lường “gần” hay “xa” sẽ dựa vào tổng bình phương chênh lệch giữa quan sát và giá trị đo lường trên đường hồi quy. Gọi \(b_0\)\(b_1\) là một estimator nào đó của các hệ số \(\beta\) tương ứng thì tổng bình phương chênh lệch sẽ được tính theo công thức:

\[\sum_{i=1}^n (Y_i - b_0 - b_1 X_i)^2\]

Phương pháp OLS sẽ tìm kiếm các hệ số \(b_0\)\(b_1\) sao cho công thức trên đạt giá trị nhỏ nhất. Dựa vào ý tưởng như vậy, ta tìm ra công thức tổng quát như sau:

\[\begin{aligned} &\hat{\beta_1} = \frac{\sum (X_i - \bar{X})(Y_i - \bar{Y})}{\sum (X_i -\bar{X})^2}, \\ &\hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X} \end{aligned}\]

Trong R ta có thể tính toán các hệ số như sau:

attach(CASchools) # allows to use the variables contained in CASchools directly
## The following object is masked _by_ .GlobalEnv:
## 
##     STR
## The following objects are masked from CASchools (pos = 3):
## 
##     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 = 6):
## 
##     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 = 9):
## 
##     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 = 17):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, STR, students,
##     teachers
## The following objects are masked from CASchools (pos = 22):
## 
##     calworks, computer, county, district, english, expenditure,
##     grades, income, lunch, math, read, school, Score, STR,
##     students, teachers
beta_1 <- sum((STR - mean(STR)) * (Score - mean(Score))) / sum((STR - mean(STR))^2)
beta_0 <- mean(Score) - beta_1 * mean(STR)
data.frame(beta_1, beta_0)
##      beta_1   beta_0
## 1 -5.548841 766.3224

Ngoài ra, để tối ưu hoá năng lực của R ta có thể sử dụng trực tiếp hàm dựng sẵn lm như sau và đạt được kết quả tương tự:

linear_model <- lm(Score~STR, data = CASchools)
linear_model
## 
## Call:
## lm(formula = Score ~ STR, data = CASchools)
## 
## Coefficients:
## (Intercept)          STR  
##      698.93        -2.28

Hình hoạ hoá mô hình như sau:

plot(Score~STR, data = CASchools, xlim = c(10,30), ylim = c(600,700))
abline(linear_model)

6.3 Đo lường năng lực mô hình

Sau khi tìm ra được mô hình, một câu hỏi đặt ra là mô hình mô tả dữ liệu tốt như thế nào. Về mặt hình học, ta có thể xem xét bằng mắt và đánh giá mức độ phù hợp của mô hình. Nhưng để lượng hoá, ta có thể sử dụng hệ số xác định \(R^2\)sai số hồi quy chuẩn \(SER\).

6.3.1 Hệ số xác định

Ta định nghĩa:

  • \(ESS\) là tổng bình phương chênh lệch giữa giá trị tuyến tính và giá trị trung bình của \(Y\).
  • \(TSS\) là tổng bình phương chênh lệch giữa giá trị quan sát và giá trị trung bình của \(Y\).

Khi đó \(R^2\) được định nghĩa là tỷ số giữa \(ESS\)\(TSS\), thể hiện phần trăm phương sai của \(Y\) được giải thích bởi \(X\).

\[R^2 = \frac{ESS}{TSS}\]

Nếu định nghĩa \(SSR = TSS - ESS = \sum \hat{u}_i^2\) thì:

\[R^2 = 1-\frac{SSR}{TSS}\]

\(R^2\) nằm giữa \(0\)\(1\). Nếu mô hình hoàn toàn trùng khớp dữ liệu, rõ ràng \(R^2 = 1\), ngược lại, \(R^2\) sẽ tiến gần về \(0\).

6.3.2 Sai số hồi quy chuẩn

\(SER\) là một estimator cho độ lệch chuẩn của \(\hat{u}_i\). Theo đó:

\[SER = s_{\hat{u}} = \frac{SSR}{n-2}\]

Trong R các tham số kể trên đều có thể tìm thấy nhờ hàm summary.

summary(linear_model)
## 
## Call:
## lm(formula = Score ~ STR, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.727 -14.251   0.483  12.822  48.540 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 698.9329     9.4675  73.825  < 2e-16 ***
## STR          -2.2798     0.4798  -4.751 2.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared:  0.05124,    Adjusted R-squared:  0.04897 
## F-statistic: 22.58 on 1 and 418 DF,  p-value: 2.783e-06

Theo đó, \(R^2\) bằng \(0.051\) tức \(5.1\%\) phương sai \(Y\) được giải thích bởi \(X\), khá thấp. \(SER = 18.58\) nghĩa là về mặt trung bình, chênh lệch giữa điểm thực tế và điểm hồi quy là 18.58 điểm (\(SER\) cùng đơn vị với biến \(Y\)).

6.4 Giả định hồi quy

Để OLS hoạt động tốt, ta cần đảm bảo dữ liệu thoả mãn các giả định trọng yếu sau:

  • Sai số \(u\): \(E(u_i|X_i)=0\)
  • \(X,Y \sim i.i.d\)
  • \(X\)\(Y\) cần có moment bậc \(4\) khác \(0\), nghĩa là các giá trị outliers không thường xuyên.

6.5 Phân phối lấy mẫu cho OLS

6.5.1 Phân phối các hệ số hồi quy

Bởi vì các hệ số \(\hat{\beta}\) đều được tính toán từ mẫu, bản thân các estimators này là biến ngẫu nhiên với một phân phối xác suất. Miễn là các giả định hồi quy được thoả mãn, thì estimators được định nghĩa theo OLS chính là ước lượng không thiên lệch của giá trị \(\beta\) đúng.

Nếu quy mô mẫu đủ lớn, nhờ vào CLT ta biết rằng phân phối xấp xỉ các estimator sẽ tiệm cận phân phối Normal đa biến. Khi đó các phân phối biên sẽ được định nghĩa như sau.

\[\begin{aligned} &\hat{\beta_1} \sim N(\beta_1, \sigma^2_{\hat{\beta}_1}): \sigma^2_{\hat{\beta}_1} = \frac{1}{n}\frac{\text{Var}[(X_i - \mu_X)u_i]}{[\text{Var}(X_i)]^2} \\ &\hat{\beta_0} \sim N(\beta_0, \sigma^2_{\hat{\beta_0}}):\sigma^2_{\hat{\beta_0}} = \frac{1}{n}\frac{\text{Var}(H_iu_i)}{[E(H_i)^2]^2} \end{aligned}\]

trong đó: \(H_i= 1- \left[\frac{\mu_X}{E(X_i)^2} \right]X_i\).

Ta có thể kiểm tra điều này thông qua R. Ta xây dựng một tổng thể \(100000\) quan sát, các biến \(X\) được lấy ngẫu nhiên từ phân phối Uniform \(U(0,20)\), các error term được lấy từ phân phối Normal \(N(0,100)\) và mô hình giả định như sau:

\[Y_i = -2+3.5 X_i + u_i\]

N <- 100000
X <- runif(N,0,20)
u <- rnorm(N,mean = 0, sd = 10)
Y <- -2 + 3.5*X+u
population <- data.frame(X,Y)

Trước tiên, ta cần tính phương sai đúng của các hệ số hồi quy đối với mẫu có quy mô \(n=100\).

n <- 100
H_i <- 1-mean(X)/mean(X^2) * X
var_b0 <- var(H_i *u)/ (n*mean(H_i^2)^2)
var_b1 <- var((X-mean(X))*u) / (100*var(X)^2)
data.frame(var_b0,var_b1)
##    var_b0     var_b1
## 1 3.95442 0.02953373

Bây giờ, ta tiến hành lấy mẫu ngẫu nhiên nhằm xác định phân phối lấy mẫu.

reps <- 10000
fit <- matrix(ncol = 2, nrow = reps) # Lưu trữ kết quả trong ma trận này
for (i in 1:reps){
  sample <- population[sample(1:N,n),]
  fit[i,] <- lm(Y~X, data = sample)$coefficients
}
data.frame(var(fit[,1]),var(fit[,2]))
##   var.fit...1.. var.fit...2..
## 1      4.051558    0.03039434

Hình hoạ hoá phân phối lấy mẫu:

par(mfrow = c(1, 2))
hist(fit[, 1],cex.main = 1,main = bquote(The ~ Distribution  ~ of ~ 10000 ~ beta[0] ~ Estimates), xlab=bquote(hat(beta)[0]), freq = F)
curve(dnorm(x, -2, sqrt(var_b0)),add = T, col = "darkred")
hist(fit[, 2],cex.main = 1,main = bquote(The ~ Distribution  ~ of ~ 10000 ~ beta[1] ~ Estimates), xlab = bquote(hat(beta)[1]), freq = F)
curve(dnorm(x, 3.5, sqrt(var_b1)), add = T, col = "darkred")

Ta thấy rõ, ước lượng phương sai rất gần với giá trị đúng. Đồ thị histogram cho thấy phân phối của estimator có thể xấp xỉ với phân phối Normal với quy cách đã cho.

6.5.2 Tính consistency của hệ số hồi quy

Các estimators đều consistent, nghĩa là chúng hội tụ về mặt xác suất về các giá trị đúng. Lý do là vì chúng xấp xỉ không thiên lệch và phương sai hội tụ về \(0\) khi \(n\) tăng.

Chẳng hạn, ta thí nghiệm với hệ số \(\beta_1\) như sau.

# set seed for reproducibility
set.seed(1)

# set repetitions and the vector of sample sizes
reps <- 1000
n <- c(100, 250, 1000, 3000)

# initialize the matrix of outcomes
fit <- matrix(ncol = 2, nrow = reps)

# divide the plot panel in a 2-by-2 array
par(mfrow = c(2, 2))

# loop sampling and plotting

# outer loop over n
for (j in 1:length(n)) {
  
  # inner loop: sampling and estimating of the coefficients
  for (i in 1:reps){
    
    sample <- population[sample(1:N, n[j]), ]
    fit[i, ] <- lm(Y ~ X, data = sample)$coefficients
    
  }
  
  # draw density estimates
  plot(density(fit[ ,2]), xlim=c(2.5, 4.5), 
       col = j, 
       main = paste("n=", n[j]), 
       xlab = bquote(hat(beta)[1]))
  
}

Ta có thể thấy rõ, khi \(n\) tăng, phân phối \(\hat{\beta}_1\) càng tập trung xung quanh trung bình của nó, nghĩa là phương sai giảm, đồng thời phân phối càng giống Normal hơn.

6.6 Bài tập

Bài 1: Quy mô lớp và Điểm thi

Một nhà nghiên cứu đo lường quy mô lớp (bằng tỷ lệ giảng viên-sinh viên) và điểm thi trung bình của \(10\) lớp học và có được kết quả như sau.

Quy mô lớp 23 19 30 22 23 29 35 36 33 25
Điểm thi 430 430 333 410 390 377 325 310 328 375
  1. Tạo một vector cs (class size - quy mô lớp học) và ts (test score - điểm thi) gồm những quan sát trên.
  2. Vẽ đồ thị điểm rải.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNyZWF0ZSBib3RoIHZlY3RvcnNcblxuIyBkcmF3IHRoZSBzY2F0dGVycGxvdCIsInNvbHV0aW9uIjoiIyBjcmVhdGUgYm90aCB2ZWN0b3JzXG5jcyA8LSBjKDIzLCAxOSwgMzAsIDIyLCAyMywgMjksIDM1LCAzNiwgMzMsIDI1KVxudHMgPC0gYyg0MzAsIDQzMCwgMzMzLCA0MTAsIDM5MCwgMzc3LCAzMjUsIDMxMCwgMzI4LCAzNzUpXG5cbiMgZHJhdyB0aGUgc2NhdHRlcnBsb3RcbnBsb3QoY3MsdHMpIn0=

Bài 2: Trung bình, phương sai, Hiệp phương sai, Hệ số tương quan

  1. Tính trung bình, phương sai, độ lệch chuẩn mẫu của ts.
  2. Tính hiệp phương sai, hệ số tương quan mẫu của tscs
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbXB1dGUgbWVhbiwgdmFyaWFuY2UgYW5kIHN0YW5kYXJkIGRldmlhdGlvbiBvZiB0ZXN0IHNjb3Jlc1xuXG5cbiMgY29tcHV0ZSB0aGUgY292YXJpYW5jZSBhbmQgdGhlIGNvcnJlbGF0aW9uIGNvZWZmaWNpZW50Iiwic29sdXRpb24iOiIjIGNvbXB1dGUgbWVhbiwgdmFyaWFuY2UgYW5kIHN0YW5kYXJkIGRldmlhdGlvbiBvZiB0ZXN0IHNjb3Jlc1xubWVhbih0cylcbnNkKHRzKVxudmFyKHRzKVxuXG4jIGNvbXB1dGUgdGhlIGNvdmFyaWFuY2UgYW5kIHRoZSBjb3JyZWxhdGlvbiBjb2VmZmljaWVudFxuY292KHRzLCBjcylcbmNvcih0cywgY3MpIn0=

Bài 3: Hồi quy đơn giản

Thực hiện hồi quy \(\text{Điểm thi}_i = \beta_0 + \beta_1 \text{Quy mô lớp} + u_i\) bằng hàm lm() trong gói AER, lưu trữ vào biến mod.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGF0dGFjaCB0aGUgcGFja2FnZSBBRVJcblxuIyBlc3RpbWF0ZSB0aGUgbW9kZWwiLCJzb2x1dGlvbiI6IiMgYXR0YWNoIHRoZSBwYWNrYWdlIEFFUlxubGlicmFyeShBRVIpXG5cbiMgZXN0aW1hdGUgdGhlIG1vZGVsXG5tb2QgPC0gbG0odHMgfiBjcykifQ==

Bài 4: Vẽ đường hồi quy

Vẽ đường hồi quy vào đồ thị điểm rải đã được cho trước.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGFkZCB0aGUgcmVncmVzc2lvbiBsaW5lIHRvIHRoZSBzY2F0dGVycGxvdFxucGxvdChjcywgdHMpIiwic29sdXRpb24iOiIjIGFkZCB0aGUgcmVncmVzc2lvbiBsaW5lIHRvIHRoZSBzY2F0dGVycGxvdFxucGxvdChjcywgdHMpXG5hYmxpbmUobW9kKSJ9

Bài 5: Đọc kết quả hồi quy

Xem xét kết quả hồi quy và hệ số xác định \(R^2\) bằng hàm summary.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIHN1bW1hcnkgbW9kZWwiLCJzb2x1dGlvbiI6IiMgc3VtbWFyeSBtb2RlbFxuc3VtbWFyeShtb2QpIn0=

Bài 6: \(TSS\)\(SSR\)

  1. Tính \(SSR\) và lưu vào ssr
  2. Tính \(TSS\) và lưu vào tss.

Chú ý đến bậc tự do.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbXB1dGUgdGhlIFNTUiBhbmQgc2F2ZSBpdCB0byBgc3NyYFxuXG5cbiMgY29tcHV0ZSB0aGUgVFNTIGFuZCBzYXZlIGl0IHRvIGB0c3NgIiwic29sdXRpb24iOiIjIGNvbXB1dGUgdGhlIFNTUiBhbmQgc2F2ZSBpdCB0byBgc3NyYFxuc3NyIDwtIHN1bShtb2QkcmVzaWR1YWxzXjIpXG5cbiMgY29tcHV0ZSB0aGUgVFNTIGFuZCBzYXZlIGl0IHRvIGB0c3NgXG50c3MgPC0gOSp2YXIodHMpXG4jIE5vdGU6IGB2YXIoKWAgY29tcHV0ZXMgdGhlIHVuYmlhc2VkIHNhbXBsZSB2YXJpYW5jZSEgPT4gY29ycmVjdCBmb3IgdGhpcyBieSBtdWx0aXBseWluZyB3aXRoIChuLTEpID0gOSJ9


  1. Kleiber, C., & Zeileis, A. (2008). Applied Econometrics with R. Springer.