Cara menghitung kesalahan standar yang kuat di r

Salah satu asumsi regresi linier adalah bahwa residu model tersebar secara merata di setiap level variabel prediktor.

Jika asumsi ini tidak terpenuhi maka dikatakan terjadi heteroskedastisitas dalam model regresi.

Jika hal ini terjadi, kesalahan standar koefisien regresi model menjadi tidak dapat diandalkan.

Untuk memperhitungkan hal ini, kita dapat menghitung kesalahan standar yang kuat , yang “kuat” terhadap heteroskedastisitas dan dapat memberi kita gambaran yang lebih baik tentang nilai kesalahan standar yang sebenarnya untuk koefisien regresi.

Contoh berikut menunjukkan cara menghitung kesalahan standar yang kuat untuk model regresi di R.

Contoh: menghitung kesalahan standar yang kuat di R

Misalkan kita memiliki kerangka data berikut di R yang berisi informasi tentang jam belajar dan nilai ujian yang diperoleh 20 siswa dalam satu kelas:

 #create data frame
df <- data. frame (hours=c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4,
                         4, 5, 5, 5, 6, 6, 7, 7, 8),
                 score=c(67, 68, 74, 70, 71, 75, 80, 70, 84, 72,
                         88, 75, 95, 75, 99, 78, 99, 65, 96, 70))

#view head of data frame
head(df)

  hours score
1 1 67
2 1 68
3 1 74
4 1 70
5 2 71
6 2 75

Kita dapat menggunakan fungsi lm() untuk menyesuaikan model regresi di R yang menggunakan jam sebagai variabel prediktor dan skor sebagai variabel respons:

 #fit regression model
fit <- lm(score ~ hours, data=df)

#view summary of model
summary(fit)

Call:
lm(formula = score ~ hours, data = df)

Residuals:
    Min 1Q Median 3Q Max 
-19,775 -5,298 -3,521 7,520 18,116 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 71.158 4.708 15.11 1.14e-11 ***
hours 1.945 1.075 1.81 0.087 .  
---
Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.48 on 18 degrees of freedom
Multiple R-squared: 0.154, Adjusted R-squared: 0.107 
F-statistic: 3.278 on 1 and 18 DF, p-value: 0.08696

Cara termudah untuk memeriksa secara visual apakah heteroskedastisitas merupakan masalah dalam model regresi adalah dengan membuat plot sisa:

 #create residual vs. fitted plot
plot(fitted(fit), reside(fit))

#add a horizontal line at y=0 
abline(0,0) 

Sumbu x menunjukkan nilai yang sesuai dari variabel respons dan sumbu y menunjukkan residu yang sesuai.

Dari grafik kita dapat melihat bahwa varians dari residu meningkat seiring dengan meningkatnya nilai yang dipasang.

Hal ini menunjukkan bahwa heteroskedastisitas kemungkinan besar menjadi masalah dalam model regresi dan kesalahan standar ringkasan model tidak dapat diandalkan.

Untuk menghitung kesalahan standar yang kuat, kita dapat menggunakan fungsi coeftest() dari paket lmtest dan fungsi vcovHC() dari paket sandwich sebagai berikut:

 library (lmtest)
library (sandwich)

#calculate robust standard errors for model coefficients
coeftest(fit, vcov = vcovHC(fit, type = ' HC0 '))

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 71.1576 3.3072 21.5160 2.719e-14 ***
hours 1.9454 1.2072 1.6115 0.1245    
---
Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perhatikan bahwa kesalahan standar untuk variabel prediktor jam kerja meningkat dari 1,075 pada ringkasan model sebelumnya menjadi 1,2072 dalam ringkasan model ini.

Karena terdapat heteroskedastisitas dalam model regresi asli, estimasi kesalahan standar ini lebih dapat diandalkan dan sebaiknya digunakan saat menghitung interval kepercayaan untuk variabel prediktor jam .

Catatan : Jenis perkiraan yang paling umum untuk dihitung dalam fungsi vcovHC() adalah ‘HC0’, namun Anda dapat merujuk ke dokumentasi untuk menemukan jenis perkiraan lainnya.

Sumber daya tambahan

Tutorial berikut menjelaskan cara melakukan tugas umum lainnya di R:

Bagaimana melakukan uji White untuk heteroskedastisitas di R
Bagaimana menginterpretasikan keluaran regresi linier di R
Cara membuat plot sisa di R