R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Robusti regressio ja robustit keskivirheet

Robustilla regressiolla tarkoitetaan yleensä menetelmää, joka ei ole esimerkiksi tavanomaisen lineaarisen regression tapaan yhtä herkkä poikkeaville havainnoille. R:ssä robustin regression voi sovittaa esimerkiksi MASS-paketin funktiolla rlm(). Funktio sovittaa lineaarisen regression käyttäen M-estimaattoria, joka on selitetty hyvin Cross Validated -palstan keskustelussa. M-estimaattori on robusti vastemuuttujassa oleville poikkeaville havainnoille, muttei välttämättä selittävien muuttujien poikkeaville havainnoille. Jos mallin residuaaleissa on lisäksi heteroskedastisuutta, ei rlm() ratkaise ongelmaa, vaan tuolloin tulee käyttää robusteja keskivirhe-estimaatteja, jotka voidaan tuottaa paketilla sandwich. Lineaarinen regressio robusteilla keskivirheillä ei siis ole sama asia kuin robusti lineaarinen regressio.

Otetaanpa esimerkki.

Robusti regressio

Sovitetaan ensin tavanomainen lineaarinen regressio ja tutkitaan sen tulokset:

library(MASS)
fit1 <- lm(stack.loss ~ ., stackloss)
summary(fit1)
 
Call:
lm(formula = stack.loss ~ ., data = stackloss)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-7.2377 -1.7117 -0.4551  2.3614  5.6978 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.9197    11.8960  -3.356  0.00375 ** 
Air.Flow      0.7156     0.1349   5.307  5.8e-05 ***
Water.Temp    1.2953     0.3680   3.520  0.00263 ** 
Acid.Conc.   -0.1521     0.1563  -0.973  0.34405    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared:  0.9136,    Adjusted R-squared:  0.8983 
F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

Robusti regressio voidaan sovittaa seuraavasti:

library(MASS)
fit2 <- rlm(stack.loss ~ ., stackloss)
summary(fit2)
 
Call: rlm(formula = stack.loss ~ ., data = stackloss)
Residuals:
     Min       1Q   Median       3Q      Max 
-8.91753 -1.73127  0.06187  1.54306  6.50163 
 
Coefficients:
            Value    Std. Error t value 
(Intercept) -41.0265   9.8073    -4.1832
Air.Flow      0.8294   0.1112     7.4597
Water.Temp    0.9261   0.3034     3.0524
Acid.Conc.   -0.1278   0.1289    -0.9922
 
Residual standard error: 2.441 on 17 degrees of freedom

Robustin regression antamat tulokset eroavat tavanomaisen lineaarisen regression antamista tuloksista, sillä menetelmät optimoivat eri kriteeriä, joten tämä on odotettavissakin.

Robusti regressio ei anna muuttujille myöskään p-arvoja, mutta ne voidaan estimoida esimerkiksi normaaliapproksimaatiotilla. P-arvot voidaan estimoida myös robustilla F-testillä (Wald test), jonka tarjoaa paketti sfsmisc:

library(sfsmisc)
f.robftest(fit2, var = "Air.Flow")
 
        robust F-test (as if non-random weights)
 
data:  from rlm(formula = stack.loss ~ ., data = stackloss)
F = 50.879, p-value = 1.677e-06
alternative hypothesis: true Air.Flow is not equal to 0

Robustit keskivirheet

Jos mallin residuaaleissa on heteroskedastisuutta, voidaan keskivirheiden arvioinnissa käyttää robusteja estimaattoreita, kuten sandwich-estimaattoria.

Tavanomaiselle lineaariselle regressiolle tämä voidaan tehdä seuraavasti:

library(sandwich)
library(lmtest)
 
# Sandwich
coeftest(fit1, vcov=sandwich)
 
t test of coefficients:
 
              Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -39.919674   6.411649 -6.2261 9.215e-06 ***
Air.Flow      0.715640   0.158944  4.5025 0.0003141 ***
Water.Temp    1.295286   0.446528  2.9008 0.0099458 ** 
Acid.Conc.   -0.152123   0.086429 -1.7601 0.0963744 .  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
# Heteroscedasticity consistent (HC1)
coeftest(fit1, vcov=vcovHC(fit1, type="HC1"))
 
t test of coefficients:
 
              Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -39.919674   7.126150 -5.6019 3.175e-05 ***
Air.Flow      0.715640   0.176657  4.0510 0.0008302 ***
Water.Temp    1.295286   0.496288  2.6099 0.0182994 *  
Acid.Conc.   -0.152123   0.096061 -1.5836 0.1317090    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Tags: , ,