読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

Stackingに関して(2)

前回の続き。
highschoolstudent.hatenablog.com


サンプルデータの作成。あまり細かいことは考えず適当なパラメータ設定で適当に作成。

library(MASS)		# 乱数の作成(多変量正規分布)
library(glmnet)		# Ridge回帰
library(lme4)		# 混合モデル
library(lattice)	# lme4による分析結果の可視化
n_obs <- 2000	# サンプル数
n_r <- 20		# rの水準数
n_x <- 30		# Xの列数

n_r_each <- n_obs/n_r	# 各rの水準のサンプル数


###  r(random effect)の作成  ###

s <- 5	# random effectの標準偏差

# rのmeanの生成(20水準)
set.seed(123)
mean_r <- floor(rnorm(n=n_r, mean=0, sd=15))
mean_r

# mean_rのそれぞれの要素を20回繰り返す
mean_r_vec <- rep(mean_r, times=1, each=n_r_each)	# n_obs行のベクトル

# 乱数(sd=sの正規分布)を加える
set.seed(123)
r_vec <- mean_r_vec + rnorm(n=n_obs, mean=0, sd=s)

# 各rのID列の作成
r_id=rep(1:n_r, times=1, each=n_r_each)


###  X(2000×30)の作成  ###

mu <- rep(0, n_x)		# 各変数の平均は0
sig <- diag(n_x)		# 相関係数行列(30×30)
sig[sig==0] <- 0.1		# 30の変数間に0.1の相関を与える
sig_var <- sig*4		# 分散共分散行列。分散を4とする
set.seed(123)
x_mtrx <- mvrnorm(n_obs, mu=mu, Sigma=sig_var)	# 1000*30の行列

X <- data.frame(x_mtrx)


###  目的変数yの作成  ###

# すべての係数(β)は1とする。重み付け無しで30列を足し合わせる
logit_p <- apply(X, MARGIN=1, sum)

# randon effectを加える
logit_r_p <- logit_p + r_vec

# scaleしておく
p <- 1 / (1+exp(-1*scale(logit_r_p)))		

# yの生成
set.seed(123)
y <- rbinom(n=n_obs, size=1, prob=p)


### data.frameにまとめる ###

data1 <- data.frame(r_id, mean_r_vec, X, p, y)
# 列IDを付けておく
data1 <- data.frame(OBS_ID=1:nrow(data1), data1)
# データ型の変更
data1[, "r_id"] <- as.factor(data1[, "r_id"])
data1[, "y"] <- as.factor(data1[, "y"])
data1[, "OBS_ID"] <- as.character(data1[, "OBS_ID"])
> str(data1)
'data.frame':   2000 obs. of  35 variables:
 $ OBS_ID    : chr  "1" "2" "3" "4" ...
 $ r_id      : Factor w/ 20 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ mean_r_vec: num  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
 $ X1        : num  -0.267 2.431 -0.354 -1.632 -2.469 ...
 $ X2        : num  0.0797 -0.6252 -3.987 -0.219 0.9425 ...
 $ X3        : num  -2.678 1.573 0.518 -3.111 0.12 ...
 $ X4        : num  3.29 -0.87 -7.05 1.1 -1.74 ...
 $ X5        : num  3.235 -0.832 2.129 -1.876 2.369 ...
 $ X6        : num  -1.846 -0.7481 0.0862 0.16 -2.193 ...
 $ X7        : num  2.383 -2.76 2.649 0.34 -0.901 ...
 $ X8        : num  -1.725 -1.1587 -1.2911 0.325 0.0011 ...
 $ X9        : num  0.0625 1.6856 0.11 -0.6957 3.046 ...
 $ X10       : num  -0.0845 2.1115 -2.3013 3.282 1.1797 ...
 $ X11       : num  0.6461 1.1761 -3.7507 1.4679 0.0739 ...
 $ X12       : num  1.024 2.317 -1.685 -0.367 -2.538 ...
 $ X13       : num  -1.27 -3.483 -0.253 0.399 2.68 ...
 $ X14       : num  -2.569 -0.559 -3.942 -3.517 -0.387 ...
 $ X15       : num  2.5588 -2.0188 0.0141 0.8221 0.2959 ...
 $ X16       : num  2.227 2.08 -0.565 0.555 0.261 ...
 $ X17       : num  -3.413 1.206 1.578 0.226 -1.649 ...
 $ X18       : num  3.189 -0.286 -1.177 1.409 -1.41 ...
 $ X19       : num  4.053 0.355 0.173 -2.395 0.607 ...
 $ X20       : num  -0.356 -0.632 -1.924 0.464 1.393 ...
 $ X21       : num  2.0457 -0.4395 0.7572 -2.1578 -0.0979 ...
 $ X22       : num  2.877 1.097 -1.977 1.85 0.151 ...
 $ X23       : num  -0.447 -0.177 -1.276 -1.001 -0.904 ...
 $ X24       : num  -0.506 0.924 -0.718 0.292 1.549 ...
 $ X25       : num  2.343 0.305 -3.97 -3.443 -0.222 ...
 $ X26       : num  -1.731 -0.769 2.971 2.385 -0.441 ...
 $ X27       : num  1.6675 -0.0147 -2.1689 1.4975 -2.7424 ...
 $ X28       : num  0.0978 3.6201 -2.9536 1.036 1.4886 ...
 $ X29       : num  -2.901 0.138 -3.443 -1.862 -1.06 ...
 $ X30       : num  0.1366 -0.6662 0.0819 3.1395 -0.1976 ...
 $ p         : num  0.491 0.429 0.162 0.374 0.364 ...
 $ y         : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 1 1 1 ...
plot(data1[,"r_id"],data1[,"p"], xlab="r_id", ylab="p")

f:id:High_School_Student:20151230135232j:plain
rによって、適当に散らばったデータ。


Train(Stacking用に分割), Testセットの作成。

Train:Test=1600:400。Trainは8分割し、各Stacking用データは200とする。

ID列(dataSet_ID)を作成し、管理。
Stacking用に8分割したTrainセットは1,2,3,4,5,6,7,8とIDを付け、Testセットを99とする。

# Train(Stacking), Test IDの作成
# Train -> 1,2,3,4,5,6,7,8
# Test -> 99
each_ID <- c(rep(1,10),rep(2,10),rep(3,10),rep(4,10),rep(5,10),rep(6,10),rep(7,10),rep(8,10),rep(99,20))
dataSet_ID <- rep(each_ID, n_r)

data2 <- data.frame(data1, dataSet_ID)

# Train, Test dataの分割
trainData <- data2[data2[,"dataSet_ID"]!=99, ]
testData <- data2[data2[,"dataSet_ID"]==99, ]


では、Ridge回帰で、1stステージモデルを作成する。
glmnetパッケージを使用。
highschoolstudent.hatenablog.com

# Trainデータの説明変数
exp_vars <- as.matrix(trainData[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
					"X11","X12","X13","X14","X15","X16","X17","X18","X19","X20",
					"X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")]
)

# Trainデータの目的変数
#target_var <- trainData[,"y"]
target_var <- as.vector(trainData[,"y"])

# Testデータの説明変数
exp_vars_test <- as.matrix(testData[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
					"X11","X12","X13","X14","X15","X16","X17","X18","X19","X20",
					"X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")]
)

まず、Cross Validationで最適なλを見つける。

cv.glmnetは、k+1回数glmnetを走らせる。
k回でλを決定し、残り1回でパラメータを推定する。

前回述べた通り、この分割はStackingの分割と同とする必要はない。
しかし、foldid引数にIDが指定できる、8分割で十分だろうという事で、Cross ValidatioとStackingの分割は同じとして処理する。

set.seed(123)
fit1st_ridgeCV <- cv.glmnet(
	x=exp_vars,
	y=target_var,
	family="binomial",
	#nfold=8,
	foldid=trainData[,"dataSet_ID"],	# CVのIDを指定
	keep=TRUE,	# CVのIDを保存
	alpha=0
)

plot(fit1st_ridgeCV)

# CV(k=8)による、最適なλ
lambdaMin <- fit1st_ridgeCV$lambda.min

f:id:High_School_Student:20151230140019j:plain

> lambdaMin
[1] 0.1540261

見つけたλ(lambdaMin)によって、2ndモデル用の変数(Ridge回帰による予測値)を作成。

Trainデータの処理。
8個に分けた各Stackingセットを一つ抜いて、残りのデータで固定したλ(lambdaMin)でglmnetを実行し、パラメータを推定、抜いておいたセットにこの推定したモデルをあてはめて、予測値を計算する。
これをStackingに分けた数(8回)実行する。

# Stackingセット(= CVセット)のID
cv_id <- as.integer(levels(as.factor(trainData[,"dataSet_ID"])))
# 1 2 3 4 5 6 7 8

# 8回のあてはめ結果オブジェクトの格納用リスト
list_res <- list()

# 予測値格納用のdf(col1:OBS_ID、col2:予測obsのID、col3:予測値)
pred_y <- NULL

for(id in cv_id){
	#print(paste("cv_id:",id))

	# パラメータ推定用データ
	train <- trainData[trainData[,"dataSet_ID"]!=id, ]
	#print("train - cv_id")
	#print(table(train[,"dataSet_ID"]))

	x <- as.matrix(train[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
				"X11","X12","X13","X14","X15","X16","X17","X18","X19","X20",
				"X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")])
	y <- as.vector(train[,"y"])
	
	# あてはめ用データ
	target <- trainData[trainData[,"dataSet_ID"]==id, ]
	#print("target - cv_id")
	#print(table(target[,"dataSet_ID"]))
	
	x_target <- as.matrix(target[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
					"X11","X12","X13","X14","X15","X16","X17","X18","X19","X20",
					"X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")])

	fitRidge <- glmnet(	# 指定したλ(lambdaMin)でリッジ回帰
		x=x, 
		y=y, 
		family="binomial", 
		lambda=lambdaMin, 
		alpha=0 )
		
	list_res[[id]] <- fitRidge	# あてはめ結果オブジェクトの保存
	
	# あてはめように抜いてあったデータに、パラメータ推定用データで推定したモデルをあてはめて、予測値を推定。あてはめデータも保存。
	pred_y <- rbind(pred_y, data.frame(target[,"OBS_ID"], mean(target[,"dataSet_ID"]), predict(fitRidge, newx=x_target)))
	
	#print("---------------------")
}
names(pred_y) <- c("OBS_ID", "target_id", "pred_input")


# リッジ回帰による予測で作成した列を含めたdata.frame(Trainデータ)の作成
trainData2 <- merge(x=trainData, y=pred_y, by.x="OBS_ID", by.y="OBS_ID", sort=FALSE, all=TRUE)

# 列の絞り込みと順序の入れ替え
trainData2 <- trainData2[,
c("OBS_ID","r_id","mean_r_vec",
"X1","X2","X3","X4","X5","X6","X7","X8","X9","X10",
"X11","X12","X13","X14","X15","X16","X17","X18","X19","X20",
"X21","X22","X23","X24","X25","X26","X27","X28","X29","X30",
"p","y","dataSet_ID","pred_input")]

Trainデータの処理。
Trainデータ全体で、固定したλ(lambdaMin)でglmnetを実行し、パラメータを推定、それにTestデータをあてはめ予測値を作成。

fitRidge_test <- glmnet(	# 指定したλ(lambdaMin)でリッジ回帰
	x=exp_vars,
	y=target_var,
	family="binomial", 
	lambda=lambdaMin, 
	alpha=0
)

# Testデータでの予測値
ypred <- predict(fitRidge_test, newx=exp_vars_test)

# リッジ回帰による予測で作成した列を含めたdata.frame(Testデータ)の作成
testData2 <- data.frame( testData, pred_input=as.vector(ypred) ) 

TrainとTestデータを一つにまとめる。

data3 <- rbind(trainData2, testData2)

# pred_input(線形結合)をpに直す
ypredp <- 1 / (1 + exp(-1*data3[,"pred_input"]))

data3 <- data.frame(data3, pred_input_p=ypredp)

以上、1stステージモデル。

では、2ndステージモデルに入る。

# Train, Test dataの分割
Sec_trainData <- data3[data3[,"dataSet_ID"]!=99, ]
Sec_testData <- data3[data3[,"dataSet_ID"]==99, ]

以前、混合モデルを使用した際は、glmmMRを使ったが、実際はlme4パッケージの方がメジャーなようなので、こちらを使う。
highschoolstudent.hatenablog.com

predict関数もサポートしており、変量効果の各水準に対する予測値(BLUP(Best Linior Unbiased Predictor))を考慮した予測が可能。

詳細に関しては、lme4: Mixed-effects modeling with R 参照。

ランダム切片モデルとする。

# Intercept model
fit_mixed <- glmer(
	y ~
	pred_input +
	(1|r_id),	# 1は切片の意味
	
	data=Sec_trainData,
	family=binomial, 
	verbose=1
)
> summary(fit_mixed)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: y ~ pred_input + (1 | r_id)
   Data: Sec_trainData

     AIC      BIC   logLik deviance df.resid 
  1978.0   1994.2   -986.0   1972.0     1597 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8686 -0.8164 -0.3042  0.8195  3.2091 

Random effects:
 Groups Name        Variance Std.Dev.
 r_id   (Intercept) 0.3299   0.5743  
Number of obs: 1600, groups:  r_id, 20

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.003755   0.139643  -0.027    0.979    
pred_input   1.168448   0.094078  12.420   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr)
pred_input 0.017 


ranef関数に、結果オブジェクトを引数で渡すと、変量効果の各水準に対する予測値(BLUP(Best Linior Unbiased Predictor))を出力できる。

condVar=TRUEとすると、各変量効果の水準の予測値の95%予測区間も出力する。

dotplot(ranef(fit_mixed, condVar=TRUE))

f:id:High_School_Student:20151230140647j:plain

では、2ndステージモデルの予測値を出力する。

# Trainデータの予測値
predTrain_fit_mixed <- 1 / ( 1 + exp(-1*predict(fit_mixed, newdata=Sec_trainData)) )

# Testデータの予測値
predTest_fit_mixed <- 1 / ( 1 + exp(-1*predict(fit_mixed, newdata=Sec_testData)) )

# data.frameにまとめる
Sec_trainData2 <- data.frame(Sec_trainData, pred_y=predTrain_fit_mixed )
Sec_testData2 <- data.frame(Sec_testData, pred_y=predTest_fit_mixed)

data4 <- rbind(Sec_trainData2, Sec_testData2)


上で使用したprediction関数であるが、引数re.form=NULL(デフォルト)としておくと、変量効果の予測値(BLUP)を含めた予測となる。

allow.new.levelsという引数があり、allow.new.levels=TRUE(デフォルトはFALSE)とすると、新しいデータに対する予測の際、モデル推定のデータになかった変量効果の水準が含まれた場合、変量効果の予測を考慮しない値が返される。
デフォルト(FALSE)だと、モデル推定のデータになかった変量効果の水準が含まれた場合、Errorが返るとのこと。

> head(data4)
  OBS_ID r_id mean_r_vec         X1          X2         X3         X4         X5
1      1    1         -9 -0.2670289  0.07965359 -2.6781329  3.2920218  3.2347046
2      2    1         -9  2.4309568 -0.62516429  1.5726213 -0.8702877 -0.8316664
3      3    1         -9 -0.3538252 -3.98699188  0.5181232 -7.0503276  2.1287953
4      4    1         -9 -1.6324691 -0.21901088 -3.1107837  1.1008922 -1.8763152
5      5    1         -9 -2.4687447  0.94246812  0.1196235 -1.7423328  2.3693425
6      6    1         -9  1.2053775 -2.38738177  0.7406409 -2.9720772 -1.2437253

### 省略 ###

         X27         X28        X29         X30         p y dataSet_ID  pred_input
1  1.6675496  0.09780506 -2.9014770  0.13657163 0.4909926 0          1  0.49382763
2 -0.0146797  3.62010783  0.1383445 -0.66618073 0.4288412 1          1  0.05697817
3 -2.1689243 -2.95355951 -3.4428897  0.08190659 0.1619845 0          1 -0.89058592
4  1.4974864  1.03595553 -1.8618701  3.13950267 0.3741005 1          1  0.02697507
5 -2.7424436  1.48861771 -1.0603866 -0.19760381 0.3637205 1          1 -0.26863659
6  0.1702958  1.57345489 -2.8409584 -1.65512597 0.1465225 0          1 -1.05824906
  pred_input_p    pred_y
1    0.6210077 0.5308346
2    0.5142407 0.4044539
3    0.2909889 0.1833040
4    0.5067434 0.3960387
5    0.4332418 0.3170407
6    0.2576442 0.1557724

pred_yが最終予測値である。


なお、パフォーマンスに関してはAUCで確認している。
結果だけだが、Trainデータで0.7404149、Testデータで0.7289979となる。

Stackingを行わなかった場合のモデルとの比較等に関しては省略。


以上。