GLMMを勉強してみる。 - R glmmML
GLMM（リンク関数＝ロジスティック関数、分布＝二項分布）を勉強してみる。RのglmmML使用。
久保(2012) の7章を参考。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (21件) を見る
ロジスティック回帰モデルの復習。
これにランダム効果を加えると。
文献で取り上げている種子の例でなく、自分の事例に合ったサンプルデータを作成してみる。
「インサイドセールス20人が顧客からインバウンドコールを受ける。セールスを行い、契約につながったかどうかを記録する。取得できる顧客情報は性別の年齢区分。顧客情報と契約率との関係を分析する。インサイドセールスによって、パフォーマンスは異なる。（インサイドセールスによる影響をランダム効果とする。）」
1200行のデータ。1200件のインバウンドコールのデータとなる。
n_row <- 1200
行番号。
Row_ID <- 1:n_row
インサイドセールスの作成。
n1 <- n_row * 0.7
n2 <- n_row * 0.3
salesID <- sort( c( rpois(n=n1, lambda=7)+1, sample(1:20, n2, replace=TRUE)) )
salesID_count <- table( salesID )
salesID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 14 28 26 86 88 133 144 149 116 111 63 48 39 36 18 18 20 24 14 25
各20人のインサイドセールスが受けたインバウンドコールの数。
インサイドセールスの違いより発生するランダム効果。
s <- 2 ## rand effect真のパラメータ
randEff <- rnorm(n=20, mean=0, sd=s)
mean(randEff)
[1] -0.007712367
sd(randEff) ## rand effect実際のパラメータ
[1] 2.163217
各インサイドセールスによるランダム効果の列の作成。
randEff_col <- c(NULL)
for(i in 1:length(randEff)){
randEff_col <- append( randEff_col, rep(randEff[i], salesID_count[i]) )
}
顧客の性別の作成。男女比、M:F = 3:7、となるように作成。
rand_num_sex <- runif(n=n_row, min = 0, max = 1)
fn_sex <- function(i){
if(i<0.3){return("M")}
else{return("F")}
}
sex_col <- sapply(rand_num_sex, fn_sex)
顧客の年代の作成。年代比、40以下:60以下:60以上 = 15:30:55で作成。
rand_num_age <- runif(n=n_row, min = 0, max = 1)
fn_age <- function(i){
if(i<0.15){return("-40")}
else if(i<0.45){return("-60")}
else{return("60-")}
}
age_col <- sapply(rand_num_age, fn_age)
二項分布のパラメータ（q）の作成。
真のパラメータ定義を行う。
切片
beta0 <- -4 # intercept
性別Mに対するFの効果
beta1 <- 2 # sex="F"
年代40以下に対する、60以下の効果
beta2_1 <- 0.5 # age="-60"
年代40以下に対する、60以上の効果
beta2_2 <- 1.5 # age="60-"
性別、年代、インサイドセールスの効果を入れて、ロジスティック関数によりqを出力する。
fn_q <- function(id){
if(sex_col[id]=="F"){effect_sex <- beta1}
else{effect_sex <- 0}
if(age_col[id]=="-60"){effect_age <- beta2_1}
else if(age_col[id]=="60-"){effect_age <- beta2_2}
else{effect_age <- 0}
lin <- beta0 + effect_sex + effect_age + randEff_col[id]
return( 1 / (1+exp(-lin)) )
}
q_col <- sapply(Row_ID, fn_q)
このqを二項分布のパラメータとし、確率変数yを出力。
fn_y <- function(q){
return( rbinom(n=1, size=1, prob=q) )
}
y_col <- sapply(q_col, fn_y)
Train, Testデータを1000:200で分ける。
td1 <- n_row * 1000/1200
td2 <- n_row * 200/1200
test_col0 <- c( rep("Train", td1), rep("Test", td2) )
# ランダム化
test_col <- sample(test_col0, n_row, replace=FALSE)
データフレームにまとめる。
data0 <- data.frame(Row_ID, Sales_ID=salesID, SalesRandEffect=randEff_col, Sex=sex_col, Age=age_col, q=q_col, y=y_col, Test=test_col)
# モデルをあてはめた際の係数の解釈上、カテゴリ変数Sexの順序を入れ替えておく
data0$Sex <- factor(data0$Sex, levels=c("M","F"))
data0
Row_ID Sales_ID SalesRandEffect Sex Age q y Test 1 1 1 -2.3580976 M 60- 0.0077054077 0 Train 2 2 1 -2.3580976 F -60 0.0206717743 0 Train 3 3 1 -2.3580976 M -40 0.0017296628 0 Train 4 4 1 -2.3580976 F 60- 0.0542642460 0 Train 5 5 1 -2.3580976 F 60- 0.0542642460 0 Train 6 6 1 -2.3580976 F 60- 0.0542642460 0 Train 7 7 1 -2.3580976 M 60- 0.0077054077 0 Train 8 8 1 -2.3580976 F 60- 0.0542642460 0 Train 9 9 1 -2.3580976 M 60- 0.0077054077 0 Train 10 10 1 -2.3580976 F 60- 0.0542642460 0 Train 11 11 1 -2.3580976 M 60- 0.0077054077 0 Train 12 12 1 -2.3580976 F 60- 0.0542642460 1 Train .............................. skip 1191 1191 20 1.2576302 F 60- 0.6808390084 0 Test 1192 1192 20 1.2576302 F 60- 0.6808390084 1 Train 1193 1193 20 1.2576302 F 60- 0.6808390084 1 Train 1194 1194 20 1.2576302 F 60- 0.6808390084 0 Test 1195 1195 20 1.2576302 F 60- 0.6808390084 1 Train 1196 1196 20 1.2576302 F 60- 0.6808390084 0 Train 1197 1197 20 1.2576302 M -40 0.0605190257 0 Train 1198 1198 20 1.2576302 M 60- 0.2240237610 0 Train 1199 1199 20 1.2576302 F -60 0.4397024392 1 Train 1200 1200 20 1.2576302 F -60 0.4397024392 1 Test
SalesRandEffectはランダム効果なので、もちろん観測されない。
qは真の二項分布のパラメータ（Sex, Age, SalesRandEffectに依存）となるが、これももちろん観測はされない。
TrainデータとTestデータを分ける。
data1 <- data0[data0$Test=="Train",]
data1_test <- data0[data0$Test=="Test",]
分布の確認。
summary(data1)
Row_ID Sales_ID SalesRandEffect Sex Age q
Min. : 1.0 Min. : 1.00 Min. :-3.1455 M:301 -40:164 Min. :0.0007878
1st Qu.: 299.8 1st Qu.: 6.00 1st Qu.:-1.6508 F:699 -60:306 1st Qu.:0.0380138
Median : 604.5 Median : 8.00 Median :-0.6385 60-:530 Median :0.1209089
Mean : 603.0 Mean : 8.75 Mean : 0.5669 Mean :0.3367441
3rd Qu.: 904.2 3rd Qu.:11.00 3rd Qu.: 1.6061 3rd Qu.:0.6569576
Max. :1199.0 Max. :20.00 Max. : 5.7100 Max. :0.9945678
y Test
Min. :0.000 Test : 0
1st Qu.:0.000 Train:1000
Median :0.000
Mean :0.341
3rd Qu.:1.000
Max. :1.000
summary(data1_test)
Row_ID Sales_ID SalesRandEffect Sex Age q
Min. : 14.0 Min. : 1.00 Min. :-3.1455 M: 61 -40: 25 Min. :0.003233
1st Qu.: 316.8 1st Qu.: 6.00 1st Qu.:-1.6508 F:139 -60: 58 1st Qu.:0.041059
Median : 580.0 Median : 8.00 Median :-0.6385 60-:117 Median :0.177129
Mean : 588.2 Mean : 8.55 Mean : 0.7733 Mean :0.365478
3rd Qu.: 889.5 3rd Qu.:10.00 3rd Qu.: 4.5274 3rd Qu.:0.740155
Max. :1200.0 Max. :20.00 Max. : 5.7100 Max. :0.994568
y Test
Min. :0.000 Test :200
1st Qu.:0.000 Train: 0
Median :0.000
Mean :0.355
3rd Qu.:1.000
Max. :1.000
JMPでグラフを描く。
上段がTestデータで、下段がTrainデータ。
y=1、つまり契約、を強調表示してある。
table(data1$y)
0 1 659 341
table(data1_test$y)
0 1 129 71
各データセットで、1（契約）の割合の確認。
# % of 1 (data1)
table(data1$y)[2]/sum(table(data1$y))
0.341
# % of 1 (data1_test)
table(data1_test$y)[2]/sum(table(data1_test$y))
0.355
まずはランダム効果を考慮せず、通常のロジスティック回帰（GLM）を試してみる。
説明変数はSexとAge。
res_glm1 <- glm(y~Sex+Age, data=data1, family=binomial)
Call: glm(formula = y ~ Sex + Age, family = binomial, data = data1) Coefficients: (Intercept) SexF Age-60 Age60- -1.31114 0.60076 0.06984 0.36629 Degrees of Freedom: 999 Total (i.e. Null); 996 Residual Null Deviance: 1283 Residual Deviance: 1262 AIC: 1270
真のパラメータはInterceptが-4、Sex Fが2、Age -60が0.5、Age 60-が1.5と設定してあるので、+/-の方角は合っているが、値はずれている感じがする。
Trainデータへ作成したモデルをあてはめ、qを出力。
pred_glm0 <- predict(res_glm1)
pred_glm <- 1 / (1 + exp(-1*pred_glm0))
Testデータへ作成したモデルをあてはめ、qを出力。
pred_glm_test0 <- predict(res_glm1, newdata=data1_test)
pred_glm_test <- 1 / (1 + exp(-1*pred_glm_test0))
ロジスティック回帰（ランダム効果あり）のあてはめ。（GLMM）
説明変数はSexとAge、ランダム効果にSales_IDを指定。
library(glmmML)
res_glmm1 <- glmmML(y~Sex+Age, data=data1, family=binomial, cluster=Sales_ID)
Call: glmmML(formula = y ~ Sex + Age, family = binomial, data = data1, cluster = Sales_ID)
coef se(coef) z Pr(>|z|)
(Intercept) -3.1060 0.5724 -5.4267 5.74e-08
SexF 1.4090 0.2760 5.1057 3.30e-07
Age-60 0.1904 0.3411 0.5583 5.77e-01
Age60- 1.0691 0.3137 3.4079 6.55e-04
Scale parameter in mixing distribution: 1.939 gaussian
Std. Error: 0.276
LR p-value for H_0: sigma = 0: 8.066e-131
Residual deviance: 671.3 on 995 degrees of freedom AIC: 681.3
真のパラメータはInterceptが-4、Sex Fが2、Age -60が0.5、Age 60-が1.5と設定してあるの。まあ、上のGLMよりかは近い気がする。また、sは1.939と推定されており、実際の2.163217（2と設定して、ランダムに20サンプルの正規確率変数を生成し、その標準偏差が2.163217）に近い。
predict関数はglmmMLに未対応のよう。
predict(res_glmm1)
Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "glmmML"
仕方がないので、自作予測関数を定義する。今回は変量効果の予測値を、予測に含めない。（各インサイドセールスの違いの効果を予測に含めない。）
推定パラメータを渡す。
para_intercept <- res_glmm1$coefficients[1]
para_Sex_F <- res_glmm1$coefficients[2]
para_Age_under60 <- res_glmm1$coefficients[3]
para_Age_over60 <- res_glmm1$coefficients[4]
説明変数のデータフレーム（Sex, Age）を引数に取り、予測qを返す関数。
fn_predict <- function(vec){
s <- vec[1]
a <- vec[2]
if(s=="F"){eff_sex <- para_Sex_F}
else if(s=="M"){eff_sex <- 0}
else{stop("Sex should be F or M")}
if(a=="-40"){eff_age <- 0}
else if(a=="-60"){eff_age <- para_Age_under60}
else if(a=="60-"){eff_age <- para_Age_over60}
else{stop("Age should be -40, -60, or 60-")}
lin_fn <- para_intercept + eff_sex + eff_age
return( 1 / (1+exp(-lin_fn)) )
}
Trainデータへ作成したモデルをあてはめ、qを出力。
pred_glmm <- apply(data1[4:5], 1, fn_predict)
Testデータへ作成したモデルをあてはめ、qを出力。
pred_glmm_test <- apply(data1_test[4:5], 1, fn_predict)
データフレームにまとめる。
data2 <- rbind( cbind(data1, pred_glm, pred_glmm), cbind(data1_test, pred_glm=pred_glm_test, pred_glmm=pred_glmm_test) )
Row_ID Sales_ID SalesRandEffect Sex Age q y Test
1 1 1 -2.3580976 M 60- 0.0077054077 0 Train
2 2 1 -2.3580976 F -60 0.0206717743 0 Train
3 3 1 -2.3580976 M -40 0.0017296628 0 Train
4 4 1 -2.3580976 F 60- 0.0542642460 0 Train
5 5 1 -2.3580976 F 60- 0.0542642460 0 Train
6 6 1 -2.3580976 F 60- 0.0542642460 0 Train
7 7 1 -2.3580976 M 60- 0.0077054077 0 Train
8 8 1 -2.3580976 F 60- 0.0542642460 0 Train
9 9 1 -2.3580976 M 60- 0.0077054077 0 Train
10 10 1 -2.3580976 F 60- 0.0542642460 0 Train
...................................
pred_glm pred_glmm
1 0.2799215 0.11538247
2 0.3451238 0.18144787
3 0.2122954 0.04285945
4 0.4148172 0.34798986
5 0.4148172 0.34798986
6 0.4148172 0.34798986
7 0.2799215 0.11538247
8 0.4148172 0.34798986
9 0.2799215 0.11538247
10 0.4148172 0.34798986
....................................
結果を集計してみる。
GLMの方が、気持ちパフォーマンスが良いように見える。ただし、上記したように、変量効果（インサイドセールス）の予測値をこの予測結果に含めてはいない。
glmmMLはprediction関数をサポートしてない。もともと医療系の統計っぽいことをやってきてた私自身、GLMMを予測に使うイメージがあまりなかった。
RでGLMMを実行する場合、glmmMLでなくlme4パッケージの方が標準だと聞いたこともある。
lme4のモデル関数（lmer, glmer）はpredict関数をサポートしているようである。
lme4: Mixed-effects modeling with R
今後、lme4の使用メモを残す予定。
以上。
