Mereplikasi opsi "kuat" Stata di R

39

Saya telah mencoba untuk mereplikasi hasil opsi Stata robustdi R. Saya telah menggunakan rlmperintah dari paket MASS dan juga perintah lmrobdari paket "robustbase". Dalam kedua kasus hasilnya sangat berbeda dari opsi "kuat" di Stata. Adakah yang bisa menyarankan sesuatu dalam konteks ini?

Berikut adalah hasil yang saya peroleh ketika saya menjalankan opsi robust di Stata:

. reg yb7 buildsqb7 no_bed no_bath rain_harv swim_pl pr_terrace, robust

Linear regression                                      Number of obs =    4451
                                                       F(  6,  4444) =  101.12
                                                       Prob > F      =  0.0000
                                                       R-squared     =  0.3682
                                                       Root MSE      =   .5721

------------------------------------------------------------------------------
             |               Robust
         yb7 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   buildsqb7 |   .0046285   .0026486     1.75   0.081    -.0005639     .009821
      no_bed |   .3633841   .0684804     5.31   0.000     .2291284    .4976398
     no_bath |   .0832654   .0706737     1.18   0.239    -.0552904    .2218211
   rain_harv |   .3337906   .0395113     8.45   0.000     .2563289    .4112524
     swim_pl |   .1627587   .0601765     2.70   0.007     .0447829    .2807346
  pr_terrace |   .0032754   .0178881     0.18   0.855    -.0317941    .0383449
       _cons |   13.68136   .0827174   165.40   0.000     13.51919    13.84353

Dan ini yang saya dapatkan di R dengan opsi lmrob:

> modelb7<-lmrob(yb7~Buildsqb7+No_Bed+Rain_Harv+Swim_Pl+Gym+Pr_Terrace, data<-bang7)
> summary(modelb7)

Call:
lmrob(formula = yb7 ~ Buildsqb7 + No_Bed + Rain_Harv + Swim_Pl + Gym + Pr_Terrace, 
    data = data <- bang7)
 \--> method = "MM"
Residuals:
      Min        1Q    Median        3Q       Max 
-51.03802  -0.12240   0.02088   0.18199   8.96699 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.648261   0.055078 229.641   <2e-16 ***
Buildsqb7    0.060857   0.002050  29.693   <2e-16 ***
No_Bed       0.005629   0.019797   0.284   0.7762    
Rain_Harv    0.230816   0.018290  12.620   <2e-16 ***
Swim_Pl      0.065199   0.028121   2.319   0.0205 *  
Gym          0.023024   0.014655   1.571   0.1162    
Pr_Terrace   0.015045   0.013951   1.078   0.2809    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Robust residual standard error: 0.1678 
Multiple R-squared:  0.8062,    Adjusted R-squared:  0.8059 
pengguna56579
sumber
3
Selamat Datang di Cross Divalidasi! Saya membuat judul Anda sedikit lebih deskriptif dan menambahkan beberapa format. Secara umum, pertanyaan pemrograman tidak pada topik di sini, tapi saya pikir Anda adalah karena melibatkan beberapa masalah statistik. Berharap untuk melihat Anda sekitar ....
Matt Krause
3
Akan sangat membantu jika Anda setidaknya menempelkan kode yang digunakan untuk memperkirakan model di Stata dan R (bahkan lebih baik jika Anda memberikan contoh yang sepenuhnya dapat direproduksi). Ketika Anda mengatakan "hasil berbeda" - jika Anda memperkirakan model yang sama hanya kesalahan standar yang harus berbeda, bukan estimasi koefisien.
Andy W
oke ... Ini adalah hasil yang saya peroleh dengan opsi yang kuat di STATA:
user56579
5
Sepertinya lmrobtidak sama dengan reg y x, robust. Google "kesalahan standar konsisten heteroskedastisitas R". Anda akan mendapatkan halaman yang menunjukkan kepada Anda bagaimana menggunakan lmtestdan sandwichperpustakaan.
generic_user
3
Stata menggunakan faktor koreksi sampel kecil n / (nk). R biasanya melakukan sesuatu yang lain, jadi pastikan Anda menyesuaikan untuk itu.
Dimitriy V. Masterov

Jawaban:

47

Charles hampir ada di jawabannya, tetapi robustopsi regressperintah (dan perintah estimasi regresi lainnya) di Stata memungkinkan untuk menggunakan beberapa jenis heteroskedastisitas dan autokorelasi kuat varians-kovarians matriks penduga, seperti halnya coeftestfungsi dalam lmtestpaket, yang pada gilirannya tergantung pada masing-masing varians-kovarians matriks yang dihasilkan oleh vcovHCfungsi disandwich paket.

Namun, matriks varians-kovarians default yang digunakan oleh keduanya berbeda:
1. Matriks varians-kovarians default yang dikembalikan oleh vcocHCadalah apa yang disebut HC3dengan alasan yang dijelaskan dalam halaman manual vcovHC.
2. sandwichOpsi yang digunakan oleh Charles coeftestmenggunakan HC0matriks varian-kovarian yang kuat.
3. Untuk mereproduksi perilaku default Stata menggunakan robustopsi dalam panggilan yang harus regressAnda mintavcovHC untuk menggunakan HC1matriks varians-kovarians yang kuat.

Baca lebih lanjut tentang itu sini .

Contoh berikut yang menunjukkan semua poin yang dibuat di atas didasarkan pada contoh di sini .

library(foreign)
library(sandwich)
library(lmtest)

dfAPI = read.dta("http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2.dta")
lmAPI = lm(api00 ~ acs_k3 + acs_46 + full + enroll, data= dfAPI)
summary(lmAPI)                                  # non-robust

# check that "sandwich" returns HC0
coeftest(lmAPI, vcov = sandwich)                # robust; sandwich
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC0"))    # robust; HC0 

# check that the default robust var-cov matrix is HC3
coeftest(lmAPI, vcov = vcovHC(lmAPI))           # robust; HC3 
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC3"))    # robust; HC3 (default)

# reproduce the Stata default
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC1"))    # robust; HC1 (Stata default)

Baris kode terakhir di atas mereproduksi hasil dari Stata:

use http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2
regress api00 acs_k3 acs_46 full enroll, robust
tchakravarty
sumber
Tautan ke data sudah mati. Bisakah Anda memperbarui tautannya? Apakah file ini sama: fakultas.smu.edu/tfomby/eco5350/data/Examples/elemapi2.dta ?
vasili111
Bagaimana cara mereproduksi juga interval kepercayaan?
vasili111
29

Saya menemukan deskripsi di situs web berikut yang mereplikasi opsi '' robust '' Stata di R.

https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r

Mengikuti instruksi, yang perlu Anda lakukan adalah memuat fungsi ke dalam sesi R Anda dan kemudian mengatur parameter '' kuat '' di fungsi ringkasan Anda ke TRUE.

summary(lm.object, robust=TRUE)
Alex Rato
sumber
10

Pada April 2018 saya yakin Anda menginginkan estimatrpaket itu , yang menyediakan pengganti hampir untuk lm. Beberapa contoh menarik hampir dari dokumentasi:

library(estimatr)
library(car)

# HC1 robust standard errors
model <- lm_robust(GPA_year2 ~ gpa0 + ssp, data = alo_star_men,
                   se_type = "stata")
summary(model)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     se_type = "stata")
#> 
#> Standard error type:  HC1 
#> 
#> Coefficients:
#>             Estimate Std. Error  Pr(>|t|) CI Lower CI Upper  DF
#> (Intercept) -3.60625    1.60084 0.0258665 -6.77180  -0.4407 137
#> gpa0         0.06814    0.02024 0.0009868  0.02812   0.1082 137
#> ssp          0.31917    0.18202 0.0817589 -0.04077   0.6791 137
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

# HC1 cluster robust standard errors
model2 <- lm_robust(GPA_year2 ~ gpa0 + ssp, cluster = ssp,
                   data = alo_star_men, se_type = "stata")
summary(model2)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     clusters = ssp, se_type = "stata")
#> 
#> Standard error type:  stata 
#> 
#> Coefficients:
#>             Estimate Std. Error Pr(>|t|) CI Lower CI Upper DF
#> (Intercept) -3.60625   1.433195 0.240821 -21.8167  14.6042  1
#> gpa0         0.06814   0.018122 0.165482  -0.1621   0.2984  1
#> ssp          0.31917   0.004768 0.009509   0.2586   0.3798  1
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

The carpaket kemudian memudahkan untuk melakukan tes hipotesis omnibus untuk model ini:

linearHypothesis(model, c("gpa0 = ssp"))
#> Linear hypothesis test
#> 
#> Hypothesis:
#> gpa0 - ssp = 0
#> 
#> Model 1: restricted model
#> Model 2: GPA_year2 ~ gpa0 + ssp
#> 
#>   Res.Df Df  Chisq Pr(>Chisq)
#> 1    138                     
#> 2    137  1 1.8859     0.1697
alexpghayes
sumber
4

Saya akan mengedit pertanyaan. Anda mengacaukan regresi yang kuat dengan perintah Stata yang kuat. Tampaknya tidak ada manfaatnya menimbulkan kebingungan ini.

Saya pikir ada beberapa pendekatan. Saya belum melihat semuanya dan tidak yakin mana yang terbaik:

Paket sandwich:

library(sandwich)    
coeftest(model, vcov=sandwich)

Tetapi ini tidak memberi saya jawaban yang sama yang saya dapatkan dari Stata karena beberapa alasan. Saya belum pernah mencoba mencari tahu mengapa - tetapi di atas dalam komentar ada jawaban yang disarankan - Saya hanya tidak menggunakan paket ini.

Paket rms:

Saya merasa ini agak merepotkan, tetapi biasanya mendapatkan jawaban yang baik dengan sedikit usaha. Dan itu yang paling bermanfaat buat saya.

model = ols(a~b, x=TRUE)    
robcov(model)

Anda dapat mengkodekannya dari awal

Lihat posting blog ini ( http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/ ). Sepertinya pilihan yang paling menyakitkan, tetapi sangat mudah dan opsi ini sering bekerja paling baik.

charles
sumber
4
Charles benar pada poin utama, tetapi untuk membuat eksplisit apa yang tersirat di tempat lain perhatikan bahwa Stata tidak memiliki robustperintah! (Ada perintah programmer _robust, tidak langsung relevan di sini.) Sebaliknya, untuk mendapatkan kesalahan standar yang kuat (Huber-Eicker-White-sandwich), pendekatan modern di Stata adalah menentukan vce(robust)sebagai opsi. Pendekatan lama untuk menentukan robustopsi masih berfungsi. Lebih luas lagi, kebingungan yang disebabkan oleh perbedaan antara regresi yang kuat (dll.) Dan SI yang "kuat" sangat disayangkan.
Nick Cox
Hei. Terima kasih banyak. Kode berfungsi dan memang memberikan hasil yang dilakukan Stata. Hanya sebuah pertanyaan. Saya mengerti bahwa regresi yang kuat berbeda dari kesalahan standar yang kuat, dan bahwa regresi yang kuat digunakan ketika data Anda mengandung outlier. Tetapi itu juga memecahkan masalah heteroskedastisitas. Adakah yang bisa memberi tahu saya apakah estimasi jenis MM yang disediakan oleh perintah "lmrob" dari paket "robustbase" dapat digunakan sebagai solusi untuk masalah pencilan dan heteroskedastisitas secara simultan?
user56579
@ user56579 Dugaan saya adalah Anda ingin mengajukan pertanyaan terpisah tentang ini.
tchakravarty