Vẽ đồ thị hệ số trong hồi quy phân vị

Ở bài minh họa hồi quy phân vị, tôi đã trình bày cách ước lượng hồi quy phân vị, hồi quy trung vị bằng cách sử dụng các câu lệnh như qreg, iqreg, bsqreg và sqreg trên Stata. Tiếp tục chuỗi bài viết về hồi quy phân vị, bài viết hôm nay sẽ trình bày cách trực quan hóa kết quả hồi quy phân vị trên Stata qua câu lệnh qregplot. Cụ thể, trên Stata bạn có thể sử dụng lệnh qregplot để vẽ các đồ thị hệ số ước lượng từ kết quả ước lượng của các câu lệnh qreg, qreg2, bsqreg, mmqreg, qrprocess và cả rifhdreg. Đồ thị trong hồi quy phân vị Nội dung bên dưới sẽ trình bày cách sử dụng câu lệnh qregplot để trực quan hóa đồ thị các hệ số ước lượng từ kết quả hồi quy phân vị. Câu lệnh qregplot hỗ trợ cho nhiều câu lệnh ước lượng hồi quy phân vị. Để sử dụng câu lệnh, bạn cần cài đặt gói lệnh qregplot và các lệnh liên quan vào máy thông qua câu lệnh ssc install. Sau đây là một đoạn chương trình sử dụng lệnh vòng lặp foreach để cài đặt những gói lệnh cần thiết cho hồi quy phân vị và các công cụ trực quan hóa kết quả ước lượng. local packages = "qreg sqreg bsqreg mmqreg ivqreg2 qreg2 xtqreg /// ftools gtools qregplot qrprocess" foreach i of local packages { capture which `i' if _rc != 0 { display _rc ssc install `i', replace } } Trực quan hóa kết quả Dữ liệu thực hành Sử dụng mẫu dữ liệu nhỏ sẵn có womenwk trên Stata . webuse womenwk, clear . set scheme white Kết quả từ lệnh qreg Giả sử, chúng ta muốn ước lượng một mô hình hồi quy phân vị của các mức lương, trong đó lương (wage) là một hàm của tuổi (age), trình độ học vấn (education), tình trạng hôn nhân (married) và khu vực sinh sống (county). Trung vị của lương trong mô hình được ước lượng như sau: . qreg wage age education married i.county Iteration 1: WLS sum of weighted deviations = 2882.0578 Iteration 1: sum of abs. weighted deviations = 2881.6394 note: alternate solutions exist. Iteration 2: sum of abs. weighted deviations = 2880.7377 Iteration 3: sum of abs. weighted deviations = 2880.6178 Iteration 4: sum of abs. weighted deviations = 2879.5666 Iteration 5: sum of abs. weighted deviations = 2879.3427 Iteration 6: sum of abs. weighted deviations = 2878.958 Iteration 7: sum of abs. weighted deviations = 2878.8413 Iteration 8: sum of abs. weighted deviations = 2878.8092 Iteration 9: sum of abs. weighted deviations = 2878.6626 Iteration 10: sum of abs. weighted deviations = 2878.4107 note: alternate solutions exist. Iteration 11: sum of abs. weighted deviations = 2878.352 Iteration 12: sum of abs. weighted deviations = 2878.2839 Iteration 13: sum of abs. weighted deviations = 2878.2243 Iteration 14: sum of abs. weighted deviations = 2878.2232 Iteration 15: sum of abs. weighted deviations = 2878.1574 Iteration 16: sum of abs. weighted deviations = 2878.1558 Iteration 17: sum of abs. weighted deviations = 2878.1395 Iteration 18: sum of abs. weighted deviations = 2878.1377 Median regression Number of obs = 1,343 Raw sum of deviations 3394.897 (about 23.511223) Min sum of deviations 2878.138 Pseudo R2 = 0.1522 ------------------------------------------------------------------------------ wage | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- age | .1910125 .0279023 6.85 0.000 .1362752 .2457498 education | .9196805 .0698232 13.17 0.000 .7827049 1.056656 married | -.23724 .5409577 -0.44 0.661 -1.298463 .8239834 | county | 1 | -2.047806 1.025437 -2.00 0.046 -4.059456 -.0361563 2 | -2.53552 1.027812 -2.47 0.014 -4.551829 -.5192117 3 | -2.915096 1.042238 -2.80 0.005 -4.959705 -.8704865 4 | -3.315867 1.044933 -3.17 0.002 -5.365764 -1.26597 5 | -1.505194 1.055567 -1.43 0.154 -3.575951 .5655644 6 | -2.539745 1.084961 -2.34 0.019 -4.668167 -.4113222 7 | -2.373452 1.102168 -2.15 0.031 -4.535629 -.211274 8 | -3.471071 1.133828 -3.06 0.002 -5.695357 -1.246786 9 | -2.857953 1.143551 -2.50 0.013 -5.101313 -.6145929 | _cons | 6.752149 1.369943 4.93 0.000 4.064664 9.439635 ------------------------------------------------------------------------------ Bây giờ, chúng ta sử dụng câu lệnh qregplot để biểu diễn các hệ số ước lượng theo phân phối của các biến age, education và married. Bạn có thể đặt tên biến để lưu lại các kết quả đồ thị và sử dụng lại khi cần qua tùy chọn estore() của câu lệnh. Qua tùy chọn quantile() bạn có thể thiết lập các phân vị muốn hiển thị ở trục hoành. . qregplot age education married, estore(gr1) quantile(5(5)95) name(qregplot0, replace) macros: r(yvar) : "wage" r(xvar) : "age education married i.county" . graph export qregplot0.png, as(png) replace name("qregplot0") Nếu muốn hiển thị nhãn biến thay cho tên biến ở đồ thị thì bạn có thể sử dụng tùy chọn label. Hoặc thực hiện thêm một số trang trí cho đồ thị như số cột hiển thị, tiêu đề (tùy chọn mtitles)... Bạn có thể tham khảo thêm ở phần help qregplot. Ở đây, tôi gọi lại đồ thị gr1 đã lưu ở bước trên. . qregplot age education married, from(gr1) /// label name(qregplot1, replace) /// mtitles("Age in years since 1980" /// "Years of education, including Highschool and college") /// labelopt(lines(2)) //hiển thị trên 2 dòng khi tiêu đề dài (results gr1 are active now) . graph export qregplot1.png, as(png) replace name("qregplot1") Sau các lệnh khác Ngoài lệnh qreg, qregplot còn hỗ trợ trực quan hóa các hệ số ước lượng được ước lượng từ kết quả hồi quy phân vị của nhiều câu lệnh khác như bsqreg, qreg2, qrprocess và mmqreg hay cả rifhdreg. Kết quả từ lệnh bsqreg . bsqreg wage age education married i.county, reps(25) (fitting base model) Bootstrap replications (25) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 ......................... Median regression, bootstrap(25) SEs Number of obs = 1,343 Raw sum of deviations 3394.897 (about 23.511223) Min sum of deviations 2878.138 Pseudo R2 = 0.1522 ------------------------------------------------------------------------------ wage | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- age | .1910125 .037505 5.09 0.000 .1174372 .2645879 education | .9196805 .0694333 13.25 0.000 .7834697 1.055891 married | -.23724 .5904457 -0.40 0.688 -1.395546 .9210664 | county | 1 | -2.047806 1.321545 -1.55 0.121 -4.640346 .544734 2 | -2.53552 .881511 -2.88 0.004 -4.264824 -.8062166 3 | -2.915096 1.054119 -2.77 0.006 -4.983014 -.8471784 4 | -3.315867 1.294298 -2.56 0.011 -5.854956 -.7767784 5 | -1.505194 1.468713 -1.02 0.306 -4.386439 1.376052 6 | -2.539745 1.354207 -1.88 0.061 -5.196359 .1168697 7 | -2.373452 1.336936 -1.78 0.076 -4.996184 .2492807 8 | -3.471071 1.33976 -2.59 0.010 -6.099343 -.8427988 9 | -2.857953 1.382681 -2.07 0.039 -5.570427 -.1454788 | _cons | 6.752149 1.300033 5.19 0.000 4.201811 9.302488 ------------------------------------------------------------------------------ . qregplot age education, q(5(5)95) seed(101) /// label title("Kết quả từ lệnh bsqreg") /// name(bsqreg, replace) ysize(4) xsize(7) macros: r(oth) : "reps(25)" r(yvar) : "wage" r(xvar) : "age education married i.county" . graph export bsqreg.png, as(png) replace name("bsqreg") Kết quả từ lệnh qreg2 . qreg2 wage age education married i.county Median regression R-squared = .26371107 Number of obs = 1343 Objective function = 2.1430661 Heteroskedasticity robust standard errors ------------------------------------------------------------------------------ wage | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- age | .1910125 .030997 6.16 0.000 .1302042 .2518209 education | .9196805 .0760262 12.10 0.000 .7705362 1.068825 married | -.23724 .5815003 -0.41 0.683 -1.377998 .9035178 | county | 1 | -2.047806 1.201217 -1.70 0.088 -4.404292 .3086801 2 | -2.53552 1.185086 -2.14 0.033 -4.860362 -.2106785 3 | -2.915096 1.126806 -2.59 0.010 -5.125606 -.7045859 4 | -3.315867 1.175865 -2.82 0.005 -5.622619 -1.009115 5 | -1.505194 1.203101 -1.25 0.211 -3.865376 .8549893 6 | -2.539745 1.274159 -1.99 0.046 -5.039325 -.0401641 7 | -2.373452 1.271166 -1.87 0.062 -4.86716 .1202565 8 | -3.471071 1.313158 -2.64 0.008 -6.047159 -.8949836 9 | -2.857953 1.337014 -2.14 0.033 -5.480838 -.2350673 | _cons | 6.752149 1.469769 4.59 0.000 3.868832 9.635467 ------------------------------------------------------------------------------ Machado-Santos Silva test for heteroskedasticity Ho: Constant variance Variables: Fitted values of wage and its squares chi2(2) = 0.854 Prob > chi2 = 0.652 . qregplot age education, q(5(5)95) seed(101) /// label title("Kết quả từ lệnh qreg2") /// name(qreg2, replace) ysize(4) xsize(7) macros: r(yvar) : "wage" r(xvar) : "age education married i.county" . graph export qreg2.png, as(png) replace name("qreg2") Kết quả từ lệnh qrprocess . qrprocess wage age education married i.county Quantile regression No. of obs. 1343 Algorithm: qreg. Variance: kernel estimate of the sandwich as proposed by Powell(1990). ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- Quant. 0.5 | age | .1910125 .0294435 6.49 0.000 .1332518 .2487733 education | .9196805 .0740152 12.43 0.000 .7744812 1.06488 married | -.23724 .5548447 -0.43 0.669 -1.325706 .8512261 1.county | -2.047806 1.058309 -1.93 0.053 -4.123943 .0283304 2.county | -2.53552 1.070409 -2.37 0.018 -4.635393 -.4356471 3.county | -2.915096 1.031819 -2.83 0.005 -4.939267 -.8909249 4.county | -3.315867 1.068975 -3.10 0.002 -5.412929 -1.218806 5.county | -1.505194 1.076199 -1.40 0.162 -3.616427 .6060397 6.county | -2.539745 1.141561 -2.22 0.026 -4.779201 -.3002885 7.county | -2.373452 1.169033 -2.03 0.043 -4.666802 -.0801013 8.county | -3.471071 1.187918 -2.92 0.004 -5.801468 -1.140674 9.county | -2.857953 1.207965 -2.37 0.018 -5.227677 -.4882282 _cons | 6.752149 1.378011 4.90 0.000 4.048838 9.455461 ------------------------------------------------------------------------------ . qregplot age education, q(5(5)95) seed(101) /// label title("Kết quả từ lệnh qrprocess") /// name(qrprocess, replace) ysize(4) xsize(7) macros: r(yvar) : "wage" r(xvar) : "age education married i.county" . graph export qrprocess.png, as(png) replace name("qrprocess") Kết quả từ lệnh mmqreg . mmqreg wage age education married, abs(county) MM-qreg Estimator Number of obs = 1343 Quantile: 50 ------------------------------------------------------------------------------ wage | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- location | age | .1749932 .0216822 8.07 0.000 .1324969 .2174896 education | .9527413 .0529194 18.00 0.000 .8490212 1.056461 married | .2235909 .3988231 0.56 0.575 -.558088 1.00527 _cons | 4.093984 1.124378 3.64 0.000 1.890244 6.297724 -------------+---------------------------------------------------------------- scale | age | .0178025 .0130512 1.36 0.173 -.0077773 .0433823 education | -.0174842 .0318538 -0.55 0.583 -.0799164 .0449481 married | .5608084 .2400634 2.34 0.019 .0902928 1.031324 _cons | 3.451374 .6767961 5.10 0.000 2.124878 4.77787 -------------+---------------------------------------------------------------- qtile | age | .1748125 .0216847 8.06 0.000 .1323113 .2173136 education | .9529188 .0529112 18.01 0.000 .8492148 1.056623 married | .2178976 .399112 0.55 0.585 -.5643476 1.000143 _cons | 4.058946 1.128454 3.60 0.000 1.847217 6.270674 ------------------------------------------------------------------------------ . qregplot age education, q(5(5)95) seed(101) /// label title("Kết quả từ lệnh mmqreg") /// name(mmqreg, replace) ysize(4) xsize(7) macros: r(oth) : "abs(county)" r(yvar) : "wage" r(xvar) : "age education married" . graph export mmqreg.png, as(png) replace name("mmqreg") Kết quả từ lệnh rifhdreg . rifhdreg wage age education married, abs(county) rif(q(50)) HDFE Linear regression Number of obs = 1,343 Absorbing 1 HDFE group F( 3, 1330) = 73.52 Prob > F = 0.0000 R-squared = 0.1772 Adj R-squared = 0.1698 Within R-sq. = 0.1422 Root MSE = 7.7354 ------------------------------------------------------------------------------ wage | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- age | .1985129 .0304391 6.52 0.000 .1387989 .2582269 education | 1.026049 .0761715 13.47 0.000 .8766194 1.175478 married | -.1943258 .5901412 -0.33 0.742 -1.352035 .9633833 _cons | 2.348411 1.610248 1.46 0.145 -.8104913 5.507312 ------------------------------------------------------------------------------ Absorbed degrees of freedom: -----------------------------------------------------+ Absorbed FE | Categories - Redundant = Num. Coefs | -------------+---------------------------------------| county | 10 0 10 | -----------------------------------------------------+ Distributional Statistic: q(50) Sample Mean RIF q(50) : 23.518 . qregplot age education, q(5(5)95) seed(101) /// label title("Kết quả từ lệnh rifhdreg") /// name(rifhdreg, replace) ysize(4) xsize(7) . graph export rifhdreg.png, as(png) replace name("rifhdreg") Bên cạnh câu lệnh qregplot, Stata còn có các câu lệnh khác để vẽ đồ thị các hệ số của hồi quy phân vị như grqreg hay rqrplot. Ở đây, câu lệnh grqreg chỉ hoạt động sau các lệnh ước lượng như qreg, bsqreg và sqreg. Xem thêm help grqreg để hiểu hơn về câu lệnh và thực hành qua nhiều ví dụ. . sysuse auto, clear (1978 automobile data) . qreg price weight length foreign, quantile(0.5) Iteration 1: WLS sum of weighted deviations = 56397.829 Iteration 1: sum of abs. weighted deviations = 55950.5 Iteration 2: sum of abs. weighted deviations = 55264.718 Iteration 3: sum of abs. weighted deviations = 54762.283 Iteration 4: sum of abs. weighted deviations = 54734.152 Iteration 5: sum of abs. weighted deviations = 54552.638 note: alternate solutions exist. Iteration 6: sum of abs. weighted deviations = 54465.511 Iteration 7: sum of abs. weighted deviations = 54443.699 Iteration 8: sum of abs. weighted deviations = 54411.294 Median regression Number of obs = 74 Raw sum of deviations 71102.5 (about 4934) Min sum of deviations 54411.29 Pseudo R2 = 0.2347 ------------------------------------------------------------------------------ price | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- weight | 3.933588 1.328718 2.96 0.004 1.283543 6.583632 length | -41.25191 45.46469 -0.91 0.367 -131.9284 49.42456 foreign | 3377.771 885.4198 3.81 0.000 1611.857 5143.685 _cons | 344.6489 5182.394 0.07 0.947 -9991.31 10680.61 ------------------------------------------------------------------------------ . quietly grqreg foreign, ci ols olsci graphregion(color(white)) name(g1, replace) . bsqreg price weight length foreign, quantile(0.5) reps(100) (fitting base model) Bootstrap replications (100) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 50 .................................................. 100 Median regression, bootstrap(100) SEs Number of obs = 74 Raw sum of deviations 71102.5 (about 4934) Min sum of deviations 54411.29 Pseudo R2 = 0.2347 ------------------------------------------------------------------------------ price | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- weight | 3.933588 2.742986 1.43 0.156 -1.537125 9.4043 length | -41.25191 72.94675 -0.57 0.574 -186.7396 104.2358 foreign | 3377.771 1065.78 3.17 0.002 1252.139 5503.403 _cons | 344.6489 6226.547 0.06 0.956 -12073.81 12763.11 ------------------------------------------------------------------------------ . quietly grqreg foreign, ci ols olsci graphregion(color(white)) name(g2, replace) . graph combine g1 g2, name(grcbine, replace) . graph export grcbine.png, as(png) replace name("grcbine") . grc1leg g1 g2, name(grc1leg, replace) . graph export grc1leg.png, as(png) replace name("grc1leg") Tôi sử dụng lệnh grc1leg để kết hợp các đồ thị riêng rẽ vào cùng một đồ thị với chung một legend để khắc phục điều mà lệnh graph combine không làm được (2 legend riêng). Trong trường hợp bạn lệnh grc1leg báo lỗi r(4018) rằng Graph.graphs[1].legend.draw_view.set_false: class type not found thì bạn có thể thêm đoạn lệnh sau vào trước dòng #delimit cr (ở gần cuối file) vào file C:\ado\plus\g\grqreg.ado, lưu lại và khởi động lại Stata. *-> combine all graphs; if ("`nodraw'"=="nodraw") {; grc1leg `graphlist', nodraw `options'; // graph combine `graphlist', nodraw `options'; }; if ("`nodraw'"=="") {; grc1leg `graphlist', `options'; // graph combine `graphlist', `options'; }; Sau đó chạy lại đoạn lệnh trên. Sở dĩ phải điều chỉnh vậy là câu lệnh bản thân trong grqreg đã kết hợp các đồ thị nên grc1leg không thể thực hiện.