{mlogit}パッケージ(ネスティッドロジットモデル)の学習記録

【学習動機】

前回の続きで少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.

moratoriamuo.hatenablog.com


Kenneth Train’s exercises using the mlogit package for R
http://www.idg.pl/mirrors/CRAN/web/packages/mlogit/vignettes/Exercises.pdf
これに沿って学習していく. 今回は, ネスティッドロジットモデル(Nested logit model)の節を扱う.
ドキュメントの全てを扱うわけではない. また変数名やコードなど自分の理解しやすい形に変えているところもあるので注意されたい.

Exercise 2: Nested logit model

ググって探索していたら, 二週間前にこれがアップロードされていたことに気付いた. 英語が読めるならこっちを読んだ方が絶対にいい気がする.

加えてこれも発見した. 理論面はTrainのテキストに沿っていて, 補完的に読めてとてもよい.
{mlogit}パッケージ作った人の解説
Yves Croissant “Estimation of multinomial logit models in R : The mlogit Packages”  
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.303.7401&rep=rep1&type=pdf

【学習内容】

2. ネスティッドロジットモデル

ライブラリ読み込み.

library(mlogit)
library(tidyverse)

数式などのメモ書き

ネストkに含まれる選択肢の集合が\(B_k\)と表記する.

ネスティッドロジットモデルでは, \(\epsilon_n = <\epsilon_{n1}, \epsilon_{n2},..., \epsilon_{nJ}>\)は以下の分布関数を持つ.
\(exp(- \Sigma_{k=1}^{K}{(\Sigma_{j\in{B_k}}e^{-\epsilon_{nj}/\lambda_k}})^{ \lambda_{k} })\)

(なんか表示が上手くいかなかったのでスクショ追加)

f:id:moratoriamuo271:20190619222519p:plain

 同じネストに含まれる選択肢\(j,m\in B_{k}\)の誤差項\(\epsilon_{nj}\)\(\epsilon_{nm}\)は相関している.
同じネストに含まれない選択肢\(j\in B_{k}, m\in B_{l}\)の誤差項\(\epsilon_{nj}\)\(\epsilon_{nm}\)は相関しない.

パラメータ\(\lambda_k\)はネストkにおける誤差項の独立の度合いとなっている.
逆に, \(1 - \lambda_k\)は相関の強さの指標になる. McFadden(1978)では実際の相関は\(1 - \lambda_k\)より複雑であるが, 相関の指標としては有用であることが指摘されているようだ.
\(\lambda_k\)が大きいと相関が低い. \(\lambda_k = 1\)のとき, ネスト内で完全に独立になり, 標準的なロジットモデルに等しくなる.

\(P_{ni} = \frac{e^{V_{ni}/ \lambda_k} (\Sigma_{j \in B_k}{e^{V_{nj}/\lambda_k } )^{\lambda_k - 1}}} {\Sigma_{l=1}^{K}{(\Sigma_{j \in B_l}{e^{V_{ni}/\lambda_l})^{\lambda_l}}}}\)

f:id:moratoriamuo271:20190619222602p:plain

 

 (なんか表示が上手くいかなかったのでスクショ追加(再度))

データの確認

カリフォルニア州の250家計の冷暖房システムの選択データ.
新築でエアコンの設置された戸建て住宅.

data("HC", package = "mlogit")
head(HC)
##   depvar ich.gcc ich.ecc ich.erc ich.hpc ich.gc ich.ec ich.er  icca
## 1    erc    9.70    7.86    8.79   11.36  24.08  24.50   7.37 27.28
## 2    hpc    8.77    8.69    7.09    9.37  28.00  32.71   9.33 26.49
## 3    gcc    7.43    8.86    6.94   11.70  25.71  31.68   8.14 22.63
## 4    gcc    9.18    8.93    7.22   12.13  29.72  26.73   8.04 25.33
## 5    gcc    8.05    7.02    8.44   10.51  23.90  28.35   7.15 25.45
## 6    gcc    9.32    8.03    6.22   12.57  27.02  21.37   8.60 19.93
##   och.gcc och.ecc och.erc och.hpc och.gc och.ec och.er occa income
## 1    2.26    4.09    3.85    1.73   2.26   4.09   3.85 2.95     20
## 2    2.30    2.69    3.45    1.65   2.30   2.69   3.45 1.63     50
## 3    2.28    5.25    4.35    1.44   2.28   5.25   4.35 2.18     50
## 4    2.62    4.89    4.85    1.93   2.62   4.89   4.85 2.70     50
## 5    2.52    3.71    3.64    1.63   2.52   3.71   3.64 2.77     60
## 6    1.99    3.20    4.30    1.30   1.99   3.20   4.30 2.68     30

選択肢は7つ:
gcc (gas central heat with cooling),
ecc (electric central resistence heat with cooling),
erc (electric room resistence heat with cooling),
hpc (electric heat pump which provides cooling also),
gc (gas central heat without cooling),
ec (electric central resistence heat without cooling),
er (electric room resistence heat without cooling)
ヒートポンプは低温部分から高温部分に熱を移動させるもの.

変数
depvar:選択された選択肢
ich.alt:暖房の設置コスト
icca:冷房の設置コスト
och.alt:暖房の運用コスト
occa:冷房の運用コスト
income:世帯の年間収入

選択の割合を確認

table(HC$depvar)/250    #250人で割って割合を出す
## 
##   gcc   ecc   erc   hpc    gc    ec    er 
## 0.744 0.016 0.004 0.104 0.096 0.004 0.032

gccが74%, hpcが10%, gcが9%

wide形式から変換

HC <- mlogit.data(HC, varying = c(2:8, 10:16), choice = "depvar", shape = "wide")
head(HC)
##       depvar  icca occa income alt   ich  och chid
## 1.ec   FALSE 27.28 2.95     20  ec 24.50 4.09    1
## 1.ecc  FALSE 27.28 2.95     20 ecc  7.86 4.09    1
## 1.er   FALSE 27.28 2.95     20  er  7.37 3.85    1
## 1.erc   TRUE 27.28 2.95     20 erc  8.79 3.85    1
## 1.gc   FALSE 27.28 2.95     20  gc 24.08 2.26    1
## 1.gcc  FALSE 27.28 2.95     20 gcc  9.70 2.26    1

各家計の選択関係なく平均

HC %>% group_by(alt) %>% summarise(ave_income = mean(income),
                                   ave_icca = mean(icca),
                                   ave_occa = mean(occa),
                                   ave_ich = mean(ich),
                                   ave_och = mean(och))
## # A tibble: 7 x 6
##   alt   ave_income ave_icca ave_occa ave_ich ave_och
##   <chr>      <dbl>    <dbl>    <dbl>   <dbl>   <dbl>
## 1 ec          44.7     25.4     2.49   25.5     4.08
## 2 ecc         44.7     25.4     2.49    8.35    4.08
## 3 er          44.7     25.4     2.49    7.76    4.13
## 4 erc         44.7     25.4     2.49    7.89    4.13
## 5 gc          44.7     25.4     2.49   25.0     2.25
## 6 gcc         44.7     25.4     2.49    7.82    2.25
## 7 hpc         44.7     25.4     2.49   10.6     1.51

選択を条件付けて平均

HC %>% filter(depvar == TRUE) %>% group_by(alt) %>% summarise(ave_income = mean(income),
                                                              ave_icca = mean(icca),
                                                              ave_occa = mean(occa),
                                                              ave_ich = mean(ich),
                                                              ave_och = mean(och))
## # A tibble: 7 x 6
##   alt   ave_income ave_icca ave_occa ave_ich ave_och
##   <chr>      <dbl>    <dbl>    <dbl>   <dbl>   <dbl>
## 1 ec          30       26.1     2.57   18.3     3.18
## 2 ecc         45       25.4     2.42    7.73    3.24
## 3 er          20       25.0     2.61    8.39    3.93
## 4 erc         20       27.3     2.95    8.79    3.85
## 5 gc          24.2     25.8     2.51   21.5     2.20
## 6 gcc         47.8     25.4     2.48    7.79    2.22
## 7 hpc         50.4     24.8     2.47    9.79    1.50

モデルのあてはめ

準備

#coolのありなしindex:三文字目にcがついてるとcoolingあり
cooling.modes <- index(HC)$alt %in% c("gcc", "ecc", "erc", "hpc")
#roomのありなしindex:二文字目がrだとroom(中央制御と部屋ごとってこと?よく分からない.)
room.modes <- index(HC)$alt %in% c("erc", "er")

#冷房あり(gcc, ecc, erc, hpc)と冷房なし (gc, ec, er)でネストを分ける
#冷房じゃないケースの冷房設置コストと冷房運用コストを0に置き換える
HC$icca[!cooling.modes] <- 0
HC$occa[!cooling.modes] <- 0

#冷房用のincomeと部屋用のincomeの変数を作って
HC$inc.cooling <- 0
HC$inc.room <- 0
HC$inc.cooling[cooling.modes] <- HC$income[cooling.modes]
HC$inc.room[room.modes] <- HC$income[room.modes]

HC$int.cooling <- as.numeric(cooling.modes)

head(HC)
##       depvar  icca occa income alt   ich  och chid inc.cooling inc.room
## 1.ec   FALSE  0.00 0.00     20  ec 24.50 4.09    1           0        0
## 1.ecc  FALSE 27.28 2.95     20 ecc  7.86 4.09    1          20        0
## 1.er   FALSE  0.00 0.00     20  er  7.37 3.85    1           0       20
## 1.erc   TRUE 27.28 2.95     20 erc  8.79 3.85    1          20       20
## 1.gc   FALSE  0.00 0.00     20  gc 24.08 2.26    1           0        0
## 1.gcc  FALSE 27.28 2.95     20 gcc  9.70 2.26    1          20        0
##       int.cooling
## 1.ec            0
## 1.ecc           1
## 1.er            0
## 1.erc           1
## 1.gc            0
## 1.gcc           1

mlogit()でネスティッドロジットモデルするときの引数:
nestsにネストのリストを入れる.
un.nest.el:TRUEにするとネストごとに違うλにしないで唯一のλを与える.
返り値の係数のivの変数がλにあたる. λはlog-sum coefficientとも呼ばれる.

ネスト2つのネスティッドロジットモデル(ネスト内の相関がネスト間で等しい)

各設置コストと各運用コスト, 部屋ありの収入と冷房ありの収入と冷房の有無で説明.
ネストは冷房ありネストと冷房なしネスト.

model_nl <- mlogit(depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 0, data = HC,
             nests = list(cooling = c("gcc", "ecc","erc", "hpc"), other = c("gc", "ec", "er")),
             un.nest.el = TRUE)
summary(model_nl)
## 
## Call:
## mlogit(formula = depvar ~ ich + och + icca + occa + inc.room + 
##     inc.cooling + int.cooling | 0, data = HC, nests = list(cooling = c("gcc", 
##     "ecc", "erc", "hpc"), other = c("gc", "ec", "er")), un.nest.el = TRUE)
## 
## Frequencies of alternatives:
##    ec   ecc    er   erc    gc   gcc   hpc 
## 0.004 0.016 0.032 0.004 0.096 0.744 0.104 
## 
## bfgs method
## 11 iterations, 0h:0m:0s 
## g'(-H)^-1g = 7.26E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##              Estimate Std. Error z-value  Pr(>|z|)    
## ich         -0.554878   0.144205 -3.8478 0.0001192 ***
## och         -0.857886   0.255313 -3.3601 0.0007791 ***
## icca        -0.225079   0.144423 -1.5585 0.1191212    
## occa        -1.089458   1.219821 -0.8931 0.3717882    
## inc.room    -0.378971   0.099631 -3.8038 0.0001425 ***
## inc.cooling  0.249575   0.059213  4.2149 2.499e-05 ***
## int.cooling -6.000415   5.562423 -1.0787 0.2807030    
## iv           0.585922   0.179708  3.2604 0.0011125 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -178.12

log-sum coefficient = 1の帰無仮説をt検定

(coef(model_nl)['iv'] - 1) / sqrt(vcov(model_nl)['iv', 'iv'])
##        iv 
## -2.304171

帰無仮説は棄却された. ネスト内の相関はありで良さそう.

ネストなしの多項ロジットモデル

model_ml <- update(model_nl, nests = NULL)
summary(model_ml)
## 
## Call:
## mlogit(formula = depvar ~ ich + och + icca + occa + inc.room + 
##     inc.cooling + int.cooling | 0, data = HC, un.nest.el = TRUE, 
##     method = "nr")
## 
## Frequencies of alternatives:
##    ec   ecc    er   erc    gc   gcc   hpc 
## 0.004 0.016 0.032 0.004 0.096 0.744 0.104 
## 
## nr method
## 9 iterations, 0h:0m:0s 
## g'(-H)^-1g = 0.000276 
## successive function values within tolerance limits 
## 
## Coefficients :
##               Estimate Std. Error  z-value  Pr(>|z|)    
## ich          -0.851583   0.078790 -10.8083 < 2.2e-16 ***
## och          -1.356336   0.147398  -9.2019 < 2.2e-16 ***
## icca         -0.257236   0.126969  -2.0260   0.04277 *  
## occa         -1.413791   1.149068  -1.2304   0.21855    
## inc.room     -0.580337   0.063257  -9.1743 < 2.2e-16 ***
## inc.cooling   0.314117   0.053994   5.8177 5.967e-09 ***
## int.cooling -10.628463   5.129315  -2.0721   0.03826 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -180.29

ネストありとなしで, 尤度比検定

lrtest(model_nl,model_ml)
## Likelihood ratio test
## 
## Model 1: depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 
##     0
## Model 2: depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 
##     0
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1   8 -178.12                       
## 2   7 -180.29 -1 4.3234    0.03759 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ネストありで良さそう.

ネスト2つのネスティッドロジットモデル(ネスト内の相関がネスト間で異なる)

un.nest.el = FALSEにする. つまり, λをそれぞれ別のものにする.

model_nl2 <- mlogit(depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 0, data = HC,
                    nests = list(cooling = c("gcc", "ecc","erc", "hpc"), other = c("gc", "ec", "er")),
                    un.nest.el = FALSE)
summary(model_nl2)
## 
## Call:
## mlogit(formula = depvar ~ ich + och + icca + occa + inc.room + 
##     inc.cooling + int.cooling | 0, data = HC, nests = list(cooling = c("gcc", 
##     "ecc", "erc", "hpc"), other = c("gc", "ec", "er")), un.nest.el = FALSE)
## 
## Frequencies of alternatives:
##    ec   ecc    er   erc    gc   gcc   hpc 
## 0.004 0.016 0.032 0.004 0.096 0.744 0.104 
## 
## bfgs method
## 4 iterations, 0h:0m:0s 
## g'(-H)^-1g =  1.18 
## last step couldn't find higher value 
## 
## Coefficients :
##              Estimate Std. Error z-value  Pr(>|z|)    
## ich         -0.562283   0.146145 -3.8474 0.0001194 ***
## och         -0.895493   0.271861 -3.2939 0.0009880 ***
## icca        -0.267062   0.150310 -1.7767 0.0756103 .  
## occa        -1.338514   1.264215 -1.0588 0.2897042    
## inc.room    -0.381441   0.096658 -3.9463 7.937e-05 ***
## inc.cooling  0.259932   0.062085  4.1867 2.830e-05 ***
## int.cooling -4.821927   5.528796 -0.8721 0.3831277    
## iv:cooling   0.611529   0.188736  3.2401 0.0011947 ** 
## iv:other     0.378394   0.133617  2.8319 0.0046270 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -178.04

ivの値がそれぞれ0.61と0.37. 相関の指標は0.39と0.63.

lrtest(model_nl, model_nl2)
## Likelihood ratio test
## 
## Model 1: depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 
##     0
## Model 2: depvar ~ ich + och + icca + occa + inc.room + inc.cooling + int.cooling | 
##     0
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -178.12                     
## 2   9 -178.04  1 0.1758      0.675

2つのネストのlog-sum coefficient が等しいという帰無仮説を棄却できない.

ネスト3つのネスティッドロジットモデル(ネスト内の相関がネスト間で等しい)

model_nl3 <- update(model_nl, nests=list(n1 = c('gcc', 'ecc', 'erc'), n2 = c('hpc'),
                                         n3 = c('gc', 'ec', 'er')))
summary(model_nl3)
## 
## Call:
## mlogit(formula = depvar ~ ich + och + icca + occa + inc.room + 
##     inc.cooling + int.cooling | 0, data = HC, nests = list(n1 = c("gcc", 
##     "ecc", "erc"), n2 = c("hpc"), n3 = c("gc", "ec", "er")), 
##     un.nest.el = TRUE)
## 
## Frequencies of alternatives:
##    ec   ecc    er   erc    gc   gcc   hpc 
## 0.004 0.016 0.032 0.004 0.096 0.744 0.104 
## 
## bfgs method
## 8 iterations, 0h:0m:0s 
## g'(-H)^-1g = 3.71E-08 
## gradient close to zero 
## 
## Coefficients :
##               Estimate Std. Error z-value  Pr(>|z|)    
## ich          -0.838394   0.100546 -8.3384 < 2.2e-16 ***
## och          -1.331598   0.252069 -5.2827 1.273e-07 ***
## icca         -0.256131   0.145564 -1.7596   0.07848 .  
## occa         -1.405656   1.207281 -1.1643   0.24430    
## inc.room     -0.571352   0.077950 -7.3297 2.307e-13 ***
## inc.cooling   0.311355   0.056357  5.5247 3.301e-08 ***
## int.cooling -10.413384   5.612445 -1.8554   0.06354 .  
## iv            0.956544   0.180722  5.2929 1.204e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -180.26

対数尤度が上がらない.

ネスト2つのネスティッドロジットモデル(ネスト内の相関がネスト間で等しい)のが良さそう.

【学習予定】

IIA特性の話はよく見かけるけどネスティッドロジットモデルやってみたよっていう日本語ブログはあまり見かけない気がする. 

それはそれとして, 次は, ミックスドロジットモデル(Mixed logit model)の節を扱うかも.  

{mlogit}パッケージ(多項ロジットモデル)の学習記録

 

【学習動機】

少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.
調べていたら, Train先生が{mlogit}の使い方のエクササイズを書いているのを発見した.
Kenneth Train’s exercises using the mlogit package for R
http://www.idg.pl/mirrors/CRAN/web/packages/mlogit/vignettes/Exercises.pdf
これに沿って学習していく. 今回は, 多項ロジットモデルの節を扱う.
ドキュメントの全てを扱うわけではない. また変数名やコードなど自分の理解しやすい形に変えているところもあるので注意されたい.  

【追記:2019/06/16】  

Exercise 1: Multinomial logit model  

ググって探索していたら, 二週間前にこれがアップロードされていたことに気付いた. 英語が読めるならこっちを読んだ方が絶対にいい気がする.  

 

【学習内容】

1. 多項ロジットモデル

library(mlogit)
library(tidyverse)
ライブラリ読み込み 
 

カリフォルニア州の900家計の暖房システムの選択データ
新築でエアコンの設置された戸建て住宅
選択肢は5つ:gc (gas central), gr (gas room), ec (electric central), er (electric room), hp (heat pump)

data("Heating",package = "mlogit")
head(Heating)
##   idcase depvar  ic.gc  ic.gr  ic.ec  ic.er   ic.hp  oc.gc  oc.gr  oc.ec
## 1      1     gc 866.00 962.64 859.90 995.76 1135.50 199.69 151.72 553.34
## 2      2     gc 727.93 758.89 796.82 894.69  968.90 168.66 168.66 520.24
## 3      3     gc 599.48 783.05 719.86 900.11 1048.30 165.58 137.80 439.06
## 4      4     er 835.17 793.06 761.25 831.04 1048.70 180.88 147.14 483.00
## 5      5     er 755.59 846.29 858.86 985.64  883.05 174.91 138.90 404.41
## 6      6     gc 666.11 841.71 693.74 862.56  859.18 135.67 140.97 398.22
##    oc.er  oc.hp income agehed rooms region
## 1 505.60 237.88      7     25     6 ncostl
## 2 486.49 199.19      5     60     5 scostl
## 3 404.74 171.47      4     65     2 ncostl
## 4 425.22 222.95      2     50     4 scostl
## 5 389.52 178.49      2     25     6 valley
## 6 371.04 209.27      6     65     7 scostl

変数は以下のようになっている.
idcase:id
depvar:暖房システム
ic.alt:暖房システムaltの設置コスト
oc.alt:暖房システムaltの年間運用コスト
pb.alt:比率oc.alt / ic.alt
income:世帯の年収
agehed:世帯主の年齢
rooms:部屋の数
region:ncostl(北部沿岸地域)、scostl(南部沿岸地域)、mountn(山岳地帯)、valley(中央渓谷地域)

選択の割合を確認

table(Heating$depvar)/900    #900人で割って割合を出す
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

gc : 64%, gr : 14%, ec : 7%, er : 9%, hp : 6%

wide形式のデータになっているので, mlogit.data()で{mlogit}パッケージで推定できるように形式を変更.

H_long <- mlogit.data(Heating, shape = "wide", choice = "depvar",
                       varying = c(3:12))

各家計の選択関係なく平均

H_long %>% group_by(alt) %>% summarise(ave_income = mean(income),
                                       ave_age = mean(agehed),
                                       ave_rooms = mean(rooms),
                                       ave_ic = mean(ic),
                                       ave_oc = mean(oc))
## # A tibble: 5 x 6
##   alt   ave_income ave_age ave_rooms ave_ic ave_oc
##   <chr>      <dbl>   <dbl>     <dbl>  <dbl>  <dbl>
## 1 ec          4.64    42.9      4.42   825.   477.
## 2 er          4.64    42.9      4.42   984.   430.
## 3 gc          4.64    42.9      4.42   777.   172.
## 4 gr          4.64    42.9      4.42   922.   154.
## 5 hp          4.64    42.9      4.42  1046.   219.

選択を条件付けて平均

H_long %>% filter(depvar == TRUE) %>% 
  group_by(alt) %>% summarise(ave_income = mean(income),
                              ave_age = mean(agehed),
                              ave_rooms = mean(rooms),
                              ave_ic = mean(ic),
                              ave_oc = mean(oc))
## # A tibble: 5 x 6
##   alt   ave_income ave_age ave_rooms ave_ic ave_oc
##   <chr>      <dbl>   <dbl>     <dbl>  <dbl>  <dbl>
## 1 ec          4.70    44.3      4.59   789.   451.
## 2 er          4.62    39.3      4.52   928.   407.
## 3 gc          4.68    43.3      4.40   781.   173.
## 4 gr          4.36    43.9      4.37   933.   155.
## 5 hp          4.88    40.5      4.46  1027.   213.

ベースモデル

設置コストと年間運用コストで説明するモデル. コストなので効用を下げるはずなので, 係数は負になると予想される.

model_0 <- mlogit(depvar ~ ic + oc | 0, data = H_long)
summary(model_0)
## 
## Call:
## mlogit(formula = depvar ~ ic + oc | 0, data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 4 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.56E-07 
## gradient close to zero 
## 
## Coefficients :
##       Estimate  Std. Error z-value  Pr(>|z|)    
## ic -0.00623187  0.00035277 -17.665 < 2.2e-16 ***
## oc -0.00458008  0.00032216 -14.217 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1095.2

全選択肢の確率

model_0 %>% fitted(outcome=FALSE) %>% head()
##           ec         er        gc        gr         hp
## 1 0.09545811 0.05094155 0.4642482 0.3166757 0.07267644
## 2 0.05831236 0.03698376 0.4482688 0.3696122 0.08682283
## 3 0.08396029 0.03195210 0.6220979 0.2250557 0.03693402
## 4 0.11457049 0.09663155 0.2883760 0.4375632 0.06285874
## 5 0.07692237 0.03737185 0.4188408 0.2806727 0.18619232
## 6 0.13353597 0.05281490 0.5279799 0.1725121 0.11315710

900人の平均選択確率

model_0 %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         ec         er         gc         gr         hp 
## 0.10413057 0.05141477 0.51695653 0.24030898 0.08718915

データの割合と比較

table(Heating$depvar)/900
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

フィット的に見て, あまりいいモデルとは言えなそう.

coef(model_0)["oc"]/coef(model_0)["ic"]
##        oc 
## 0.7349453

年間運用コストを1ドル削減するのに支払える設置コストを0.73ドル多く支払えるという結果になっている.
年間運用コストは毎年払うものであるのに対して, 設置コストは1回のみ. この結果は合理的でない.
やはりモデルが良くなさそう.

生涯コストモデル

割引率を使って年間運用コストと設置コストを生涯コストに変換したモデル.

#H_long <- H_long %>% mutate(lcc = ic + oc/0.12)
#ここでdplyr::mutateを使うとmlogit.data型が消えてしまうので使えない
H_long$lcc <- H_long$ic + H_long$oc/0.12
model_lcc <- mlogit(depvar ~ lcc | 0, data = H_long)
summary(model_lcc)
## 
## Call:
## mlogit(formula = depvar ~ lcc | 0, data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 9.32E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##        Estimate  Std. Error z-value  Pr(>|z|)    
## lcc -7.1585e-04  4.2761e-05 -16.741 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1248.7

ベースモデルとの比較で, 尤度比検定を行う.

lrtest(model_0,model_lcc)
## Likelihood ratio test
## 
## Model 1: depvar ~ ic + oc | 0
## Model 2: depvar ~ lcc | 0
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -1095.2                         
## 2   1 -1248.7 -1 306.93  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#年間運用コストの係数と設置コストの係数との関係が制約条件となっている
qchisq(0.05, df = 1, lower.tail = FALSE)
## [1] 3.841459

model_0が支持される.

選択肢別切片モデル

各選択ごとに異なる切片を設定したモデル. hpの選択肢を0として基準にしている.

model_con <- mlogit(depvar ~ ic + oc, data = H_long, reflevel = "hp")
summary(model_con)
## 
## Call:
## mlogit(formula = depvar ~ ic + oc, data = H_long, reflevel = "hp", 
##     method = "nr")
## 
## Frequencies of alternatives:
##       hp       ec       er       gc       gr 
## 0.055556 0.071111 0.093333 0.636667 0.143333 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 9.58E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##                   Estimate  Std. Error z-value  Pr(>|z|)    
## ec:(intercept)  1.65884594  0.44841936  3.6993 0.0002162 ***
## er:(intercept)  1.85343697  0.36195509  5.1206 3.045e-07 ***
## gc:(intercept)  1.71097930  0.22674214  7.5459 4.485e-14 ***
## gr:(intercept)  0.30826328  0.20659222  1.4921 0.1356640    
## ic             -0.00153315  0.00062086 -2.4694 0.0135333 *  
## oc             -0.00699637  0.00155408 -4.5019 6.734e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1008.2
## McFadden R^2:  0.013691 
## Likelihood ratio test : chisq = 27.99 (p.value = 8.3572e-07)

このモデルから得られる900人の平均選択確率

model_con %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         hp         ec         er         gc         gr 
## 0.05555556 0.07111111 0.09333333 0.63666667 0.14333333

データの割合と比較

table(Heating$depvar)/900
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

データでの割合とモデルで予測される平均選択確率が一致していてフィットは良さそう.

支払い意思額

coef(model_con)["oc"]/coef(model_con)["ic"]
##       oc 
## 4.563385

割引率

1/(coef(model_con)["oc"]/coef(model_con)["ic"])
##        oc 
## 0.2191356

割引率が22%と計算でき, 前のモデルに比べたら妥当.

収入あたり設置コストモデル

設置コストを収入あたりにしたモデル.

model_income <- mlogit(depvar ~ oc + I(ic/income), data = H_long)
summary(model_income)
## 
## Call:
## mlogit(formula = depvar ~ oc + I(ic/income), data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.03E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##                  Estimate Std. Error z-value  Pr(>|z|)    
## er:(intercept)  0.0639934  0.1944893  0.3290  0.742131    
## gc:(intercept)  0.0563481  0.4650251  0.1212  0.903555    
## gr:(intercept) -1.4653063  0.5033845 -2.9109  0.003604 ** 
## hp:(intercept) -1.8700773  0.4364248 -4.2850 1.827e-05 ***
## oc             -0.0071066  0.0015518 -4.5797 4.657e-06 ***
## I(ic/income)   -0.0027658  0.0018944 -1.4600  0.144298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1010.2
## McFadden R^2:  0.011765 
## Likelihood ratio test : chisq = 24.052 (p.value = 5.9854e-06)

対数尤度が選択肢別切片モデルより低くなっている.

 

model_income2 <- mlogit(depvar ~ oc + ic | income, data = H_long, reflevel = "hp")
summary(model_income2)
## 
## Call:
## mlogit(formula = depvar ~ oc + ic | income, data = H_long, reflevel = "hp", 
##     method = "nr")
## 
## Frequencies of alternatives:
##       hp       ec       er       gc       gr 
## 0.055556 0.071111 0.093333 0.636667 0.143333 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 6.27E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##                   Estimate  Std. Error z-value  Pr(>|z|)    
## ec:(intercept)  1.95445797  0.70353833  2.7780 0.0054688 ** 
## er:(intercept)  2.30560852  0.62390478  3.6954 0.0002195 ***
## gc:(intercept)  2.05517018  0.48639682  4.2253 2.386e-05 ***
## gr:(intercept)  1.14158139  0.51828845  2.2026 0.0276231 *  
## oc             -0.00696000  0.00155383 -4.4792 7.491e-06 ***
## ic             -0.00153534  0.00062251 -2.4664 0.0136486 *  
## ec:income      -0.06362917  0.11329865 -0.5616 0.5743846    
## er:income      -0.09685787  0.10755423 -0.9005 0.3678281    
## gc:income      -0.07178917  0.08878777 -0.8085 0.4187752    
## gr:income      -0.17981159  0.10012691 -1.7958 0.0725205 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1005.9
## McFadden R^2:  0.01598 
## Likelihood ratio test : chisq = 32.67 (p.value = 1.2134e-05)
#尤度比検定
lrtest(model_con, model_income2)
## Likelihood ratio test
## 
## Model 1: depvar ~ ic + oc
## Model 2: depvar ~ oc + ic | income
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   6 -1008.2                     
## 2  10 -1005.9  4 4.6803     0.3217
#ワルド検定
waldtest(model_con, model_income2)
## Wald test
## 
## Model 1: depvar ~ ic + oc
## Model 2: depvar ~ oc + ic | income
##   Res.Df Df  Chisq Pr(>Chisq)
## 1    894                     
## 2    890  4 4.6456     0.3256
#スコア検定
scoretest(model_con, model_income2)
## 
##  score test
## 
## data:  depvar ~ oc + ic | income
## chisq = 4.6761, df = 4, p-value = 0.3222
## alternative hypothesis: unconstrained model

収入は選択に影響しないのかもしれない.

ヒートポンプの設置コスト割引するケース

hp(ヒートポンプ)の設置費用の選択に対して払い戻しを10%行うとき, シェアはどのようになるか予測する.

X <- model.matrix(model_con)
head(X)
##      ec:(intercept) er:(intercept) gc:(intercept) gr:(intercept)      ic
## 1.ec              1              0              0              0  859.90
## 1.er              0              1              0              0  995.76
## 1.gc              0              0              1              0  866.00
## 1.gr              0              0              0              1  962.64
## 1.hp              0              0              0              0 1135.50
## 2.ec              1              0              0              0  796.82
##          oc
## 1.ec 553.34
## 1.er 505.60
## 1.gc 199.69
## 1.gr 151.72
## 1.hp 237.88
## 2.ec 520.24
alt <- index(H_long)$alt
chid <- index(H_long)$chid
#hpの設置コストを1割引に
X[alt == "hp", "ic"] <- X[alt == "hp", "ic"] * 0.9
eXb <- as.numeric(exp(X %*% coef(model_con)))
SeXb <- tapply(eXb, chid, sum)
P <- eXb/SeXb[chid]
P <- matrix(P, ncol = 5, byrow = TRUE)
P_pred <- apply(P, 2, mean) %>% t()
colnames(P_pred) <- c("ec","er","gc","gr","hp")
P_pred 
##              ec         er        gc        gr        hp
## [1,] 0.07045486 0.09247026 0.6306444 0.1419681 0.0644623

割引前と比較

model_con %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         hp         ec         er         gc         gr 
## 0.05555556 0.07111111 0.09333333 0.63666667 0.14333333

割引によって, ヒートポンプの割合が5.5%から6.4%に上昇することがこのモデルから予測される.

【学習予定】

次は, ネスティッドロジットモデルの節を扱う.

moratoriamuo.hatenablog.com

ロジットモデルの人工データ生成と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の章に書かれている)。

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

ライブラリ読み込み。

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}\)の符号から違うという結果に。

失敗例

推定が上手くいった例を上に載せたが, 失敗例も役立つかもしれないので載せておく. 失敗例を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`.