霞と側杖を食らう

ほしいものです。なにかいただけるとしあわせです。[https://www.amazon.jp/hz/wishlist/ls/2EIEFTV4IKSIJ?ref_=wl_share]

ロジットモデルの人工データ生成とOVBの学習記録

 

【学習動機】

傾向スコアをそのうち仕事で使いそうな雰囲気があるので、傾向スコアを学び直しておこう、というかちゃんと学ぼうってことになったわけなんだが、その前にロジットモデルを学び直しておこうと思った。そんなこんなで、Trainの“Discrete Choice Methods with Simulation”をパラパラと読んでいる。

https://www.amazon.co.jp/Discrete-Choice-Methods-Simulation-Kenneth/dp/0521747384

ロジットモデルを学んですぐのチュートリアルでは、実データでやることが多いが、それらは本当に推定できているのかという疑問を持った。
今回はシミュレーションデータを作りながら、ちゃんと推定できるのか確認していき、学び直す。ついでに、ロジットモデルでもOVB(Omitted Variable Bias : 欠落変数バイアス)が生じることをシミュレーションで確認する。

 

【学習内容】

ロジットモデルに関して、ざっくりとメモ。離散選択モデルはランダム効用モデルから導出される。個人は選択によって効用が最大になるような選択肢を選ぶものとし、効用(U)を代表的効用(V)と誤差項を分解し、誤差項にガンベル分布を当てはめると、選択肢iの選択確率がexp(Vi)/sum(exp(Vi))となってロジットモデルが導かれる。Vを線形モデルで当てはめ、選択確率の積に対数をとって対数尤度関数を導出し、最尤法でもって、パラメータを推定する。対数尤度関数は線形のパラメータに関して大局的に凹であることはMcFadden(1974)で証明されているようだ(標準的なロジットモデルではなく、ネスティッドロジットモデルになるとそうではないとTrainのテキストのGEVの章に書かれている)。

 

【追加 : 2019/11/09】

欠落変数バイアスのシミュレーションで良い資料を見つけたのでリンクを貼っておく。

yukiyanai.github.io

 

 

以下でロジットモデルの人工データを生成して、ロジットモデルで推していく。

ライブラリ読み込み。

library(ggplot2)    #visualization

関数を先に定義しておく。

#ロジット関数
logit <- function(p){
  return( log(p) - log(1-p) )
}
#ロジスティク関数
logistic <- function(y){
  return( 1 / (1 + exp(-y)))
}

モデルの説明

真のモデルのパラメータ設定\(\beta_{1}=-4, \beta_{2}=8\)でサンプルサイズは1000, シミュレーション回数は1000としておく。

# true parameter
beta1 <- -4
beta2 <- 8
# sample size
n <- 1000
# # of simulation
m <- 1000
#seed
set.seed(1234)

説明変数x1とx2を作っておく。図示もしておく。

x1 <- c(runif(n = n/2, min = 0, max = 800), runif(n = n/2, min = 400, max = 1600))/1000
x2 <- rbinom(n = n, size = 1, prob = 0.6)
plot(x1)

plot(x2)

今回の真のモデルについての説明。
選択確率は以下のように、
\(\beta_{1}=-4, \beta_{2}=8\)がパラメータとなり、
\(y^{*} = x1*\beta_{1} + x2*\beta_{2}\)
\(P = \frac{1}{1+exp(-y^{*})}\)とする。
選択は以下のように、
\(y\) ~ \(Bernoulli(P)\)
データが発生するものとする。

このデータに対して、\(y^{*} = x1*\beta_{1} + x2*\beta_{2}\)で真のモデルに等しいモデル(model0と呼ぶ)で推定する場合と、\(y^{*} = x1*\beta_{1}\)でx2が脱落している過少定式のモデル(modelmisと呼ぶ)で推定する場合を見ていく。

err <- evd::rgumbel(n = n, loc = 0, scale = 1)
ystar <- x1*beta1 + x2*beta2
utility <- ystar + err
p <- logistic(ystar)
y <- rbinom(n = n, size = 1, prob = p)
  
df <- data.frame(
    ystar = ystar,
    utility = utility,
    p = p,
    y = y,
    x1 = x1,
    x2 = x2
  )
  
model0 <- glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), data = df)
modelmis <- glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), data = df)

summary(model0)
## 
## Call:
## glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -3.11148  -0.18597   0.06014   0.12358   2.68804  
## 
## Coefficients:
##    Estimate Std. Error z value Pr(>|z|)    
## x1  -3.8564     0.3471  -11.11   <2e-16 ***
## x2   7.8470     0.5388   14.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.29  on 1000  degrees of freedom
## Residual deviance:  323.78  on  998  degrees of freedom
## AIC: 327.78
## 
## Number of Fisher Scoring iterations: 7
summary(modelmis)
## 
## Call:
## glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5064  -1.2777   0.9916   1.0831   1.1771  
## 
## Coefficients:
##    Estimate Std. Error z value Pr(>|z|)    
## x1  0.46689    0.07947   5.875 4.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.3  on 1000  degrees of freedom
## Residual deviance: 1350.5  on  999  degrees of freedom
## AIC: 1352.5
## 
## Number of Fisher Scoring iterations: 4

当然、model0はよく推定できているし、modelmisよりAICが低い。modelmisの\(\beta_{1}\)の推定値が全然違うし、符号も逆という結果になっている。

ちなみに、今回は使わないが、{evd}パッケージのrgumbel()関数でガンベル分布からの乱数を発生させられる。
utility出してるけど、ここでは、なんの意味もない。utilityのmax取る選択肢をとる実験を何回も繰り返すと、選択確率がざっくり出てくるだろうけど。

シミュレーション

上の人工データ生成と推定を1000回繰り返す。model0で推定されたパラメータはb1_0とb2_0に保存し、modelmisで推定されたパラメータはb1_misに保存する。

#for estimates
b1_0 <- numeric(m)
b1_mis <- numeric(m)
b2_0 <- numeric(m)

#seed
set.seed(1234)

for(i in 1:m){
  err <- evd::rgumbel(n = n, loc = 0, scale = 1)
  ystar <- x1*beta1 + x2*beta2
  utility <- ystar + err
  p <- logistic(ystar)
  y <- rbinom(n = n, size = 1, prob = p)
  
  df <- data.frame(
    ystar = ystar,
    utility = utility,
    p = p,
    y = y,
    x1 = x1,
    x2 = x2
  )
  
  model0 <- glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), data = df)
  modelmis <- glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), data = df)
  
  b1_0[i] <- model0$coefficients[[1]]
  b1_mis[i] <- modelmis$coefficients[[1]]
  b2_0[i] <- model0$coefficients[[2]]
}
df_est <- data.frame(
  b1_0 = b1_0,
  b1_mis = b1_mis,
  b2_0 = b2_0
)

#histogram of estimates of beta1
p0 <- ggplot(df_est,aes(x=b1_0)) + geom_histogram(colour="grey") 
pmis <- ggplot(df_est,aes(x=b1_mis)) + geom_histogram(colour="grey") 

#b1
gridExtra::grid.arrange(p0, pmis)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#b2
ggplot(df_est,aes(x=b2_0)) + geom_histogram(colour="grey") 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

model0は真のモデルと定式が同じだからしっかり推定できている。一方、modelmisは先ほど確認した通り、\(\beta_{1}\)の符号から違うという結果に。

 

【追加 : 2019/11/09】

以上のシミュレーションにおいて欠落変数と残っている変数の相関、共分散を考えていなかった。欠落変数バイアスの大きさを考えるなら、そこを考慮すべきだったかもしれない。

 

失敗例

推定が上手くいった例を上に載せたが, 失敗例も役立つかもしれないので載せておく. 失敗例を2つ. 映画のエンディングのNGシーン集的な.

サンプルサイズが100のケース

nを減らすと正しく収束しなかった. サンプルサイズを100として同じシミュレーションをやる

# sample size
n <- 100
set.seed(1234)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredと警告が出る.
推定が上手くいってない. 失敗というか上手くいってない推定が何回か起きてる感じか.

x1を1000で割らないケース

サンプルサイズを元の戻して, x1を1000で割らずに同じシミュレーションをする.

# sample size
n <- 1000
set.seed(1234)
x1 <- c(runif(n = n/2, min = 0, max = 800), runif(n = n/2, min = 400, max = 1600))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

最適化の問題か. 探索範囲とかそういう問題で推定が上手くいってないと思われる. 最適化の実践的な問題に詳しくないので, 事実かどうかはよく分からない.

【学習予定】

次はネスティッドロジットモデルあたりのシミュレーションデータ作りながら確認か, 傾向スコアの方をやってみたい.
ロジットモデルのOVBの理論が簡単に読めるならそれも読んでおきたい.