【学習動機】
前回の続きで少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.
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
また, Trainの書いたテキストは王道の本だと思うので, リンクを貼っておく.
【学習内容】
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} })\)
(なんか表示が上手くいかなかったのでスクショ追加)
同じネストに含まれる選択肢\(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}}}}\)
(なんか表示が上手くいかなかったのでスクショ追加(再度))
データの確認
カリフォルニア州の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.032wide形式から変換
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 1mlogit()でネスティッドロジットモデルするときの引数：
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.12log-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.04ivの値がそれぞれ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.6752つのネストの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)の節を扱うかも. (やった)
