Rで重回帰分析の変数選択を自動化してみた


# やったこと
Rを使って複数の非説明変数に対して重回帰モデルをつくり、p値に対して設定したクライテリアを満たす変数のみを使ったモデルを作成する関数を作成した。

# モチベーション
１セットの実験で、複数の評価結果が得られるときに、1つずつモデルを作るのが面倒だった。

# 実装

サンプルデータとして、`iris`を使用した。Speciesは、文字型データであるためダミー変数に変換して使用した。
まずは、すべての４つの変数をすべて使用して重回帰モデルを作成してみる。

```r script.R
iris2 <- 
  iris %>% 
  mutate(Species = as.numeric(as.factor(iris$Species)))

model <- lm(Species ~ ., iris2)

summary(model)

# Call:
#   lm(formula = Species ~ ., data = iris2)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.59215 -0.15368  0.01268  0.11089  0.55077 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   1.18650    0.20484   5.792 4.15e-08 ***
#   Sepal.Length -0.11191    0.05765  -1.941   0.0542 .  
# Sepal.Width  -0.04008    0.05969  -0.671   0.5030    
# Petal.Length  0.22865    0.05685   4.022 9.26e-05 ***
#   Petal.Width   0.60925    0.09446   6.450 1.56e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.2191 on 145 degrees of freedom
# Multiple R-squared:  0.9304,	Adjusted R-squared:  0.9285 
# F-statistic: 484.5 on 4 and 145 DF,  p-value: < 2.2e-16
```

Sepal.Width以外は、p値が0.1以下であり有意である結果となることが確認できたため、暫定的に変数選択のクライテリアを0.1として、クライテリアを満たす変数のみでモデルをつくる関数をつくってみた。

```r script.R
lm2 <- function(criterion = 0.1){
  model <- lm(Species ~ ., iris2)
  coeff <- summary(model)$coefficients[-1,]
  vars <- row.names(coeff)[which(coeff[,4] <= criterion)]
  f <- formula(paste0('Species ~', paste(vars, collapse = '+'))) 
  model2 <- lm(f, iris2)
  return(model2)
}

model2 <- lm2()

summary(model2)
# Call:
#   lm(formula = f, data = iris2)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.60753 -0.16188  0.01367  0.11217  0.54740 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   1.14469    0.19478   5.877 2.72e-08 ***
#   Sepal.Length -0.13624    0.04475  -3.044  0.00277 ** 
#   Petal.Length  0.25213    0.04473   5.637 8.67e-08 ***
#   Petal.Width   0.58689    0.08822   6.652 5.41e-10 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.2187 on 146 degrees of freedom
# Multiple R-squared:  0.9302,	Adjusted R-squared:  0.9287 
# F-statistic: 648.3 on 3 and 146 DF,  p-value: < 2.2e-16
```

`lm2()`関数の引数のcriterionに閾値を設定することで、p値が閾値以下のパラメータのみで重回帰モデルをつくることができた。わずかだが自由度調整済み決定係数の値も改善した。


# 追記
回帰問題に`iris`データセットは適切ではないというご意見を頂いたので`airquality`データセットで改めて実施したいと思います。

データを確認してみるとOzoneとSolar.Rに欠損値があったので、それらの行をとりぞきました。

``` r scpript.R
library(dplyr)

summary(airquality)

# Ozone           Solar.R           Wind             Temp           Month            Day      
# Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00   Min.   :5.000   Min.   : 1.0  
# 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00   1st Qu.:6.000   1st Qu.: 8.0  
# Median : 31.50   Median :205.0   Median : 9.700   Median :79.00   Median :7.000   Median :16.0  
# Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88   Mean   :6.993   Mean   :15.8  
# 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00   3rd Qu.:8.000   3rd Qu.:23.0  
# Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00   Max.   :9.000   Max.   :31.0  
# NA's   :37       NA's   :7    


airquality2 <-
  airquality %>% 
  filter(Ozone != 'NA' &Solar.R != 'NA')

summary(airquality2)

# Ozone          Solar.R           Wind            Temp           Month            Day       
# Min.   :  1.0   Min.   :  7.0   Min.   : 2.30   Min.   :57.00   Min.   :5.000   Min.   : 1.00  
# 1st Qu.: 18.0   1st Qu.:113.5   1st Qu.: 7.40   1st Qu.:71.00   1st Qu.:6.000   1st Qu.: 9.00  
# Median : 31.0   Median :207.0   Median : 9.70   Median :79.00   Median :7.000   Median :16.00  
# Mean   : 42.1   Mean   :184.8   Mean   : 9.94   Mean   :77.79   Mean   :7.216   Mean   :15.95  
# 3rd Qu.: 62.0   3rd Qu.:255.5   3rd Qu.:11.50   3rd Qu.:84.50   3rd Qu.:9.000   3rd Qu.:22.50  
# Max.   :168.0   Max.   :334.0   Max.   :20.70   Max.   :97.00   Max.   :9.000   Max.   :31.00
````

欠損値を取り除けたので重回帰分析をしてみます。

``` r script.R
model <- lm(Ozone ~ ., airquality2)

summary(model)

# Call:
#   lm(formula = Ozone ~ ., data = airquality2)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -37.014 -12.284  -3.302   8.454  95.348 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -64.11632   23.48249  -2.730  0.00742 ** 
#   Solar.R       0.05027    0.02342   2.147  0.03411 *  
#   Wind         -3.31844    0.64451  -5.149 1.23e-06 ***
#   Temp          1.89579    0.27389   6.922 3.66e-10 ***
#   Month        -3.03996    1.51346  -2.009  0.04714 *  
#   Day           0.27388    0.22967   1.192  0.23576    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 20.86 on 105 degrees of freedom
# Multiple R-squared:  0.6249,	Adjusted R-squared:  0.6071 
# F-statistic: 34.99 on 5 and 105 DF,  p-value: < 2.2e-16


lm3 <- function(criterion = 0.1){
  model <- lm(Ozone ~ ., airquality2)
  coeff <- summary(model)$coefficients[-1,]
  vars <- row.names(coeff)[which(coeff[,4] <= criterion)]
  f <- formula(paste0('Ozone ~', paste(vars, collapse = '+'))) 
  model2 <- lm(f, airquality2)
  return(model2)
}

model3 <- lm3()

summary(model3)

# #Call:
# lm(formula = f, data = airquality2)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -35.870 -13.968  -2.671   9.553  97.918 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -58.05384   22.97114  -2.527   0.0130 *  
#   Solar.R       0.04960    0.02346   2.114   0.0368 *  
#   Wind         -3.31651    0.64579  -5.136 1.29e-06 ***
#   Temp          1.87087    0.27363   6.837 5.34e-10 ***
#   Month        -2.99163    1.51592  -1.973   0.0510 .  
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 20.9 on 106 degrees of freedom
# Multiple R-squared:  0.6199,	Adjusted R-squared:  0.6055 
# F-statistic: 43.21 on 4 and 106 DF,  p-value: < 2.2e-16
```

あと、変数選択の方法としてp値を指標とした選択のほかに一般的な赤池情報量基準（AIC）を指標とした方がいいのではというご意見もあったので、それも試してみます。

今回のデータセットに対しては、`step()`を用いて変数選択した場合とp値を指標に選択した場合とで結果は同様になりました。もちろん、p値のクライテリアを厳しくすれば違いは生じると思います。


``` r script
model3_1 <- step(model)

# Start:  AIC=680.21
# Ozone ~ Solar.R + Wind + Temp + Month + Day
# 
# Df Sum of Sq   RSS    AIC
# - Day      1     618.7 46302 679.71
# <none>                 45683 680.21
# - Month    1    1755.3 47438 682.40
# - Solar.R  1    2005.1 47688 682.98
# - Wind     1   11533.9 57217 703.20
# - Temp     1   20845.0 66528 719.94
# 
# Step:  AIC=679.71
# Ozone ~ Solar.R + Wind + Temp + Month
# 
# Df Sum of Sq   RSS    AIC
# <none>                 46302 679.71
# - Month    1    1701.2 48003 681.71
# - Solar.R  1    1952.6 48254 682.29
# - Wind     1   11520.5 57822 702.37
# - Temp     1   20419.5 66721 718.26

summary(model3_1)

# Call:
#   lm(formula = Ozone ~ Solar.R + Wind + Temp + Month, data = airquality2)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -35.870 -13.968  -2.671   9.553  97.918 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -58.05384   22.97114  -2.527   0.0130 *  
#   Solar.R       0.04960    0.02346   2.114   0.0368 *  
#   Wind         -3.31651    0.64579  -5.136 1.29e-06 ***
#   Temp          1.87087    0.27363   6.837 5.34e-10 ***
#   Month        -2.99163    1.51592  -1.973   0.0510 .  
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 20.9 on 106 degrees of freedom
# Multiple R-squared:  0.6199,	Adjusted R-squared:  0.6055 
# F-statistic: 43.21 on 4 and 106 DF,  p-value: < 2.2e-16
```