やったこと
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値のクライテリアを厳しくすれば違いは生じると思います。
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