【学習動機】
少し複雑なロジットモデルを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.05555556gc : 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.11315710900人の平均選択確率
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.841459model_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.24alt <- 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%に上昇することがこのモデルから予測される.
