Stackingに関して（２）
前回の続き。
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")
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
> 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))
では、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を行わなかった場合のモデルとの比較等に関しては省略。
以上。
