1 THỰC HÀNH

1.1 Chạy hồi quy tuyến tính đa biến

cụ thể chi tiết:

  • Biến phụ thuộc: thu nhập (income)

  • Biến độc lập: women, education, điểm uy tín nghề nghiệp lần lượt là trình độ học vấn, tỷ lệ người đi làm là phụ nữ

1.1.1 Ước lượng mô hình

library(AER)
data("Prestige")
pob <- lm(income ~ women + education, data = Prestige)
summary(pob)
## 
## Call:
## lm(formula = income ~ women + education, data = Prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7257.6 -1160.1  -238.6   681.1 16044.3 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1491.998   1162.299  -1.284    0.202    
## women         -64.056      8.921  -7.180 1.31e-10 ***
## education     944.881    103.731   9.109 9.60e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2839 on 99 degrees of freedom
## Multiple R-squared:  0.5618, Adjusted R-squared:  0.5529 
## F-statistic: 63.46 on 2 and 99 DF,  p-value: < 2.2e-16

Nhận xét: Hệ số R bình phương hiệu chỉnh bằng 0.5549 nghĩa là 55.49% sự biến thiên trong thu nhập được giải thích bởi các biến trình độ học vấn, tỷ lệ người đi làm là phụ nữ tức là thu nhập của 1 người còn phụ thuộc vào các yếu tố trên để quyết định đến mức lương của người đó.

1.2 Tìm khoảng tin cậy cho các hệ số hồi quy

Sử dụng hàm confint()

  • độ tin cậy là 95%
confint(pob) 
##                   2.5 %     97.5 %
## (Intercept) -3798.25140  814.25546
## women         -81.75738  -46.35407
## education     739.05564 1150.70689
  • độ tin cậy là 99%
confint(pob, level = 0.99)
##                   0.5 %     99.5 %
## (Intercept) -4544.66648 1560.67054
## women         -87.48649  -40.62496
## education     672.44052 1217.32201

1.3 Kiểm định và lựa chọn mô hình

1.3.1 Dựa vào đồ thị phần dư

par(mfrow = c(1, 1))
plot(pob)

plot(pob, 1)

1.3.2 Dựa vào các kiểm định

re <- resid(pob)
  • Kiểm định 1: sai số ngẫu nhiên có phân phối chuẩn

Giả thiết kiểm định:

  • H0: Dữ liệu có phân phối chuẩn
  • H1: Dữ liệu không có phân phối chuẩn
shapiro.test(re)
## 
##  Shapiro-Wilk normality test
## 
## data:  re
## W = 0.75852, p-value = 1.174e-11
  • Kiểm định 2: Kỳ vọng của sai số ngẫu nhiên tại mỗi giá trị bằng 0
t.test(re, mu = 0)
## 
##  One Sample t-test
## 
## data:  re
## t = 6.2875e-16, df = 101, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -552.0689  552.0689
## sample estimates:
##    mean of x 
## 1.749799e-13
  • Kiểm định 3: Phương sai của sai số ngẫu nhiên không đổi H0: phương sai sai số không đổi

H1: Phương sai sai số thay đổi

ncvTest(pob)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 49.8013, Df = 1, p = 1.7013e-12
  • Kiểm định 4: Giữa các biến độc lập không có mối quan hệ đa cộng tuyến hoàn hảo
vif(pob)
##     women education 
##   1.00384   1.00384

1.4 Dự báo

women <- c(0.78, 1.65, 2.46, 3.69, 4.14, 5.13, 6.01, 7.83, 8.13, 9.11)
education <- c(6.38, 7.11, 8.24, 9.05, 10.05, 11.09, 12.39, 13.62, 14.15, 15.21)
new <- data.frame(women, education)

1.4.1 Dự báo giá trị trung bình

predict(pob , newdata = new, interval = "confidence")
##          fit       lwr       upr
## 1   4486.381  3342.059  5630.703
## 2   5120.416  4091.796  6149.036
## 3   6136.247  5261.056  7011.437
## 4   6822.812  6039.893  7605.731
## 5   7738.868  7020.000  8457.736
## 6   8658.129  7952.193  9364.066
## 7   9830.106  9049.646 10610.566
## 8  10875.729  9964.367 11787.090
## 9  11357.299 10371.124 12343.473
## 10 12296.098 11146.291 13445.906

1.4.2 Dự báo giá trị cá biệt

predict(pob, newdata = new, interval = "prediction")
##          fit        lwr      upr
## 1   4486.381 -1261.7156 10234.48
## 2   5120.416  -605.7697 10846.60
## 3   6136.247   435.6237 11836.87
## 4   6822.812  1135.6241 12510.00
## 5   7738.868  2060.1434 13417.59
## 6   8658.129  2981.0272 14335.23
## 7   9830.106  4143.2562 15516.96
## 8  10875.729  5169.4405 16582.02
## 9  11357.299  5638.5854 17076.01
## 10 12296.098  6546.9071 18045.29
LS0tDQp0aXRsZTogInRodWNoYW5oIg0KYXV0aG9yOiAicHRwaHVvbmciDQpkYXRlOiAiMjAyMy0wNy0yNSINCm91dHB1dDogDQogaHRtbF9kb2N1bWVudDogDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KICAgIGNvZGVfZm9sZGluZzogaGlkZQ0KICAgIGhpZ2hsaWdodDogcHlnbWVudHMNCiAgICBudW1iZXJfc2VjdGlvbnM6IHllcw0KICAgIHRoZW1lOiAiZmxhdGx5Ig0KICAgIHRvYzogVFJVRQ0KICAgIHRvY19mbG9hdDogVFJVRQ0KLS0tDQoNCiMgVEjhu7BDIEjDgE5IDQoNCiMjIENo4bqheSBo4buTaSBxdXkgdHV54bq/biB0w61uaCDEkWEgYmnhur9uDQoNCiBj4bulIHRo4buDIGNoaSB0aeG6v3Q6DQoNCiogQmnhur9uIHBo4bulIHRodeG7mWM6IHRodSBuaOG6rXAgKGluY29tZSkgDQoNCiogQmnhur9uIMSR4buZYyBs4bqtcDogd29tZW4sIGVkdWNhdGlvbiwgxJFp4buDbSB1eSB0w61uIG5naOG7gSBuZ2hp4buHcCBs4bqnbiBsxrDhu6N0IGzDoCB0csOsbmggxJHhu5kgaOG7jWMgduG6pW4sIHThu7cgbOG7hyBuZ8aw4budaSDEkWkgbMOgbSBsw6AgcGjhu6UgbuG7rw0KDQojIyMgxq/hu5tjIGzGsOG7o25nIG3DtCBow6xuaA0KDQpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQ0KbGlicmFyeShBRVIpDQpkYXRhKCJQcmVzdGlnZSIpDQpwb2IgPC0gbG0oaW5jb21lIH4gd29tZW4gKyBlZHVjYXRpb24sIGRhdGEgPSBQcmVzdGlnZSkNCnN1bW1hcnkocG9iKQ0KYGBgDQpOaOG6rW4geMOpdDogSOG7hyBz4buRIFIgYsOsbmggcGjGsMahbmcgaGnhu4d1IGNo4buJbmggYuG6sW5nIDAuNTU0OSBuZ2jEqWEgbMOgIDU1LjQ5JSBz4buxIGJp4bq/biB0aGnDqm4gdHJvbmcgdGh1IG5o4bqtcCDEkcaw4bujYyBnaeG6o2kgdGjDrWNoIGLhu59pIGPDoWMgYmnhur9uIHRyw6xuaCDEkeG7mSBo4buNYyB24bqlbiwgdOG7tyBs4buHIG5nxrDhu51pIMSRaSBsw6BtIGzDoCBwaOG7pSBu4buvIHThu6ljIGzDoCB0aHUgbmjhuq1wIGPhu6dhIDEgbmfGsOG7nWkgY8OybiBwaOG7pSB0aHXhu5ljIHbDoG8gY8OhYyB54bq/dSB04buRIHRyw6puIMSR4buDIHF1eeG6v3QgxJHhu4tuaCDEkeG6v24gbeG7qWMgbMawxqFuZyBj4bunYSBuZ8aw4budaSDEkcOzLg0KDQoNCiMjIFTDrG0ga2hv4bqjbmcgdGluIGPhuq15IGNobyBjw6FjIGjhu4cgc+G7kSBo4buTaSBxdXkNCg0KU+G7rSBk4bulbmcgaMOgbSAqKmNvbmZpbnQoKSoqDQoNCiogxJHhu5kgdGluIGPhuq15IGzDoCA5NSUNCg0KYGBge3J9DQpjb25maW50KHBvYikgDQpgYGANCg0KKiDEkeG7mSB0aW4gY+G6rXkgbMOgIDk5JQ0KDQpgYGB7cn0NCmNvbmZpbnQocG9iLCBsZXZlbCA9IDAuOTkpDQoNCmBgYA0KDQojIyBLaeG7g20gxJHhu4tuaCB2w6AgbOG7sWEgY2jhu41uIG3DtCBow6xuaA0KDQojIyMgROG7sWEgdsOgbyDEkeG7kyB0aOG7iyBwaOG6p24gZMawDQoNCmBgYHtyfQ0KcGFyKG1mcm93ID0gYygxLCAxKSkNCnBsb3QocG9iKQ0KYGBgDQoNCg0KYGBge3J9DQpwbG90KHBvYiwgMSkNCmBgYA0KDQojIyMgROG7sWEgdsOgbyBjw6FjIGtp4buDbSDEkeG7i25oDQoNCmBgYHtyfQ0KcmUgPC0gcmVzaWQocG9iKQ0KYGBgDQoNCiogS2nhu4NtIMSR4buLbmggMTogc2FpIHPhu5Egbmfhuqt1IG5oacOqbiBjw7MgcGjDom4gcGjhu5FpIGNodeG6qW4NCg0KR2nhuqMgdGhp4bq/dCBraeG7g20gxJHhu4tuaDoNCg0KIC0gSDA6IEThu68gbGnhu4d1IGPDsyBwaMOibiBwaOG7kWkgY2h14bqpbg0KIC0gSDE6IEThu68gbGnhu4d1IGtow7RuZyBjw7MgcGjDom4gcGjhu5FpIGNodeG6qW4NCg0KYGBge3J9DQpzaGFwaXJvLnRlc3QocmUpDQpgYGANCiogS2nhu4NtIMSR4buLbmggMjogS+G7syB24buNbmcgY+G7p2Egc2FpIHPhu5Egbmfhuqt1IG5oacOqbiB04bqhaSBt4buXaSBnacOhIHRy4buLIGLhurFuZyAwDQoNCmBgYHtyfQ0KdC50ZXN0KHJlLCBtdSA9IDApDQpgYGANCg0KKiBLaeG7g20gxJHhu4tuaCAzOiBQaMawxqFuZyBzYWkgY+G7p2Egc2FpIHPhu5Egbmfhuqt1IG5oacOqbiBraMO0bmcgxJHhu5VpDQpIMDogcGjGsMahbmcgc2FpIHNhaSBz4buRIGtow7RuZyDEkeG7lWkNCg0KSDE6IFBoxrDGoW5nIHNhaSBzYWkgc+G7kSB0aGF5IMSR4buVaQ0KDQpgYGB7cn0NCm5jdlRlc3QocG9iKQ0KYGBgDQoqIEtp4buDbSDEkeG7i25oIDQ6IEdp4buvYSBjw6FjIGJp4bq/biDEkeG7mWMgbOG6rXAga2jDtG5nIGPDsyBt4buRaSBxdWFuIGjhu4cgxJFhIGPhu5luZyB0dXnhur9uIGhvw6BuIGjhuqNvDQoNCmBgYHtyfQ0KdmlmKHBvYikNCg0KYGBgDQojIyBE4buxIGLDoW8NCg0KYGBge3J9DQp3b21lbiA8LSBjKDAuNzgsIDEuNjUsIDIuNDYsIDMuNjksIDQuMTQsIDUuMTMsIDYuMDEsIDcuODMsIDguMTMsIDkuMTEpDQplZHVjYXRpb24gPC0gYyg2LjM4LCA3LjExLCA4LjI0LCA5LjA1LCAxMC4wNSwgMTEuMDksIDEyLjM5LCAxMy42MiwgMTQuMTUsIDE1LjIxKQ0KbmV3IDwtIGRhdGEuZnJhbWUod29tZW4sIGVkdWNhdGlvbikNCmBgYA0KDQojIyMgROG7sSBiw6FvIGdpw6EgdHLhu4sgdHJ1bmcgYsOsbmgNCg0KYGBge3J9DQpwcmVkaWN0KHBvYiAsIG5ld2RhdGEgPSBuZXcsIGludGVydmFsID0gImNvbmZpZGVuY2UiKQ0KYGBgDQoNCiMjIyBE4buxIGLDoW8gZ2nDoSB0cuG7iyBjw6EgYmnhu4d0DQoNCmBgYHtyfQ0KcHJlZGljdChwb2IsIG5ld2RhdGEgPSBuZXcsIGludGVydmFsID0gInByZWRpY3Rpb24iKQ0KYGBgDQoNCg0K