Chương 9 Kiểm định giả thuyết và khoảng tin cậy hồi quy đa biến
Chương này mở rộng từ chương đối với hồi quy đa biến. Trong chương này, hai package sẽ được sử dụng:
9.1 Kiểm định và khoảng tin cậy cho một hệ số
Các bước kiểm định giả thuyết
\[\begin{aligned} &H_0: \beta_j = \beta_{j,0} \\ &H_A: \beta_j \ne \beta_{j,0} \end{aligned}\]
- Tính sai số chuẩn của \(\hat{\beta}_j\)
- Tính thống kê \(t\): \(t^{act} = \frac{\hat{\beta}_j - \beta_{j,0}}{SE(\hat{\beta}_j)}\)
- Tính p-value: \(p\text{ - value} = 2 \Phi (-|t^{act}|)\)
Như vậy, kiểm định đơn hệ số trong hồi quy đa biến thực hiện giống hệt như kiểm định trong hồi quy đơn giản.
data("CASchools")
CASchools$score <- (CASchools$read + CASchools$math)/2
CASchools$size <- CASchools$students/CASchools$teachers
model <- lm(score ~ size + english, data = CASchools)
coeftest(model, vcov. = vcovHC, type = "HC1")##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 686.032245 8.728225 78.5993 < 2e-16 ***
## size -1.101296 0.432847 -2.5443 0.01131 *
## english -0.649777 0.031032 -20.9391 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ta cũng có thể tự tính giá trị \(p\)-value đối với mô hình vừa thực hiện trong R như sau.
# compute two-sided p-value
2 * (1 - pt(abs(coeftest(model, vcov. = vcovHC, type = "HC1")[2, 3]),
df = model$df.residual))## [1] 0.01130921Khoảng tin cậy được xây dựng theo công thức:
\[\left[\hat{\beta}_j - 1.96 \times SE(\hat{\beta}_j), \hat{\beta}_j + 1.96 \times SE(\hat{\beta}_j) \right]\]
9.2 Ứng dụng vào Điểm thi và tỷ lệ \(SRT\)
Ta có thể xây khoảng tin cậy đồng thời cho các hệ số hồi quy bằng hàm confint() trong R.
model <- lm(score ~ size + english, data = CASchools)
confint(model)## 2.5 % 97.5 %
## (Intercept) 671.4640580 700.6004311
## size -1.8487969 -0.3537944
## english -0.7271113 -0.5724424Để có được khoảng tin cậy cho một mức tin cậy khác, chẳng hạn \(90\%\),
confint(model, level = 0.9)## 5 % 95 %
## (Intercept) 673.8145793 698.2499098
## size -1.7281904 -0.4744009
## english -0.7146336 -0.5849200Hàm confint() không xây dựng khoảng tin cậy vững ( robust confidence interval) dựa vào robust SE. Ta có thể điều chỉnh bằng cách:
# compute robust standard errors
rob_se <- diag(vcovHC(model, type = "HC1"))^0.5
# compute robust 95% confidence intervals
rbind("lower" = coef(model) - qnorm(0.975) * rob_se,
"upper" = coef(model) + qnorm(0.975) * rob_se)## (Intercept) size english
## lower 668.9252 -1.9496606 -0.7105980
## upper 703.1393 -0.2529307 -0.5889557# compute robust 90% confidence intervals
rbind("lower" = coef(model) - qnorm(0.95) * rob_se,
"upper" = coef(model) + qnorm(0.95) * rob_se)## (Intercept) size english
## lower 671.6756 -1.8132659 -0.7008195
## upper 700.3889 -0.3893254 -0.5987341Khi ta thêm biến expenditure vào mô hình, ta nhận được kết quả như sau.
# scale expenditure to thousands of dollars
CASchools$expenditure <- CASchools$expenditure/1000
# estimate the model
model <- lm(score ~ size + english + expenditure, data = CASchools)
coeftest(model, vcov. = vcovHC, type = "HC1")##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 649.577947 15.458344 42.0212 < 2e-16 ***
## size -0.286399 0.482073 -0.5941 0.55277
## english -0.656023 0.031784 -20.6398 < 2e-16 ***
## expenditure 3.867901 1.580722 2.4469 0.01482 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ta thấy rõ biến size từ có ý nghĩa thống kê trở thành không có ý nghĩa thống kê với \(p\)-value mới là \(0.55\). Như vậy việc thêm vào biến expenditure đã làm tăng sai số chuẩn của \(\hat{\beta}\), như vậy rất có thể có hiện tượng đa cộng tuyến phi hoàn hảo giữa expenditure và size. Ta có thể kiểm tra lại bằng cách xem xét hệ số tương quan giữa hai biến này.
# compute the sample correlation between 'size' and 'expenditure'
cor(CASchools$size, CASchools$expenditure)## [1] -0.6199822Tóm lại, ta có thể kết luận mô hình mới cung cấp bằng chứng cho thấy sự thay đổi quy mô lớp học không ảnh hưởng đến điểm thi khi giữ chi phí mỗi học sinh và phần trăm học sinh học tiếng Anh ở mức cố định.
9.3 Kiểm định gộp sử dụng thống kê \(F\)
Một câu hỏi đặt ra là: liệu giả thuyết hệ số hồi quy cho biến size và expenditure cùng bằng \(0\)?
Mặc dù việc kiểm định từng tham số cho mỗi biến là khả thi nhờ vào thống kê \(t\), cách tiếp cận như vậy không đáng tin cậy. Như vậy câu hỏi đặt ra sẽ chuyển thành: bác bỏ giả thuyết gộp \(H_0\) nếu hoặc \(t_1\) hoặc \(t_2\) đều lớn hơn \(1.96\) về mặt tuyệt đối.
Bởi vì câu hỏi liên quan đến \(t_1\) và \(t_2\), để trả lời được câu hỏi giả thuyết ta cần một phân phối lấy mẫu hợp của \(t_1\) và \(t_2\). Do \(\hat{\beta}_1\) và \(\hat{\beta}_2\) tuân theo phân phối Normal nhị biến, nên các thống kê \(t\) tương ứng cũng sẽ tuân theo phân phối Normal nhị biến, với trung bình là \(0\) và phương sai \(1\).
Nhắc lại, pdf của một phân phối hợp Normal có dạng như sau.
\[\begin{aligned} g_{X,Y}(x,y) = \frac{1}{2\pi \sigma_X \sigma_Y \sqrt{1-\rho^2}} \exp\left[ \frac{1}{-2(1-\rho^2)}\left( \left(\frac{x-\mu_X}{\sigma_X}\right)^2 \\ -2\rho\left(\frac{x-\mu_X}{\sigma_X}\right)\left(\frac{y-\mu_Y}{\sigma_Y}\right)+\left(\frac{y-\mu_Y}{\sigma_Y}\right)^2 \right) \right] \end{aligned}\]
Xét trường hợp đặc biệt nhất là tương quan giữa các thống kê \(t\) là bằng \(0\). Để không bác bỏ \(H_0\) ta cần \(|t_1| \le 1.96\) và \(|t_2| \le 1.96\). Bởi vì chúng độc lập nên xác suất \(P(|t_1| \le 1.96 \text{ và } |t_2| \le 1.96) = 0.95^2 = 0.9025\) dẫn đến xác suất bác bỏ \(H_0\) sẽ là \(1-0.9025 = 9.75\%\). Việc bác bỏ \(H_0\) xảy ra thường xuyên hơn nếu sử dụng cách này.
Xét trường hợp còn lại khi tương quan giữa các thống kê \(t\) khác \(0\). Vấn đề sẽ nghiêm trọng hơn vì còn phụ thuộc vào các xác suất có điều kiện của \(\rho\). Giả sử \(\rho = 0.6\) thì xác suất \(P(|t_1| \le 1.96 \text{ và } |t_2| \le 1.96) = 0.9125\). Khi đó vẫn dẫn đến kết luận bác bỏ \(H_0\) xảy ra thường xuyên hơn nếu sử dụng cách này. Trong R ta tính giá trị xác suất như sau.
InnerFunc <- function(x,y,rho=0.6){
1/(2*pi*sqrt(1-rho^2)) * exp(-1/(2*(1-rho^2)) * (x^2+y^2 - 2*rho*x*y))
}
InnerIntegral <- function(y){sapply(y, function(z) { integrate(InnerFunc, lower = -1.96, upper = 1.96,y=z)$value})}
integrate(InnerIntegral, lower = -1.96, upper = 1.96)$value## [1] 0.912458Như vậy ta cần một phương pháp kiểm định mới.
Trong trường hợp phương sai sai số đồng nhất, thống kê \(F\) có dạng:
\[\begin{equation} F= \frac{(SSR_{\text{restricted}}-SSR_{\text{unrestricted}})/q}{SSR_{\text{unrestricted}}/(n-k-1)} \end{equation}\]trong đó \(SSR_{\text{restricted}}\) là tổng bình phương sai số từ mô hình hồi quy giới hạn, \(q\) là số lượng giới hạn, \(SSR_{\text{unrestricted}}\) là tổng bình phương sai số từ mô hình đầy đủ.
Trong R ta có thể tính thống kê \(F\) dựa vào hàm linearHypothesis() của package car.
# estimate the multiple regression model
model <- lm(score ~ size + english + expenditure, data = CASchools)
# execute the function on the model object and provide both linear restrictions
# to be tested as strings
linearHypothesis(model, c("size=0", "expenditure=0"))## Linear hypothesis test
##
## Hypothesis:
## size = 0
## expenditure = 0
##
## Model 1: restricted model
## Model 2: score ~ size + english + expenditure
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 418 89000
## 2 416 85700 2 3300.3 8.0101 0.000386 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Kết quả cho thấy thống kê \(F\) bằng \(8.01\) và \(p\)-value tương ứng là \(0.0004\). Do đó, ta có thể bác bỏ \(H_0\).
Trong trường hợp phương sai sai số thay đổi, thống kê \(F\) được tính như sau.
# heteroskedasticity-robust F-test
linearHypothesis(model, c("size=0", "expenditure=0"), white.adjust = "hc1")## Linear hypothesis test
##
## Hypothesis:
## size = 0
## expenditure = 0
##
## Model 1: restricted model
## Model 2: score ~ size + english + expenditure
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df F Pr(>F)
## 1 418
## 2 416 2 5.4337 0.004682 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hàm summary() vẫn có thể lấy được thống kê \(F\), nhưng là phiên bản phương sai sai số đồng nhất.
# Access the overall F-statistic from the model's summary
summary(model)$fstatistic## value numdf dendf
## 107.4547 3.0000 416.0000# execute the function on the model object and provide the restrictions
# to be tested as a character vector
linearHypothesis(model, c("size=0", "english=0", "expenditure=0"))## Linear hypothesis test
##
## Hypothesis:
## size = 0
## english = 0
## expenditure = 0
##
## Model 1: restricted model
## Model 2: score ~ size + english + expenditure
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 419 152110
## 2 416 85700 3 66410 107.45 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 19.4 Khoảng tin cậy cho nhiều hệ số hồi quy
Khoảng tin cậy cho hai hệ số hồi quy là một hình ellipse mà trọng tâm ở điểm được định nghĩa bởi cả hai ước lượng. Trong R ta có thể xây dựng ellipse này bằng confidenceEllipse() từ package car.
# draw the 95% confidence set for coefficients on size and expenditure
confidenceEllipse(model,
fill = T,
lwd = 0,
which.coef = c("size", "expenditure"),
main = "95% Confidence Set")
Ta thấy bộ giá trị \((0,0)\) không nằm trong khoảng tin cậy này, do đó ta có \(95\%\) xác suất bác bỏ giả thuyết \(H_0:\beta_1 = 0, \beta_3=0\). Đoạn code trên dựa trên sai số chưa vững. Để điều chỉnh khoảng tin cậy cho phương sai sai số thay đổi, ta điều chỉnh đoạn code trên như sau.
# draw the robust 95% confidence set for coefficients on size and expenditure
confidenceEllipse(model,
fill = T,
lwd = 0,
which.coef = c("size", "expenditure"),
main = "95% Confidence Sets",
vcov. = vcovHC(model, type = "HC1"),
col = "red")
# draw the 95% confidence set for coefficients on size and expenditure
confidenceEllipse(model,
fill = T,
lwd = 0,
which.coef = c("size", "expenditure"),
add = T)
Vì sai số điều chỉnh phương sai thay đổi lớn hơn so với phiên bản chưa điều chỉnh nên khoảng tin cậy theo đó cũng rộng hơn.
9.5 Quy cách mô hình hồi quy đa biến
Chọn lựa quy cách mô hình, nghĩa là ta đang chọn biến giải thích đưa vào mô hình. Mục tiêu ở đây là: đạt được ước lượng nhân quả phi lệch và chính xác. Ta xem xét omitted variable bias đối với mô hình hồi quy đa biến. Giả sử ta xem xét mô hình hồi quy gồm hai biến giải thích là size và english. Giả sử ta quan tâm quan hệ nhân quả của size lên điểm thi. Quan hệ này có thể bị thiên lệch do ta đã bỏ qua cơ hội học tập bên ngoài. Những ba mẹ giàu có đủ khả năng về thời gian, tiền bạc cho việc học thêm của con trẻ. Do đó, ta có thể đánh giá khả năng này qua biến lunch, phần trăm học sinh được nhận buổi trưa miễn phí do thu nhập gia đình dưới ngưỡng quy định.
# estimate the model and print the summary to console
model <- lm(score ~ size + english + lunch, data = CASchools)
coeftest(model, vcov. = vcovHC, type = "HC1")##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 700.149957 5.568453 125.7351 < 2.2e-16 ***
## size -0.998309 0.270080 -3.6963 0.0002480 ***
## english -0.121573 0.032832 -3.7029 0.0002418 ***
## lunch -0.547345 0.024107 -22.7046 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ta thấy rằng, không có sự khác biệt nhiều về hệ số hồi quy biến size lên điểm thi ở cả hai mô hình.
9.5.1 Quy cách mô hình dựa vào \(R^2\) và \(\bar{R}^2\)
Cả hai chỉ tiêu đều cho ta biết liệu biến giải thích có giải thích tốt phương sai của biến được giải thích. Tuy nhiên, cả hai đều không cho ta biết:
- biến được thêm vào có ý nghĩa thống kê không.
- biến giải thích có là quan hệ nhân quả thật sự không
- liệu có omitted variable bias không
- bộ biến lựa chọn là tốt nhất chưa.
9.6 Phân tích bộ dữ liệu Điểm thi
Cho đến giờ, ta đã xem xét hai biến có ảnh hưởng đến điểm thi và có tương quan đến STR: english và lunch. Một biến khác cũng như vậy là calworks, thể hiện học sinh sống trong gia đình có tổng thu nhập dưới ngưỡng cho một bữa ăn trung bình. Như vậy cả hai biến calworks và lunch đều chỉ cùng một ý nghĩa và do đó, hệ số tương quan giữa chúng sẽ cao.
# estimate the correlation between 'calworks' and 'lunch'
cor(CASchools$calworks, CASchools$lunch)## [1] 0.7394218Khi sử dụng cả hai biến này trong mô hình, ắt hẳn sẽ mắc phải đa cộng tuyến. Ta xem xét các đồ thị sau.
# set up arrangement of plots
m <- rbind(c(1, 2), c(3, 0))
graphics::layout(mat = m)
# scatterplots
plot(score ~ english,
data = CASchools,
col = "steelblue",
pch = 20,
xlim = c(0, 100),
cex.main = 0.9,
main = "Percentage of English language learners")
plot(score ~ lunch,
data = CASchools,
col = "steelblue",
pch = 20,
cex.main = 0.9,
main = "Percentage qualifying for reduced price lunch")
plot(score ~ calworks,
data = CASchools,
col = "steelblue",
pch = 20,
xlim = c(0, 100),
cex.main = 0.9,
main = "Percentage qualifying for income assistance")
Ta thấy mọi quan hệ đều âm. Ta có 5 mô hình như sau.
- \(\text{Điểm thi} = \beta_0 + \beta_1\text{size}+u\)
- \(\text{Điểm thi} = \beta_0 + \beta_1\text{size} + \beta_2\text{english}+u\)
- \(\text{Điểm thi} = \beta_0 + \beta_1\text{size} + \beta_2\text{english} + \beta_3\text{lunch}+u\)
- \(\text{Điểm thi} = \beta_0 + \beta_1\text{size} + \beta_2\text{english} + \beta_4\text{calworks}+u\)
- \(\text{Điểm thi} = \beta_0 + \beta_1\text{size} + \beta_2\text{english} + \beta_3\text{lunch}+\beta_4\text{calworks}+u\)
| score | |||||
| (I) | (II) | (III) | (IV) | (V) | |
| (1) | (2) | (3) | (4) | (5) | |
| size | -2.280*** | -1.101** | -0.998*** | -1.308*** | -1.014*** |
| (0.519) | (0.433) | (0.270) | (0.339) | (0.269) | |
| english | -0.650*** | -0.122*** | -0.488*** | -0.130*** | |
| (0.031) | (0.033) | (0.030) | (0.036) | ||
| lunch | -0.547*** | -0.529*** | |||
| (0.024) | (0.038) | ||||
| calworks | -0.790*** | -0.048 | |||
| (0.068) | (0.059) | ||||
| Constant | 698.933*** | 686.032*** | 700.150*** | 697.999*** | 700.392*** |
| (10.364) | (8.728) | (5.568) | (6.920) | (5.537) | |
| Observations | 420 | 420 | 420 | 420 | 420 |
| R2 | 0.051 | 0.426 | 0.775 | 0.629 | 0.775 |
| Adjusted R2 | 0.049 | 0.424 | 0.773 | 0.626 | 0.773 |
| Residual Std. Error | 18.581 (df = 418) | 14.464 (df = 417) | 9.080 (df = 416) | 11.654 (df = 416) | 9.084 (df = 415) |
| F Statistic | 22.575*** (df = 1; 418) | 155.014*** (df = 2; 417) | 476.306*** (df = 3; 416) | 234.638*** (df = 3; 416) | 357.054*** (df = 4; 415) |
| Notes: | ***Significant at the 1 percent level. | ||||
| **Significant at the 5 percent level. | |||||
| *Significant at the 10 percent level. | |||||
Ta có thể rút ra kết luận như sau.
- Việc giảm STR một đơn vị dẫn đến điểm thi trung bình tăng lên 1 điểm.
- Tăng biến thì \(\bar{R}^2\) cũng tăng lên đáng kể, do đó ta có thể xem xét các biến này.
- Nếu đã bao gồm biến
lunchthì có thể bỏ biếncalworks.
9.7 Câu hỏi lập trình
9.7.1 Bài 1: Kiểm định mô hình hồi quy đa biến
Xem xét dữ liệu giá nhà Boston ở chương trước.
data('Boston')
mod <- lm(medv ~ lstat + crim + age, data = Boston)- Tính thống kê \(t\) và lưu chúng vào
tstats. - Tính \(p\)-value và lưu chúng vào
pval. - Ra quyết định thống kê.
9.7.2 Bài 2: Khoảng tin cậy mô hình hồi quy đa biến
Xây dựng khoảng tin cậy mọi hệ số hồi quy.
9.7.3 Bài 3: Kiểm định vững mô hình hồi quy đa biến
- Xây dựng sai số chuẩn vững.
- Ra quyết định thống kê vững.
9.7.4 Bài 4: Kiểm định gộp 1
Chúng ta quan tâm \(H_0: \beta_2 = \beta_3\) và \(H_1: \beta_2 \ne \beta_3\).
Ý tưởng ở đây là: thực hiện hai mô hình hồi quy và so sánh kết quả. Một mô hình giới hạn theo \(H_0\) và một mô hình đầy đủ. Từ hai kết quả, ta tính được giá trị thống kê \(F\).
Mô hình giới hạn như sau.
\[medv_i = \beta_0 + \beta_1 lstat_i + \beta_2 crim_i + \beta_2 age_i + u_i\]
Hướng dẫn:
- Ước lượng mô hình giới hạn, lưu vào biến
model_res. - Tính SSR của mô hình giới hạn, lưu vào
RSSR. - Ước lượng mô hình đầy đủ, lưu vào biến
model_unres. - Tính toán SSR cho mô hình đầy đủ, lưu vào
USSR.
9.7.5 Bài 5: Kiểm định gộp 2
- Tính thống kê \(F\), lưu vào biến
Fstat. - Tính \(p\)-value, lưu vào
pval. - Ra quyết định.
- Kiểm tra lại bằng hàm
linearHypothesis().
9.7.6 Bài 6: Khoảng tin cậy cho kiểm định gộp
Giả thuyết: \(H_0: \beta_2 = \beta_3 = 0\) và \(H_A: \beta_2 \ne 0 \text{ hoặc } \beta_3 \ne 0\).
- Xây dựng khoảng tin cậy.
- Hình hoạ hoá khoảng tin cậy.