Rによるニューラルネットワーク―neuralnet関数を利用した例―

個人的なメモとしてneuralnet関数の使い方についてまとめてみます。

関数の基本的な使い方については以下の本が詳しいのではないかと思います。

www.shoeisha.co.jp

本記事ではニューラルネットワークを構築して回帰問題に取り組みます。

目次

データセット

UCI Machine Learning Repository: Concrete Compressive Strength Data Set

  • データセットの中身は以下のとおり。
    • Number of instances 1030
    • Number of Attributes 9
    • Attribute breakdown 8 quantitative input variables, and 1 quantitative output variable
      • input variables: Cement, Slag, Ash, Water, Superplasticizer, Coarse Aggregate, Fine Aggregate, Age
      • Output Variable: Concrete compressive strength

使用するパッケージ

> library(readr)
> library(neuralnet)
> library(caret)

データの読込み

> data <- read.csv("C:/Users/admin/data/Concrete_data.csv")
> str(data)

'data.frame':	1030 obs. of  9 variables:
 $ cement          : num  540 540 332 332 199 ...
 $ slag            : num  0 0 142 142 132 ...
 $ ash             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ water           : num  162 162 228 228 192 228 228 228 228 228 ...
 $ superplasticizer: num  2.5 2.5 0 0 0 0 0 0 0 0 ...
 $ coarseagg       : num  1040 1055 932 932 978 ...
 $ fineagg         : num  676 676 594 594 826 ...
 $ age             : num  28 28 270 365 360 90 365 28 28 28 ...
 $ strength        : num  80 61.9 40.3 41 44.3 ...

特徴量の正規化

どうやら機械学習では特徴量を正規化して構築したモデルに投入することがお作法らしいです。特徴量を正規化することで変数の値(桁数)の大きさが予測精度に及ぼす影響を緩和しているということなんでしょうか?とりあえず、ここではお作法どおり各特徴量が区間{0, 1}の値を有するように正規化(置換)します。

> norm <- function(x){
>     return((x-min(x)) / (max(x)-min(x)))
> }
> data.norm <- as.data.frame(lapply(data, norm))
> summary(data.norm$strength)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.2664  0.4001  0.4172  0.5457  1.0000 

# 正規化前の要約統計量
> summary(data$strength)
  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.33   23.71   34.44   35.82   46.14   82.60 

訓練データとテストデータの作成

Yeh(1998)と同様にデータの75%を訓練データに利用して、25%をテストデータに利用する。

smp.size <- floor(0.75*nrow(data.norm))
set.seed(1030)
train.ind <- sample(seq_len(nrow(data.norm)), smp.size)
train <- data.norm[train.ind,]
test <- data.norm[-train.ind,]

neuralnet関数によるモデルの構築

neuralnet関数の引数のうちhiddenは隠れ層におけるユニット数を指定します。デフォルトでは1となっており、ここではとりあえずデフォルトどおり1を設定します。act.fctは活性化関数を指定する引数です。デフォルトではロジスティックシグモイド関数となっており、これもとりあえずデフォルトにしておきます。learningrateは学習率でありデフォルトはNULLですが、ここでは0.01に設定します。thresholdはイテレーション(学習)を止める閾値であり、デフォルトの0.01を設定します。

n <- names(train)
formula.bpn <- as.formula(paste("strength~",paste(n[!n%in%c("strength")],collapse="+")))

bpn <- neuralnet(formula = formula.bpn, 
                  data = train,
                  hidden = c(1), # 隠れ層におけるユニット数
                  act.fct="logistic", # 活性化関数にロジスティックシグモイド関数
                  learningrate = 0.01, # 学習率
                  threshold = 0.01, # 学習を止める閾値
                  linear.output=TRUE
                )

推定結果

neuralnet関数ではplot関数で推定結果を図示することができます。引数のrep="best"によって誤差(2乗誤差)が最も小さなものをプロットします。

plot(bpn, rep="best") 

f:id:joure:20170515211158p:plain
上図の左側の〇はコンクリートの圧縮強度を予測するために使用した8つの特徴量を示しており、隠れ層へ向かう黒線に書かれている数値は各特徴量のユニットが隠れ層のユニットに及ぼす影響度を表すウエイトを示しています。次に青色の線は各ステップにおけるバイアス項の値を示しており、これは回帰モデルにおける定数項みたいなものらしいです。
上図の下部にあるErrorとは損失関数(2乗誤差)の値を意味しており、ここでは正解ラベルの値と予測値の差の2乗和を2で割り算した値が5.053034であることを意味しています。この損失関数はneuralnet関数の引数(err.fct)で指定することができ、デフォルトでは"sse"(2乗誤差)が指定されています。sseの具体的な計算式は以下のように計算されます。

> bpn$err.fct
structure(function (x, y) 
{
    1/2 * (y - x)^2
}, type = "sse")

err.fctには他にも交差エントロピーを指定することも可能です。Stepsは学習が収束するまでのステップ数であり、ここでは4063回となっています。
上図で示されているErrorやStepsなどの値は以下のコマンドから取得できます。

> bpn$result.matrix
	                      1
error	                      5.053033984221
reached.threshold	      0.008933137079
steps	                      4063.000000000000
Intercept.to.1layhid1	      -3.010942534867
cement.to.1layhid1	      4.695003963273
slag.to.1layhid1	      2.944619408533
ash.to.1layhid1	              1.093491879204
water.to.1layhid1	      -1.388814333100
superplasticizer.to.1layhid1  1.251426794179
coarseagg.to.1layhid1	      0.274430837516
fineagg.to.1layhid1	      0.473816307600
age.to.1layhid1	              10.011792413222
Intercept.to.strength	      -0.013639063787
1layhid.1.to.strength	      0.691225238214

次に、目的変数の予測に利用した個々の特徴量が予測に役立つものかを確認する手軽な手法として一般化ウエイトと特徴量の散布図を確認することがあげられます*1。この散布図は横軸に特徴量、縦軸に一般化ウエイトをプロットしており、散布図の縦軸(一般化ウエイト)の値がゼロ近辺に特徴量がプロットされているとき、その特徴量は目的変数の予測に役立っていないと判断できます。ちなみに、今回のデータによる散布図は以下のとおりです。

# 一般化ウエイトと特徴量の散布図を作成
par(mfrow=c(2,4))
gwplot(bpn,selected.covariate="cement")
gwplot(bpn,selected.covariate="slag")
gwplot(bpn,selected.covariate="ash")
gwplot(bpn,selected.covariate="water")
gwplot(bpn,selected.covariate="superplasticizer")
gwplot(bpn,selected.covariate="coarseagg")
gwplot(bpn,selected.covariate="fineagg")
gwplot(bpn,selected.covariate="age")

f:id:joure:20170515211240p:plain

上図をみると縦軸の値がゼロのところで横軸と水平に値が分布している図は見当たらないので、予測に利用した特徴量はこのまま含めてもよいと判断できそうです。
この一般化ウエイトの値を特徴量ごとに確認するには以下のようにします。

> print(head(bpn$generalized.weights[[1]]))
            [,1]         [,2]         [,3]          [,4]         [,5]
895 3.5893379669 2.2511662022 0.8359762738 -1.0617507576 0.9567177664
838 3.1450228972 1.9725000310 0.7324928850 -0.9303193164 0.8382880936
482 0.5892127611 0.3695433158 0.1372308468 -0.1742931709 0.1570513343
698 5.2561940936 3.2965874500 1.2241961033 -1.5548182178 1.4010088544
955 3.9371617592 2.4693148337 0.9169863209 -1.1646394179 1.0494282341
253 3.0937397406 1.9403361864 0.7205487598 -0.9151494074 0.8246188578
             [,6]          [,7]         [,8]
895 0.20980280999 0.36223331770  7.654031159
838 0.18383185071 0.31739337136  6.706557998
482 0.03444047178 0.05946291357  1.256458120
698 0.30723333964 0.53045119811 11.208494084
955 0.23013369264 0.39733543609  8.395742909
253 0.18083426441 0.31221791335  6.597199982

訓練データによる予測とRMSE

ここではパラメタを未調整の状態で構築したニューラルネットワークによる予測値を計算し、その予測精度の指標であるRMSEを計算してみます。

> pred <- compute(bpn, test[,1:8])
> pred.strength <- pred$net.result
> rmse.bpn <- sqrt(sum((test$strength - pred.strength)^2)/nrow(test))
> print(rmse.bpn)
[1] 0.1215388084

RMSEの値は0.1215388084だと分かりました。これをチューニング後のモデルのRMSEと比べて、チューニングによって予測精度が高まっているかを以降で確認します。

neuralnet関数におけるパラメタのチューニング

caretパッケージのtrain関数を使用します。neuralnet関数を使用する場合は、隠れ層の数とそのユニット数がチューニングするパラメタとなります。
隠れ層の数やユニット数の組合せごとにRMSEの値を計算し、最も小さい値を有する組合せで再度ニューラルネットワークを構築します。
ここでは隠れ層が1つの場合と2つの場合を考え、1つ目の隠れ層のユニットが1から3つの場合と2つ目の隠れ層のユニットが0から3つの場合の合計12通りの組合せからRMSEを計算してみます*2

model <- train(form=formula.bpn,    
               data=train,           
               method="neuralnet",
               tuneGrid = expand.grid(.layer1=c(1:3), .layer2=c(0:3), .layer3=c(0)),
               learningrate = 0.01,
               threshold = 0.01
              )
> (model)
Neural Network 

772 samples
  8 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 772, 772, 772, 772, 772, 772, ... 
Resampling results across tuning parameters:

  layer1  layer2  RMSE           Rsquared    
  1       0       0.11673050316  0.6861291084
  1       1       0.12088359121  0.6591349899
  1       2       0.11698356020  0.6850522478
  1       3       0.11656958967  0.6869186251
  2       0       0.08534486470  0.8324254050
  2       1       0.08344890188  0.8401588492
  2       2       0.08380815321  0.8387608912
  2       3       0.08371587352  0.8390426413
  3       0       0.08594628891  0.8324012544
  3       1       0.08061781486  0.8499824496
  3       2       0.08074193620  0.8498935825
  3       3       0.07930272754  0.8556659702

Tuning parameter 'layer3' was held constant at a value of 0
RMSE was used to select the optimal model using  the smallest value.
The final values used for the model were layer1 = 3, layer2 = 3 and layer3 = 0. 

> plot(model)

f:id:joure:20170515211534p:plain

上図をみると、隠れ層1のユニットが3つで、隠れ層2のユニットが3つの場合のRMSEの値が約0.079と12通りの組合せの中で最も小さいことが分かります。これは、パラメタのチューニングによってモデルの精度が高くなったことを意味します。したがって、チューニング後のパラメタを用いて再度モデルを構築してみたいと思います。

チューニング後のパラメタを利用してモデルを再構築

bpn2 <- neuralnet(formula = formula.bpn, 
                  data = train,
                  hidden = c(3,3), # 隠れ層1:ニューロン数2
                  act.fct="logistic",
                  learningrate = 0.01, # learning rate
                  threshold = 0.01,
                  linear.output=TRUE
                )


plot(bpn2, rep="best") 

f:id:joure:20170515211646p:plain

上図をみると、最初のモデルと比べて、2乗誤差で表現されるErrorの値が減少していることが分かります。

テストデータによる予測とRMSE

結果は以下のとおり。

pred2 <- compute(bpn2, test[,1:8])
pred2.strength <- pred2$net.result
rmse.bpn2 <- sqrt(sum((test$strength - pred2.strength)^2)/nrow(test))
print(rmse.bpn2)

[1] 0.09096336435

RMSEの値は0.09096336435という結果になりました。

*1:一般化ウエイトについてはGünther and Fritsch(2010)を参照してください。

*2:本当はより多くの組合せを試して見たかったのですが、PCのスペックが低いので今回の12通りが限界でした