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

東京に棲む日々

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

GLMMを勉強してみる。 - R glmmML

GLMM(リンク関数=ロジスティック関数、分布=二項分布)を勉強してみる。RのglmmML使用。

 

久保(2012) の7章を参考。

 

ロジスティック回帰モデルの復習。

f:id:High_School_Student:20150418134359j:plain

 

これにランダム効果を加えると。

 

f:id:High_School_Student:20150418134410j:plain

 


文献で取り上げている種子の例でなく、自分の事例に合ったサンプルデータを作成してみる。

「インサイドセールス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でグラフを描く。

f:id:High_School_Student:20150419075911j:plain

上段が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
....................................

 

結果を集計してみる。 

f:id:High_School_Student:20150419082559j:plain

 

GLMの方が、気持ちパフォーマンスが良いように見える。ただし、上記したように、変量効果(インサイドセールス)の予測値をこの予測結果に含めてはいない。

 

glmmMLはprediction関数をサポートしてない。もともと医療系の統計っぽいことをやってきてた私自身、GLMMを予測に使うイメージがあまりなかった。

 

RでGLMMを実行する場合、glmmMLでなくlme4パッケージの方が標準だと聞いたこともある。

lme4のモデル関数(lmer, glmer)はpredict関数をサポートしているようである。 

 

Package ‘lme4’

lme4: Mixed-effects modeling with R

 

今後、lme4の使用メモを残す予定。

以上。