stackingを試してみた

つい先日，stackingについての以下の記事が話題になっていました．

このあたり，私自身は試したことがなかったので，実際に試してみましたよというお話．

コード

Rでちゃちゃっと書きました．データをk分割して，分割サンプルごとに訓練データでロジスティック回帰予測モデル構築→予測結果を訓練データに追加→RandomForestで予測モデルを構築，までが訓練フェーズ．テストフェーズでは構築したロジスティック回帰とRandomForestを使ってテストデータの分類を行いました．

使ったのは手元にあった2000サンプル，分類クラス数は2クラスで，正例負例がそれぞれ1000サンプルでした．素性ベクトルはすべてバイナリの15個の変数になります．

library(randomForest)

# load data
data = read.delim("data/sample.tsv", sep="\t")

# stacking sample
k = 10 # cross validation split number
result = c()
for (i in 1:k) {
  print(i)
  data.splitted = cv(data, k)
  # construct predict model
  data.train = cv.training(data.splitted, 1)
  model.glm = glm(y~., data=data.train, family=binomial)
  data.train = stacking(data.train, model.glm) # commentout this line when not using stacking
  model.rf = randomForest(y~., data=data.train)
  # predict with test data
  data.test = cv.test(data.splitted, 1)
  data.test = stacking(data.test, model.glm) # commentout this line when not using stacking
  model.rf.predict = predict(model.rf, newdata=data.test, type="class")
  result = rbind(result, score(model.rf.predict, data.test$y))
}

上で説明した処理はすべて関数化して，以下のようにまとめました．また判別制度の評価をちゃんとするには交差検定も必要なので，こちらも関数にしてまとめておきました*1．

# create data for k-fold cross validation
cv = function(d, k) {
  n = sample(nrow(d), nrow(d))
  d.randomized = data[n,] # randomize data
  n.residual = k-nrow(d)%%k
  d.dummy = as.data.frame(matrix(NA, nrow=n.residual, ncol=ncol(d)))
  names(d.dummy) = names(d)
  d.randomized = rbind(d.randomized, d.dummy) # append dummy for residuals
  d.splitted = split(d.randomized, 1:k)
  for (i in 1:k) {
    d.splitted[[i]] = na.omit(d.splitted[[i]])
  }
  d.splitted
}

# training data
cv.training = function(d, k) {
  d.train = as.data.frame(matrix(0, nrow=0, ncol=ncol(d[[1]])))
  names(d.train) = names(d[[1]])
  for (i in 1:length(d)) {
    if (i != k) {
      d.train = rbind(d.train, d[[i]])
    }
  }
  d.train
}

# test data
cv.test = function(d, k) {
  d[[k]]
}

# stacking with glm
stacking = function(d, m) {
  d = cbind(d, predict(m, newdata=d, type="response"))
  names(d)[length(d)] = "stacking"
  d
}

# check
score = function(p, r) {
  s = c(0, 0, 0, 0)
  for (i in 1:length(p)) {
    pi = 2-as.integer(p[[i]])
    ri = 2-as.integer(r[i])
    s[pi*2+ri+1] = s[pi*2+ri+1]+1
  }
  s
}

結果

ということで，10分割の交差検定をした結果を見ましょう．以下のように予測値と実際の結果のマトリックスで表示してみます．

# show results
m = matrix(apply(result, 2, sum), 2, 2)
dimnames(m) = list(c("res$n", "res$p"), c("pred$n", "pred$p"))
print(m)
print(m/nrow(data))

stackingなし

stacking処理を無効にするには，単にstacking()関数を読んでいる部分2箇所をコメントアウトしてあげるだけです．結果は．Precisionは850/(850+149)=83.2%で，Recallは850/(850+134)=86.3%でした．

> print(m)
       res$p res$n
pred$p   850   149
pred$n   134   867
> print(m/nrow(data))
       res$p  res$n
pred$p 0.425 0.0745
pred$n 0.067 0.4335

stackingあり

それに対してstackingはどうでしょう．Precisionは827/(827+167)=83.2%で，Recallは827/(827+152)=84.5%でした．むしろ悪くなってる...orz

> print(m)
       res$p res$n
pred$p   827   167
pred$n   152   854
> print(m/nrow(data))
        res$p  res$n
pred$p 0.4135 0.0835
pred$n 0.0760 0.4270

というわけで，とりあえず試してみましたよというお話でした．コード全体はgistにあげました．よろしければどうぞ．

*1:このあたりをまとめるにあたっては，shakezoさんの RでデータフレームをK分割する - shakezoの日記 を参考にさせていただきました．