Estimasi dan Pelaporan

Estimasi Model Regresi

Menjalankan regresi OLS pada R menggunakan perintah lm. Berikut adalah regresi tingkat morbiditas terhadap jumlah dokter, bidan, dan pendapatan domestik bruto per kapita di tiap kabupaten, menggunakan data INDO-DAPOER dari World Bank yang sudah disiapkan di bab sebelumnya dan tersimpan di file datakab.rds.

Variabel pendapatan domestik bruto per kapita tidak terdapat di data tersebut, namun dapat dihitung sendiri dengan membagi variabel pendapatan domestik bruto dengan total populasi. Perhitungan itu bisa langsung dilakukan dalam spesifikasi persamaan pada fungsi lm dengan menempatkan perhitungan dalam I().

library(tidyverse)
datakab = readRDS('datakab.rds')

lm(morbidity_rate ~ number_of_doctors + number_of_midwives + 
                    I(total_gdp_based_on_expenditure/total_population), 
   data = datakab)

Call:
lm(formula = morbidity_rate ~ number_of_doctors + number_of_midwives + 
    I(total_gdp_based_on_expenditure/total_population), data = datakab)

Coefficients:
                                       (Intercept)  
                                        28.9713794  
                                 number_of_doctors  
                                        -0.0009816  
                                number_of_midwives  
                                         0.0002814  
I(total_gdp_based_on_expenditure/total_population)  
                                         0.0135187  

Jika variabel hasil perhitungan akan digunakan lagi dalam analisis lainnya, lebih baik ditambahkan sebagai variabel baru dalam data.

datakab = datakab %>% mutate(gdp_per_capita = total_gdp_based_on_expenditure/total_population)
reg_morbid = lm(morbidity_rate ~ number_of_doctors + number_of_midwives + gdp_per_capita, datakab)
reg_morbid

Call:
lm(formula = morbidity_rate ~ number_of_doctors + number_of_midwives + 
    gdp_per_capita, data = datakab)

Coefficients:
       (Intercept)   number_of_doctors  number_of_midwives      gdp_per_capita  
        28.9713794          -0.0009816           0.0002814           0.0135187  

Output fungsi lm bisa disimpan menjadi obyek tipe/kelas lm, yang bisa digunakan untuk proses analisis berikutnya.

class(reg_morbid)

Obyek lm merupakan kumpulan dari vektor dan matriks yang menyimpan hasil estimasi, nilai residual dan prediksi variabel respon, dan bagian data dari variabel-variabel yang dilibatkan dalam estimasi.

Komponen obyek lm bisa didaftar namanya dengan perintah names dan ditampilkan dengan operator $.

names(reg_morbid) %>% t

reg_morbid$coefficients %>% t

Catatan: Kedua perintah di atas bisa dijalankan tanpa piping ke fungsi transpose ` %>% t`, yang digunakan di sini semata untuk menghemat ruang untuk menampilkan output R dan membuatnya lebih mudah dibaca.

tail(reg_morbid$model)

Pelaporan Hasil Regresi

Ada berapa library R yang memudahkan pelaporan hasil regresi beserta statistik terkait, seperti stargazer dan texreg. Library tersebut juga memfasilitasi perbandingan antar model dengan menampilkan hasil dari beberapa regresi secara berdampingan.

Untuk menampilkan laporan di console R, gunakan argumen type = 'text' di stargazer dan fungsi screenreg untuk library texreg.

Piping hasil screenreg ke print tidak diperlukan di console, hanya perlu ditambahkan untuk penampilan di jupyter notebook.

reg_morbid_2 = lm(morbidity_rate ~ number_of_doctors + number_of_midwives, datakab)

list(reg_morbid, reg_morbid_2) %>% texreg::screenreg() %>% print

===========================================
                    Model 1     Model 2    
-------------------------------------------
(Intercept)          28.97 ***    29.11 ***
                     (1.03)       (0.22)   
number_of_doctors    -0.00        -0.00 *  
                     (0.00)       (0.00)   
number_of_midwives    0.00         0.00    
                     (0.00)       (0.00)   
gdp_per_capita        0.01                 
                     (0.04)                
-------------------------------------------
R^2                   0.03         0.00    
Adj. R^2             -0.00         0.00    
Num. obs.           110         1825       
RMSE                  6.71         8.78    
===========================================
*** p < 0.001, ** p < 0.01, * p < 0.05

stargazer::stargazer(reg_morbid, reg_morbid_2, type = 'text')

============================================================
                              Dependent variable:           
                    ----------------------------------------
                                 morbidity_rate             
                            (1)                 (2)         
------------------------------------------------------------
number_of_doctors         -0.001              -0.001**      
                          (0.001)             (0.001)       
                                                            
number_of_midwives        0.0003              0.0003*       
                         (0.0002)             (0.0002)      
                                                            
gdp_per_capita             0.014                            
                          (0.039)                           
                                                            
Constant                 28.971***           29.111***      
                          (1.032)             (0.217)       
                                                            
------------------------------------------------------------
Observations                110                1,825        
R2                         0.026               0.002        
Adjusted R2               -0.002               0.001        
Residual Std. Error  6.708 (df = 106)    8.778 (df = 1822)  
F Statistic         0.928 (df = 3; 106) 2.219 (df = 2; 1822)
============================================================
Note:                            *p<0.1; **p<0.05; ***p<0.01

stargazer dan texreg juga menyediakan format output latex dan html.

• untuk output latex, gunakan argumen type = tex di fungsi stargazer, dan fungsi texreg di library texreg
• untuk output html, gunakan argumen type = html di fungsi stargazer, dan fungsi htmlreg di library texreg

Format latex dan html sulit dibaca di layar. Tampilan di layar perlu disalin dulu ke file dokumen dan ditampilkan sebagai pdf untuk latex atau melalui browser untuk html. Output latex dan html bisa juga langsung disimpan ke file dengan menyebutkan path dan nama file pada argumen out untuk stargazer dan file untuk texreg.

stargazer otomatis mengenali tipe output yang diinginkan dari ekstensi nama file, sehingga tidak perlu menyebutkan argumen type lagi.

texreg::texreg(list(reg_morbid, reg_morbid_2), file = 'reg_morbid.tex')
stargazer::stargazer(reg_morbid, reg_morbid_2, out = 'reg_morbid.html')

Untuk menampilkan output di dokumen format doc/docx atau ppt/pptx, bisa dengan menyalin tampilan output html di browser. Output html bisa juga ditampilkan di jupyter notebook dengan bantuan IRdisplay::display_html

IRdisplay::display_html(file = 'reg_morbid.html')

Robust Standar Error

Tabel hasil regresi dalam jurnal ilmiah biasanya melaporkan standar error yang telah dikoreksi dari masalah heterokesdastisitas dan atau korelasi serial. Standard error yang dikoreksi tersebut bisa dihitung dari obyek lm yang sudah ada, lalu digunakan untuk mengganti standard error bawaan OLS. Namun cara yang lebih mudah adalah dengan sejak awal estimasi menggunakan fungsi lm_robust dari library estimatr yang langsung memberikan standard error yang telah dikoreksi.

Default lm_robust menggunakan metode koreksi HC2 dari McKinnon-White (1985). Standard error yang dihasilkan opsi robust di software Stata bisa direplikasi dengan menggunakan argumen se_type = "stata.

stargazer belum bisa melaporkan hasil dari lm_robust, namun texreg sudah mendukung. Gunakan argumen include.ci = FALSE untuk menampilkan standard error, bukan confidence interval yang merupakan tampilan default untuk lm_robust.

library(estimatr)
robust_default = lm_robust(morbidity_rate ~ number_of_doctors + number_of_midwives, datakab)
robust_stata = lm_robust(morbidity_rate ~ number_of_doctors + number_of_midwives, datakab, 
                         se_type = 'stata')

list(reg_morbid_2, robust_default, robust_stata) %>% 
texreg::htmlreg(digits = 4, include.ci = FALSE) %>% 
IRdisplay::display_html()

Bisa dilihat koreksi heterokesdastisitas mengubah nilai standard error dibanding hasil OLS. Akan tetapi tidak nampak ada perbedaan antara metode koreksi bawaan default lm_robust dengan metode koreksi opsi robust Stata.