はじめに

遺伝的アルゴリズム（GA）に基づいて説明変数選択しPLS回帰を行うGAPLS法をコードに書き起こしたものです。実施環境は以下の通りです。

R version 4.0.0 (2020-04-24)
使用ライブラリ
- baseline 1.3-1
- prospectr 0.2.1
- GA 3.2
- pls 2.7-2

この記事が目的とすること

GAPLSをコードに書き起こすことで私自身が理解を深めること
コードを自作することでライブラリに依存しなくてもよくなるようにすること
GAPLSのパラメータ最適化をスペクトルの前処理と合わせて実施すること
- ベースライン補正（Asymmetry Least Squares）
- スムージング（Savitzky Golay）

この記事が目的としないこと

GAPLS手法そのものの説明¹
新しい知見・手法の提供

使用データ

pls関数にあるサンプルデータ「yarn」を使用します。PET繊維の密度をを268波長のNIRスペクトルから推定します。トレーニングセット21点、テストセット7点の合計28点のデータセットです。

実装

では早速実装してみます。

初期設定

setting

# library
library(pls)
library(baseline)
library(prospectr)
library(GA)

# data
data(yarn) # sample data
X0 <- yarn$NIR
Y0 <- as.matrix(yarn$density)

trainSel <- yarn$train

Xtrain <- X0[trainSel,]
Xtest  <- X0[!trainSel,]
Ytrain <- as.matrix(Y0[trainSel])
Ytest  <- as.matrix(Y0[!trainSel])

trainSelはTRUE/FALSEで区別されています。前21点がトレーニングセット、下側7点がテストデータセットでした。
データセットX,Yは、大文字は関数の外、小文字は関数の内側で区別させています（自己流ですが...）

オリジナルデータがどんな形をしているか確認しておきます。

plot code

# 確認
dim(X0); dim(Y0)
ycol <- heat.colors(100)[ round((Y0 - min(Y0))/diff(range(Y0))*99,0) +1 ]
matplot(t(X0), type="l", col=ycol)

密度が大きいほど赤く、小さいほど黄色で示しています。具体的な波長はわかりませんが、Indexで40~50あたりのピークの裾・75あたりの谷・90~140のなだらかなピークあたりに繊維密度とのきれいな相関が見られそうです。
データセットの詳細な説明は有料²で確証はできませんがスペクトルの形から推測するに透過スペクトル強度かと思われます。回帰の際には吸収スペクトル強度に変換してから前処理をかけていくことにします。

前処理

前処理は三段階で行います。「（吸収スペクトル強度変換 ⇒）ベースライン処理 ⇒ スムージング処理 ⇒ 選択波長の抽出」の順です。
結果確認のところでも同様の操作を行うのでfPreprocessXという関数で定義しておくことにします。

preprocess

fPreprocessX <- function(p,x, selX=NULL){
  #function preprocess X
  #p[1]~p[2]:baseline.als()
  #p[3]~p[4]:savitzky-Golay
  #p[5]~:select explanatory variable

  p[3:4] <- floor(p[3:4])
  param <- list(bl.lambda=p[1], bl.p=p[2],
               sg.m=p[3], sg.p=sum(p[3:4]), sg.w=9)
  
  xabs <- max(x)-x
  xb   <- baseline.als(xabs, lambda=param$bl.lambda, p=param$bl.p)
  xbsg <- savitzkyGolay(xb$corrected, m=param$sg.m, p=param$sg.p, w=param$sg.w)
  
  xsel <- NULL
  if(!is.null(selX)){
    selXcol <- unique(floor(p[-(1:4)]))
    selXcol <- selXcol[selXcol<ncol(xbsg)]
    xsel    <- xbsg[,selXcol]
    param   <- c(param, selXcol=NULL)
    param$selXcol <- selXcol
  }
  
  result <- list(xb = xb, xbsg = xbsg, xsel = xsel, param = param)
  return(result)
}

pは前処理のハイパーパラメータを格納したベクトルです（これをGAで最適化します）。
7行目：最適化するハイパラには整数値と実数値が混ざっています（ベースライン処理は実数、ほかは全て整数）。なので実数値で生成したあとに必要なものを整数値に変換させています³。
8行目：使用するハイパラを変換しています。Savitzky-Golay法には3つの入力（m：微分の次数、p：補間の次数、w：窓幅）があり、それぞれ制限があります。それを満たすための変換を行っています。窓幅は説明変数行列の列数を変えてしまい、後の変数選択に影響を与えてしまう恐れがあることから固定させています（今回のコードはえいやで9にしていますが根拠はありません）⁴。
- m≦p
- 「p≦w」かつ「wは奇数」
16~22行目：説明変数の選択を行っています。ランダムに生成された列番号ではエラーを吐くのでいくつかの処理をしています。
- 17行目：同じ列番号が生成されてしまうことも当然あるのでunique()で解消
- 18行目：スムージングの窓幅によって説明変数行列の列数は変わってくるので、それを解消するため。今回は9で固定しているので、GAの大元の探索範囲を制限すれば問題ありません

GAPLSで最適化する関数

GAPLSに最適化させる関数です。入力ベクトルpの値を変化させて、出力を最大化させます。今回は実装の簡便性を取って、pls関数のfitted.valuesでR2値が最大となる因子数を計算してそのときのR2値を出力としました⁵。R2を計算する関数、RMSEを計算する関数をついでに定義しておきます。

fitting_function

fR2 = function(obs,pred){
  f = (!is.na(obs) & !is.na(pred))
  res = sum((obs-pred)^2)
  Var = sum((obs-mean(obs))^2)
  1-res/Var
}

fRMSE = function(obs, pred){
  f = (!is.na(obs) & !is.na(pred))
  N = length(obs)
  RMSE = sqrt(sum((obs-pred)^2)/N)
  return(RMSE)
}

fGAPLS <- function(p){
  #function optimesed by GA
  #p[1]~p[2]:baseline.als()
  #p[3]~p[4]:savitzky-Golay
  #p[5]~:select explanatory variable

  xpre <- fPreprocessX(x=Xtrain, p=p, selX=T)
  pls  <- plsr(Ytrain ~ xpre$xsel, scale=FALSE, validation="LOO")
  
  # val R2(改良の余地大いにアリ)
  R2.comp <- apply(pls$fitted.values[,1,],2,fR2,obs=Ytrain)
  ncomp <- which.min(R2.comp)
  
  returnGA <- R2.comp[ncomp]
  return(returnGA)  
}

7行目(fGAPLS関数内)
- 変数選択せずに前処理だけの最適化を実施する場合はxpre$xselをxpre$xbsgとすればOKです。
- スペクトルなのでscale=FALSEとしています。
10~11行目(fGAPLS関数内)
- R2値を各因子数について計算しています。最もR2値の小さい因子数を自動的に持ってきています。PLSの因子数上限によっては大きすぎる値を選択してしまう可能性もあるので、「最初の極小値にする」「別の評価因子を使う」などの改良余地があります。

GAPLSの実行

GAPLSを実行します。最適化するパラメータの探索範囲を準備します。今回はGAの各パラメータ（突然変異の割合など）などはデフォルトパラメータのままです。

GAPLS

# GA input parameter
bl.lambda.range <- c(1,10)
bl.p.range      <- c(0.001,0.1)
sg.m.range      <- c(0,2.99)
sg.p.range      <- c(1,7)

nExVal          <- 4
ExVal.range     <- c(1, ncol(Xtrain)-1)

GAfitParRange <- cbind(bl.lambda.range, bl.p.range, sg.m.range, sg.p.range,
                       matrix(rep(ExVal.range,times=nExVal),2,nExVal))

fit.Param <- ga(type="real-valued", fGAPLS, maxiter=100, seed=1,
                lower=GAfitParRange[1,], upper=GAfitParRange[2,])

10行目：GAfitParRangeは1行目に探索範囲の下限、2行目に探索範囲の上限値をまとめています（※入力がベクトルであれば特にこの形である必要はないです）。ほかに最適化したいパラメータがある場合は入力pに付け足していけばOKです。
13行目：GA実行部分です。上で定義したfGAPLS()にパラメータの組合せを入力して、出力を最大化させるよう実行してくれます。

コード実装の根幹部分は以上です。

結果の確認

結果を確認していきます。最適化基準（今回はR2値）の世代推移を確認します。

result1

plot(fit.Param)

緑の点が各世代における適合度のベスト、青線が世代の平均、緑塗りつぶしが世代の中央値を示しています。ベストの値を見ると1世代目からほぼ1になっていることから、今回のケースではハイパラにあまりこだわらなくても良い（≠汎化性能がある）ことがわかりました。汎化性能を議論するには適切な交差検証方法、例えばグループK分割などの方法を取ることが良いと思います⁶。

続いて最適化したハイパラメータの確認をします。

result2

fit.Param@solution #今回はこれだとわかりにくいので↓で
Xfit <- fPreprocessX(p=fit.Param@solution[1,], X0, selX=T)
Xfit$param

ベースライン処理のハイパラメータはlambdaが6.8、pが0.05となりました。微分の次数は1、多項式の次数は4となっています。選択された波長はインデックスで165, 148, 143, 167でした。GAPLSの場合、最終的に欲しい情報の一つはこの選択の情報です。スペクトルのプロットとともに選択波長を確認してみます。

plot code

result3

matplot(t(Xfit$xbsg),col=ycol,type="l",main="ベースライン + ノイズ除去", xlab="Index", ylab="")
abline(v=Xfit$param$selXcol,lty=2,col="grey60")

![Rplot04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/619386/8c73ce62-e35f-fa7e-cdb9-666d13c2020c.png)

150付近の小さなピークを選択してきていました。意外にも、オリジナルデータの観察時に着目していた波長領域とは異なるところを引っ張ってきているようです。これが果たしてグローバルな最適解なのかは、GAの計算を複数セット回す・世代数などのGA初期設定を変える・選択波長の吸収帯がドメインから見て妥当性を考える、といった方法で確認していくのがよいと思います。

最後に、実測予測プロットを見ておきます。

result4

fitpls  <- plsr(Ytrain ~ Xfit$xbsg[trainSel,], scale=F)
plot(fitpls, "validation")

PRESSを確認すると因子数は1で十分なようです。

plot code

result5

ncomp <- 1
Yfit  <- fitpls$fitted.values[,1,ncomp]
Ypred <- predict(fitpls, newdata=Xfit$xsel[!trainSel,])[,1,ncomp]

# 実測予測プロット
{
  ranges <- range(Ytrain,Yfit,Ytest,Ypred)
  plot(Ytrain,Yfit,pch=5,
       xlab="obs", ylab="pred", xlim=ranges, ylim=ranges)
  points(Ytest,Ypred,pch=4,col=2)
  abline(a=0, b=1, lty=2)
  
  R2.train   <- round(fR2(obs=Ytrain, pred=Yfit),3)
  RMSE.train <- round(fRMSE(obs=Ytrain, pred=Yfit),3)
  RMSE.test  <- round(fRMSE(obs=Ytest,  pred=Ypred),3)
  legends    <- c(paste0("R2.train = ", R2.train),
                  paste0("RMSE.train = ",RMSE.train),
                  paste0("RMSE.test = ",RMSE.test))
  legend("topleft", legend=legends, box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
  legend("bottomright", legend=c("train","test"),col=1:2,pch=c(5,4),cex=1.5,
         box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
}

4波長だけからでも学習データ・テストデータともにかなりの精度で予測することができました。

まとめ

GAPLSを自作コードで実装しました。自分でひとつひとつアルゴリズムを確認しコードに落とし込むという作業は良い思考整理方法だと感じています。次回はGAWLSの実装について気が向いたらまとめます。

実装コード

Code All

GAPLS.R

# setting ----

# library
library(pls)
library(baseline)
library(prospectr)
library(GA)

# data
data(yarn) # sample data
X0 <- yarn$NIR
Y0 <- as.matrix(yarn$density)

trainSel <- yarn$train

Xtrain <- X0[trainSel,]
Xtest  <- X0[!trainSel,]
Ytrain <- as.matrix(Y0[trainSel])
Ytest  <- as.matrix(Y0[!trainSel])


# 確認
dim(X0); dim(Y0)
ycol <- heat.colors(100)[ round((Y0 - min(Y0))/diff(range(Y0))*99,0) +1 ]
matplot(t(X0), type="l", col=ycol, xlab="Index", ylab="Intensity")

# function 

fR2 = function(obs,pred){
  f = (!is.na(obs) & !is.na(pred))
  res = sum((obs-pred)^2)
  Var = sum((obs-mean(obs))^2)
  1-res/Var
}

fRMSE = function(obs, pred){
  f = (!is.na(obs) & !is.na(pred))
  N = length(obs)
  RMSE = sqrt(sum((obs-pred)^2)/N)
  return(RMSE)
}


fPreprocessX <- function(p,x, selX=NULL){
  
  p[3:4] <- floor(p[3:4])
  param <- list(bl.lambda=p[1], bl.p=p[2],
               sg.m=p[3], sg.p=sum(p[3:4]), sg.w=9)
  
  xabs <- max(x)-x
  xb   <- baseline.als(xabs, lambda=param$bl.lambda, p=param$bl.p)
  xbsg <- savitzkyGolay(xb$corrected, m=param$sg.m, p=param$sg.p, w=param$sg.w)
  
  xsel <- NULL
  if(!is.null(selX)){
    selXcol <- unique(floor(p[-(1:4)]))
    selXcol <- selXcol[selXcol<ncol(xbsg)]
    xsel    <- xbsg[,selXcol]
    param   <- c(param, selXcol=NULL)
    param$selXcol <- selXcol
  }
  
  result <- list(xb = xb, xbsg = xbsg, xsel = xsel, param = param)
  return(result)
}

fGAPLS <- function(p){
  #function optimesed by GA
  #p:Hyperparameter(vector)
  #p[1]~p[2]:baseline.als()
  #p[3]~p[4]:savitzky-Golay
  #p[5]~:select explanatory variable

  # preprocess
  xpre <- fPreprocessX(x=Xtrain, p=p, selX=T)
  xpls <- xpre$xsel

  # model
  pls <- plsr(Ytrain ~ xpls, scale=FALSE, validation="LOO")
  
  # val R2
  R2.comp <- apply(pls$fitted.values[,1,],2,fR2,obs=Ytrain)
  ncomp   <- which.min(R2.comp)
  
  returnGA <- R2.comp[ncomp]
  return(returnGA)
  
}

fPLS <- function(p){
  #function optimesed by GA
  #p:Hyperparameter(vector)
  #p[1]~p[2]:baseline.als()
  #p[3]~p[5]:savitzky-Golay
  xpre <- fPreprocessX(x=Xtrain, p=p, selX=T)
  xpls <- xpre$xbsg
  
  # model
  pls <- plsr(Ytrain ~ xpls, scale=FALSE, validation="LOO")
  
  # val R2
  R2.comp <- apply(pls$fitted.values[,1,],2,fR2,obs=Ytrain)
  ncomp <- which.min(R2.comp)
  
  returnGA <- R2.comp[ncomp]
  return(returnGA)
}



# GAPLS ----

# GA input parameter
bl.lambda.range <- c(1,10)
bl.p.range      <- c(0.001,0.1)
sg.m.range      <- c(0,2.99)
sg.p.range      <- c(1,7)

nExVal          <- 4
ExVal.range     <- c(1, ncol(Xtrain)-1)

GAfitParRange <- cbind(bl.lambda.range, bl.p.range, sg.m.range, sg.p.range,
                       matrix(rep(ExVal.range,times=nExVal),2,nExVal))

fit.Param <- ga(type="real-valued", fGAPLS, maxiter=100, seed=1,
                lower=GAfitParRange[1,], upper=GAfitParRange[2,])

plot(fit.Param)
fit.Param@solution



#> result ----
Xfit <- fPreprocessX(p=fit.Param@solution[1,], X0, selX=T)
matplot(t(Xfit$xbsg),col=ycol,type="l",main="ベースライン + ノイズ除去")
abline(v=Xfit$param$selXcol, lty=2, col="grey70")

fitpls  <- plsr(Ytrain ~ Xfit$xsel[trainSel,], scale=F)
plot(fitpls, "validation")
ncomp <- 1

Yfit  <- fitpls$fitted.values[,1,ncomp]
Ypred <- predict(fitpls, newdata=Xfit$xsel[!trainSel,])[,1,ncomp]

# 実測予測プロット
{
  ranges <- range(Ytrain,Yfit,Ytest,Ypred)
  plot(Ytrain,Yfit,pch=5,
       xlab="obs", ylab="pred", xlim=ranges, ylim=ranges)
  points(Ytest,Ypred,pch=4,col=2)
  abline(a=0, b=1, lty=2)
  
  R2.train   <- round(fR2(obs=Ytrain, pred=Yfit),3)
  RMSE.train <- round(fRMSE(obs=Ytrain, pred=Yfit),3)
  RMSE.test  <- round(fRMSE(obs=Ytest,  pred=Ypred),3)
  legends    <- c(paste0("R2.train = ", R2.train),
                  paste0("RMSE.train = ",RMSE.train),
                  paste0("RMSE.test = ",RMSE.test))
  legend("topleft", legend=legends, box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
  legend("bottomright", legend=c("train","test"),col=1:2,pch=c(5,4),cex=1.5,
         box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
}

おまけ

GAPLSが有効だったかどうかを調べるためには全変数ver.での精度も確かめておくとより説得性も増すでしょう。ということで前処理のみGA最適化するコードも置いておきます。

Code All(前処理ハイパラのみのGA)

preprocessGA.R


fPLS <- function(p){
  #function optimesed by GA
  #p:Hyperparameter(vector)
  #p[1]~p[2]:baseline.als()
  #p[3]~p[4]:savitzky-Golay
  
  xpre <- fPreprocessX(x=Xtrain, p=p, selX=F)
  pls <- plsr(Ytrain ~ xpre$xbsg, scale=FALSE, validation="LOO")
  
  # val R2
  R2.comp <- apply(pls$fitted.values[,1,],2,fR2,obs=Ytrain)
  ncomp <- which.min(R2.comp)
  
  returnGA <- R2.comp[ncomp]
  return(returnGA)
}

# GA input parameter
bl.lambda.range <- c(1,10)
bl.p.range      <- c(0.001,0.1)
sg.m.range      <- c(0,2.99)
sg.p.range      <- c(1,7)

GAfitParRange <- cbind(bl.lambda.range, bl.p.range, sg.m.range, sg.p.range)

fit.Param <- ga(type="real-valued", fPLS, maxiter=10, seed=1,
                lower=GAfitParRange[1,], upper=GAfitParRange[2,])
plot(fit.Param)
fit.Param@solution

#> result ----
Xfit <- fPreprocessX(p=fit.Param@solution[1,], X0, selX=F)
matplot(t(Xfit$xbsg),col=ycol,type="l",main="ベースライン + ノイズ除去")
abline(v=Xfit$param$selXcol, lty=2, col="grey70")

fitpls  <- plsr(Ytrain ~ Xfit$xsel[trainSel,], scale=F)
plot(fitpls, "validation")
ncomp <- 5 

Yfit  <- fitpls$fitted.values[,1,ncomp]
Ypred <- predict(fitpls, newdata=Xfit$xsel[!trainSel,])[,1,ncomp]

# 実測予測プロット
{
  ranges <- range(Ytrain,Yfit,Ytest,Ypred)
  plot(Ytrain,Yfit,pch=5,
       xlab="obs", ylab="pred", xlim=ranges, ylim=ranges)
  points(Ytest,Ypred,pch=4,col=2)
  abline(a=0, b=1, lty=2)
  
  R2.train   <- round(fR2(obs=Ytrain, pred=Yfit),3)
  RMSE.train <- round(fRMSE(obs=Ytrain, pred=Yfit),3)
  RMSE.test  <- round(fRMSE(obs=Ytest,  pred=Ypred),3)
  legends    <- c(paste0("R2.train = ", R2.train),
                  paste0("RMSE.train = ",RMSE.train),
                  paste0("RMSE.test = ",RMSE.test))
  legend("topleft", legend=legends, box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
  legend("bottomright", legend=c("train","test"),col=1:2,pch=c(5,4),cex=1.5,
         box.col=rgb(0,0,0,0), bg=rgb(0,0,0,0))
}

参考

金子研究室にGAPLSについての丁寧な説明があります⇒ https://datachemeng.com/gaplsgasvr/ ↩
Swierenga H., de Weijer A. P., van Wijk R. J., Buydens L. M. C. (1999) Strategy for constructing robust multivariate calibration models Chemometrics and Intelligent Laboratoryy Systems, 49(1), 1–17. ↩
GA側で制御できるのかもしれませんが、記事作成時点では調べても方法がわからなかったのでこのような処理にしています。 ↩
正直なところ多項式補間の次数と窓幅はそこまで予測精度に大きな影響を与えないと感じているのでやや雑な処理になっています。もっと正確で良い方法はあると思うので知っている方は教えてもらえると嬉しいです。 ↩
容易に過学習すると思いますので出力を何にするかはよく検討すべき項目かと思います。参考文献ではQ2値等が使用されています。 ↩
化学の世界では殆どのケースでランダム性が担保されない・点数が少ない・交絡あり、というデータセットかと思いますので、考えなしのランダムCVやLOOではまずい処理になる可能性が高いのではと感じています。本記事を書くモチベーションの9割くらいはこの問題の解消にあります。 ↩