Bagaimana cara mendapatkan tabel ANOVA dengan kesalahan standar yang kuat?

10

Saya menjalankan regresi OLS yang dikumpulkan menggunakan paket plm di R. Meskipun, pertanyaan saya lebih lanjut tentang statistik dasar, jadi saya coba posting di sini dulu;)

Karena hasil regresi saya menghasilkan residu heteroskedastik, saya ingin mencoba menggunakan kesalahan standar kuat heteroskedastisitas. Sebagai hasil dari coeftest(mod, vcov.=vcovHC(mod, type="HC0"))saya mendapatkan tabel yang berisi perkiraan, kesalahan standar, nilai-t dan nilai-p untuk setiap variabel independen, yang pada dasarnya adalah hasil regresi "kuat" saya.

Untuk membahas pentingnya variabel yang berbeda, saya ingin merencanakan pembagian varian yang dijelaskan oleh masing-masing variabel independen, jadi saya perlu jumlah kuadrat masing-masing. Namun, menggunakan fungsi aov(), saya tidak tahu bagaimana cara memberitahu R untuk menggunakan kesalahan standar yang kuat.

Sekarang pertanyaan saya adalah: Bagaimana cara saya mendapatkan tabel / jumlah kuadrat ANOVA yang mengacu pada kesalahan standar yang kuat? Apakah mungkin untuk menghitungnya berdasarkan tabel ANOVA dari regresi dengan kesalahan standar normal?

Edit:

Dengan kata lain dan mengabaikan masalah-R saya:

Jika R tidak terpengaruh dengan menggunakan kesalahan standar yang kuat, apakah kontribusi masing-masing untuk menjelaskan varians oleh variabel penjelas yang berbeda tidak akan berubah?2

Edit:

Dalam R, apakah aov(mod)benar - benar memberikan tabel ANOVA yang benar untuk panelmodel (plm)?

Aki
sumber

Jawaban:

12

ANOVA dalam model regresi linier setara dengan uji Wald (dan uji rasio kemungkinan) dari model bersarang yang sesuai. Jadi, ketika Anda ingin melakukan tes yang sesuai menggunakan kesalahan standar heteroskedasticity-konsisten (HC), ini tidak dapat diperoleh dari dekomposisi jumlah kuadrat tetapi Anda dapat melakukan tes Wald menggunakan estimasi kovarians HC. Ide ini digunakan dalam kedua Anova()dan linearHypothesis()dari carpaket dan coeftest()dan waldtest()dari lmtestpaket. Tiga yang terakhir juga dapat digunakan dengan plmobjek.

Contoh sederhana (meskipun tidak terlalu menarik / bermakna) adalah sebagai berikut. Kami menggunakan model standar dari ?plmhalaman manual dan ingin melakukan uji Wald untuk signifikansi keduanya log(pcap)dan unemp. Kami membutuhkan paket-paket ini:

library("plm")
library("sandwich")
library("car")
library("lmtest")

Model (di bawah alternatif) adalah:

data("Produc", package = "plm")
mod <- plm(log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp,
  data = Produc, index = c("state", "year"))

Pertama, mari kita lihat tes Wald marginal dengan kesalahan standar HC untuk semua koefisien individu:

coeftest(mod, vcov = vcovHC)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
log(pc)    0.2920069  0.0617425  4.7294 2.681e-06 ***
log(emp)   0.7681595  0.0816652  9.4062 < 2.2e-16 ***
log(pcap) -0.0261497  0.0603262 -0.4335   0.66480    
unemp     -0.0052977  0.0024958 -2.1226   0.03411 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dan kemudian kami melakukan tes Wald untuk keduanya log(pcap)dan unemp:

linearHypothesis(mod, c("log(pcap)", "unemp"), vcov = vcovHC)

Linear hypothesis test

Hypothesis:
log(pcap) = 0
unemp = 0

Model 1: restricted model
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp

Note: Coefficient covariance matrix supplied.

  Res.Df Df  Chisq Pr(>Chisq)  
1    766                       
2    764  2 7.2934    0.02608 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Atau, kita juga dapat menyesuaikan model di bawah hipotesis nol ( mod0katakanlah) tanpa dua koefisien dan kemudian panggil waldtest():

mod0 <- plm(log(gsp) ~ log(pc) + log(emp),
  data = Produc, index = c("state", "year"))
waldtest(mod0, mod, vcov = vcovHC)

Wald test

Model 1: log(gsp) ~ log(pc) + log(emp)
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp
  Res.Df Df  Chisq Pr(>Chisq)  
1    766                       
2    764  2 7.2934    0.02608 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Statistik uji dan nilai-p yang dihitung oleh linearHypothesis()dan waldtest()persis sama. Hanya antarmuka dan format keluaran agak berbeda. Dalam beberapa kasus satu atau yang lain lebih sederhana atau lebih intuitif.

Catatan: Jika Anda memberikan taksiran matriks kovarians (yaitu, matriks seperti vocvHC(mod)) alih-alih penduga matriks kovarians (yaitu, fungsi seperti vocvHC), pastikan Anda menyediakan taksiran matriks kovarians HC dari model di bawah alternatif, yaitu, model tidak terbatas.

Achim Zeileis
sumber
1
Jika saya memahaminya dengan benar, Wald-test menunjukkan apakah memasukkan variabel-variabel tertentu signifikan atau tidak. Tetapi apakah ada ukuran seberapa banyak mereka benar-benar memperbaiki model, seperti misalnya jumlah kuadrat?
Aki
Bagaimana saya bisa mengimplementasikan kesalahan standar HAC? Saya mencoba coeftest (mod, vcov = vcovHAC) tetapi mendapat pesan kesalahan yang mengatakan "tidak ada metode yang berlaku untuk 'estfun' diterapkan pada objek kelas" c ('plm', 'panelmodel') ". Sepertinya saya perlu menghitung bobot atau fungsi estimasi terlebih dahulu. Metode apa yang akan Anda rekomendasikan?
Aki
Sementara plmpaket memiliki metode untuk vcovHC()generik dari sandwichpaket, itu tidak menyediakan metode untuk vcovHAC(). Sebagai gantinya plmdikirimkan dengan fungsi khusus untuk menghitung matriks kovarian dalam model panel yang berpotensi menyertakan korelasi serial juga. Lihat vcovNW()atau vcovSCC()dalam paket plm.
Achim Zeileis
Terima kasih! Sejauh yang saya mengerti kedua fungsi berhubungan dengan autocorrelation-robust SE. Apakah ada fungsi yang dapat digunakan untuk SE heteroscedasticity- dan autocorrelation-robust?
Aki
Semua tiga fungsi ( vcovHAC, vcovNW, vcovSCC) adalah _H_eteroskedasticity dan _A_utocorrelation _C_onsistent ... yang ini apa berdiri HAC untuk.
Achim Zeileis
2

Ini dapat dilakukan dengan Anovafungsi dalam carpaket. Lihat jawaban yang luar biasa ini untuk lebih detail dan ulasan teknik lain untuk berurusan dengan heteroskedastisitas di ANOVA.

shadowtalker
sumber
Terima kasih. Masalahnya adalah, bahwa Anova () tampaknya tidak bekerja dengan model tipe plm (panelmodel).
Aki
@ Aki jika saya tidak salah OLS yang dikumpulkan setara dengan OLS biasa, berdasarkan apa yang tertulis di sketsa: cran.r-project.org/web/packages/plm/vignettes/plm.pdf
shadowtalker
@ Aki namun sepertinya Anda mungkin tertarik pada model ANOVA yang lebih kaya. Ada beberapa contoh R di sini: stats.stackexchange.com/questions/3874/…
shadowtalker