井出草平の研究ノート

トップ > R > MASSパッケージでロバスト推定のM推定、MM推定、P値を計算する[R]

MASSパッケージでロバスト推定のM推定、MM推定、P値を計算する[R]

R 回帰分析 計量

以前に書いたロバスト推定の方法。こちらにP値の出し方を書いていなかったので、補足。

ides.hatenablog.com

通常の重回帰分析

mod <- lm(wage ~ gender + age + education, data = CPS1985)
summary(mod)

結果：

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.84275    1.24425  -3.892 0.000112 ***
genderfemale -2.33542    0.38807  -6.018 3.30e-09 ***
age           0.11310    0.01669   6.775 3.32e-11 ***
education     0.82744    0.07462  11.089  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.454 on 530 degrees of freedom
Multiple R-squared:  0.253, Adjusted R-squared:  0.2488 
F-statistic: 59.85 on 3 and 530 DF,  p-value: < 2.2e-16

M推定

Huber(1981)によって提唱された方法。

      1. Huber (1981) Robust Statistics. Wiley.
library(AER)
data("CPS1985")

library(MASS)
mod.robust <- rlm(wage ~ gender + age + education, data = CPS1985)
summary(mod.robust)
Coefficients:
             Value   Std. Error t value
(Intercept)  -4.6060  1.0236    -4.4999
genderfemale -2.3399  0.3192    -7.3297
age           0.1171  0.0137     8.5281
education     0.7701  0.0614    12.5448

Residual standard error: 3.543 on 530 degrees of freedom

このようにt値までしか出力してくれない。

t 値から p-value を計算する

これは t 分布の累積分布関数であるpt() 関数を使用する方法がある。

# モデルの要約を取得
summary_mod <- summary(mod.robust)

# 自由度の計算
df <- nrow(CPS1985) - length(coef(mod.robust))

# 各係数の t 値から p-value を計算
p_values <- 2 * pt(-abs(summary_mod$coefficients[, "t value"]), df)

# p-value をsummary_modに追加
summary_mod$coefficients <- cbind(summary_mod$coefficients, "p value" = p_values)

# 結果を表示
print(summary_mod$coefficients)

結果。

                  Value Std. Error   t value      p value
(Intercept)  -4.6060026 1.02357607 -4.499912 8.361872e-06
genderfemale -2.3399281 0.31923960 -7.329693 8.668241e-13
age           0.1171214 0.01373363  8.528072 1.563109e-16
education     0.7700551 0.06138443 12.544797 8.596133e-32

MM推定

ちなみにMASSパッケージのrlm関数ではM推定とのほかにMM推定が実行できる。デフォルトはM推定。

mod.robust02 <- rlm(wage ~ gender + age + education, 
                  data = CPS1985, method = "MM")
summary(mod.robust02)

MM推定の結果は下記のようになる

Coefficients:
             Value   Std. Error t value
(Intercept)  -4.2019  1.0418    -4.0332
genderfemale -2.2212  0.3249    -6.8358
age           0.1162  0.0140     8.3141
education     0.7244  0.0625    11.5947

Residual standard error: 3.438 on 530 degrees of freedom

参考：

www.rdocumentation.org

ides.hatenablog.com