多項ロジット(Multinomial Logit), R - mlogit 使用メモ
Rのmlogitパッケージで多項ロジット(Multinomial Logit)を使用する際のメモ。
まず、用語の整理。参考文献(A) p.8より。
-------------------------
A model with only individual specific variables is sometimes called a multinomial logit (多項ロジット)model, one with only alternative specific variables a conditional logit(条件付ロジット) model and one with both kind of variables a mixed logit(混合ロジット) model.
This is seriously misleading: conditional logit model is also a logit model for longitudinal data in the statistical literature and mixed logit is one of the names of a logit model with random parameters. Therefore, in what follow, we'll use the name multinomial logit model for the model we've just described whatever the nature of the explanatory variables included in the model.
-------------------------
とのことなので、多項ロジット(Multinomial Logit)と呼ぶことにする。
指定する説明変数がややこしいので、整理しておく。
以下のように効用関数が定義されるとする。
i : 個人の添え字
j : 選択肢の添え字
αj は、選択肢jの基本的効用と考えられる。
β : 個人iと選択肢jによって異なる変数Xijに対して、共通の係数
γj : 個人iに固有の変数Ziに対して、選択肢j別に異なる係数
σj : 個人iと選択肢jによって異なる変数Wijに対して、選択肢j別に異なる係数
library(mlogit)
data(Fishing)
head(Fishing, 10)
mode price.beach price.pier price.boat price.charter catch.beach 1 charter 157.930 157.930 157.930 182.930 0.0678 2 charter 15.114 15.114 10.534 34.534 0.1049 3 boat 161.874 161.874 24.334 59.334 0.5333 4 pier 15.134 15.134 55.930 84.930 0.0678 5 boat 106.930 106.930 41.514 71.014 0.0678 6 charter 192.474 192.474 28.934 63.934 0.5333 7 beach 51.934 51.934 191.930 220.930 0.0678 8 charter 15.134 15.134 21.714 56.714 0.0678 9 boat 34.914 34.914 34.914 53.414 0.2537 10 boat 28.314 28.314 28.314 46.814 0.2537 catch.pier catch.boat catch.charter income 1 0.0503 0.2601 0.5391 7083.332 2 0.0451 0.1574 0.4671 1250.000 3 0.4522 0.2413 1.0266 3750.000 4 0.0789 0.1643 0.5391 2083.333 5 0.0503 0.1082 0.3240 4583.332 6 0.4522 0.1665 0.3975 4583.332 7 0.0789 0.1643 0.5391 8750.001 8 0.0789 0.0102 0.0209 2083.333 9 0.1498 0.0233 0.0219 3750.000 10 0.1498 0.0233 0.0219 2916.667
一列が1個人となっている”wide型”と呼ばれるデータ。
str(Fishing)
'data.frame': 1182 obs. of 10 variables: $ mode : Factor w/ 4 levels "beach","pier",..: 4 4 3 2 3 4 1 4 3 3 ... $ price.beach : num 157.9 15.1 161.9 15.1 106.9 ... $ price.pier : num 157.9 15.1 161.9 15.1 106.9 ... $ price.boat : num 157.9 10.5 24.3 55.9 41.5 ... $ price.charter: num 182.9 34.5 59.3 84.9 71 ... $ catch.beach : num 0.0678 0.1049 0.5333 0.0678 0.0678 ... $ catch.pier : num 0.0503 0.0451 0.4522 0.0789 0.0503 ... $ catch.boat : num 0.26 0.157 0.241 0.164 0.108 ... $ catch.charter: num 0.539 0.467 1.027 0.539 0.324 ... $ income : num 7083 1250 3750 2083 4583 ...
Mlogitによってモデルをあてはめるには、mlogit.data型に変換する必要がある。
Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice ="mode")
Fishingはサンプルデータで、実データで実施すると上手く行かないことが多いが、とりあえず続ける。
ここら辺の変換はhelpを参考にすればよい。
head(Fish, 12)
mode income alt price catch chid 1.beach FALSE 7083.332 beach 157.930 0.0678 1 1.boat FALSE 7083.332 boat 157.930 0.2601 1 1.charter TRUE 7083.332 charter 182.930 0.5391 1 1.pier FALSE 7083.332 pier 157.930 0.0503 1 2.beach FALSE 1250.000 beach 15.114 0.1049 2 2.boat FALSE 1250.000 boat 10.534 0.1574 2 2.charter TRUE 1250.000 charter 34.534 0.4671 2 2.pier FALSE 1250.000 pier 15.114 0.0451 2 3.beach FALSE 3750.000 beach 161.874 0.5333 3 3.boat TRUE 3750.000 boat 24.334 0.2413 3 3.charter FALSE 3750.000 charter 59.334 1.0266 3 3.pier FALSE 3750.000 pier 161.874 0.4522 3
4行(選択肢の総数)で1個人のデータとなる。
str(Fish)
Classes ‘mlogit.data’ and 'data.frame': 4728 obs. of 6 variables: $ mode : logi FALSE FALSE TRUE FALSE FALSE FALSE ... $ income: num 7083 7083 7083 7083 1250 ... $ alt : chr "beach" "boat" "charter" "pier" ... $ price : num 157.9 157.9 182.9 157.9 15.1 ... $ catch : num 0.0678 0.2601 0.5391 0.0503 0.1049 ... $ chid : int 1 1 1 1 2 2 2 2 3 3 ... - attr(*, "reshapeLong")=List of 4 ..$ varying:List of 2 .. ..$ price: chr "price.beach" "price.pier" "price.boat" "price.charter" .. ..$ catch: chr "catch.beach" "catch.pier" "catch.boat" "catch.charter" .. ..- attr(*, "v.names")= chr "price" "catch" .. ..- attr(*, "times")= chr "beach" "pier" "boat" "charter" ..$ v.names: chr "price" "catch" ..$ idvar : chr "chid" ..$ timevar: chr "alt" - attr(*, "index")='data.frame': 4728 obs. of 2 variables: ..$ chid: Factor w/ 1182 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ... ..$ alt : Factor w/ 4 levels "beach","boat",..: 1 2 3 4 1 2 3 4 1 2 ... - attr(*, "choice")= chr "mode"
目的変数となるmodeはlogical型になる。
もう一つ、サンプルデータ。
library(AER)
data(TravelMode)
str(TravelMode)
'data.frame': 840 obs. of 9 variables: $ individual: Factor w/ 210 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ... $ mode : Factor w/ 4 levels "air","train",..: 1 2 3 4 1 2 3 4 1 2 ... $ choice : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 2 1 1 ... $ wait : int 69 34 35 0 64 44 53 0 69 34 ... $ vcost : int 59 31 25 10 58 31 25 11 115 98 ... $ travel : int 100 372 417 180 68 354 399 255 125 892 ... $ gcost : int 70 71 70 30 68 84 85 50 129 195 ... $ income : int 35 35 35 35 30 30 30 30 40 40 ... $ size : int 1 1 1 1 2 2 2 2 1 1 ...
head(TravelMode, 10)
individual mode choice wait vcost travel gcost income size 1 1 air no 69 59 100 70 35 1 2 1 train no 34 31 372 71 35 1 3 1 bus no 35 25 417 70 35 1 4 1 car yes 0 10 180 30 35 1 5 2 air no 64 58 68 68 30 2 6 2 train no 44 31 354 84 30 2 7 2 bus no 53 25 399 85 30 2 8 2 car yes 0 11 255 50 30 2 9 3 air no 69 115 125 129 40 1 10 3 train no 34 98 892 195 40 1
“long型”と呼ばれるデータ。
mlogit.data型に変換。
TM <- mlogit.data(TravelMode, choice = "choice", shape = "long", id.var = "individual", alt.levels = c("air", "train", "bus", "car"))
str(TM)
Classes ‘mlogit.data’ and 'data.frame': 840 obs. of 9 variables: $ individual: Factor w/ 210 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ... $ mode : Factor w/ 4 levels "air","train",..: 1 2 3 4 1 2 3 4 1 2 ... $ choice : logi FALSE FALSE FALSE TRUE FALSE FALSE ... $ wait : int 69 34 35 0 64 44 53 0 69 34 ... $ vcost : int 59 31 25 10 58 31 25 11 115 98 ... $ travel : int 100 372 417 180 68 354 399 255 125 892 ... $ gcost : int 70 71 70 30 68 84 85 50 129 195 ... $ income : int 35 35 35 35 30 30 30 30 40 40 ... $ size : int 1 1 1 1 2 2 2 2 1 1 ... - attr(*, "index")='data.frame': 840 obs. of 3 variables: ..$ chid: Factor w/ 210 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ... ..$ alt : Factor w/ 4 levels "air","train",..: 1 2 3 4 1 2 3 4 1 2 ... ..$ id : Factor w/ 210 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ... - attr(*, "choice")= chr "choice"
head(TM, 12)
individual mode choice wait vcost travel gcost income size 1.air 1 air FALSE 69 59 100 70 35 1 1.train 1 train FALSE 34 31 372 71 35 1 1.bus 1 bus FALSE 35 25 417 70 35 1 1.car 1 car TRUE 0 10 180 30 35 1 2.air 2 air FALSE 64 58 68 68 30 2 2.train 2 train FALSE 44 31 354 84 30 2 2.bus 2 bus FALSE 53 25 399 85 30 2 2.car 2 car TRUE 0 11 255 50 30 2 3.air 3 air FALSE 69 115 125 129 40 1 3.train 3 train FALSE 34 98 892 195 40 1 3.bus 3 bus FALSE 35 53 882 149 40 1 3.car 3 car TRUE 0 23 720 101 40 1
4行(選択肢の総数)で1個人。
Helpにデータ項目の簡単な説明がある。
210人のデータである。
選択に関するクロス集計。
( actChoice <- table(TM$mode, TM$choice) )
FALSE TRUE air 152 58 train 147 63 bus 180 30 car 151 59
割合で出す場合。
( actChoiceRatio <- actChoice[,2]/210 ) # 210 = #individuals
air train bus car 0.2761905 0.3000000 0.1428571 0.2809524
では、多項ロジットモデルのあてはめを行う。
FitTM <- mlogit(choice ~ vcost | income | travel, data=TM)
formulaは、“|”が仕切りとなって分かれている。
これは、上で説明したように、Xij | Zi | Wij とそれぞれの変数の役割に従って指定する。
vcostは各個人iが提示された各選択肢jの値段(各乗り物の値段)である。個人iによって同じ選択肢jでも値段が異なっていることに注意。もしi間で共通の変数であれば、選択肢と全く同じとなってしまうので、モデルに含めることはできない。
vcost変数の、選択肢jによって異ならない共通の効用への影響度を推定する(β)。
incomeは収入で、選択肢jに依存しない、個人i依存の変数である。
incom変数の、効用への影響度を選択肢j別に推定する(γj)。
travelは各個人iが選択肢jを選択した場合に経験するトラベルタイムである。
travel変数の、効用への影響度を選択肢j別に推定する(σj)。トラベルタイムに対して感じる効用が、選択肢(乗り物)によって異なると仮定する場合はここに変数を入れるのだと思われる。
結果の確認。
summary( FitTM )
Call: mlogit(formula = choice ~ vcost | income | travel, data = TM, method = "nr", print.level = 0) Frequencies of alternatives: air train bus car 0.27619 0.30000 0.14286 0.28095 nr method 5 iterations, 0h:0m:0s g'(-H)^-1g = 7.87E-05 successive function values within tolerance limits Coefficients : Estimate Std. Error t-value Pr(>|t|) train:(intercept) 1.6866632 0.8477430 1.9896 0.046636 * bus:(intercept) 0.1449168 0.9201414 0.1575 0.874856 car:(intercept) -0.7699463 0.8161644 -0.9434 0.345491 vcost -0.0056710 0.0066967 -0.8468 0.397089 train:income -0.0598253 0.0129413 -4.6228 3.786e-06 *** bus:income -0.0414433 0.0136923 -3.0268 0.002472 ** car:income -0.0129009 0.0111799 -1.1539 0.248527 air:travel -0.0367774 0.0067680 -5.4340 5.509e-08 *** train:travel -0.0074593 0.0011859 -6.2898 3.180e-10 *** bus:travel -0.0068219 0.0013397 -5.0922 3.539e-07 *** car:travel -0.0066324 0.0011331 -5.8531 4.825e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log-Likelihood: -232.48 McFadden R^2: 0.18071 Likelihood ratio test : chisq = 102.56 (p.value = < 2.22e-16)
ここで出てくる“Frequencies of alternatives:”は上で計算したとおり、実測選択の割合である。予測選択の割合ではない。
切片(intercept)とγ、σは一つの水準が0に固定される。
0にする水準を入れ替える場合は、引数reflevelを指定する。
summary( mlogit(choice ~ vcost | income | travel, reflevel=4, data=TM) )
Call: mlogit(formula = choice ~ vcost | income | travel, data = TM, reflevel = 4, method = "nr", print.level = 0) Frequencies of alternatives: car air train bus 0.28095 0.27619 0.30000 0.14286 nr method 5 iterations, 0h:0m:0s g'(-H)^-1g = 7.87E-05 successive function values within tolerance limits Coefficients : Estimate Std. Error t-value Pr(>|t|) air:(intercept) 0.7699463 0.8161644 0.9434 0.3454908 train:(intercept) 2.4566096 0.6244004 3.9343 8.342e-05 *** bus:(intercept) 0.9148631 0.7662845 1.1939 0.2325191 vcost -0.0056710 0.0066967 -0.8468 0.3970894 air:income 0.0129009 0.0111799 1.1539 0.2485269 train:income -0.0469243 0.0121817 -3.8520 0.0001171 *** bus:income -0.0285424 0.0129991 -2.1957 0.0281118 * car:travel -0.0066324 0.0011331 -5.8531 4.825e-09 *** air:travel -0.0367774 0.0067680 -5.4340 5.509e-08 *** train:travel -0.0074593 0.0011859 -6.2898 3.180e-10 *** bus:travel -0.0068219 0.0013397 -5.0922 3.539e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log-Likelihood: -232.48 McFadden R^2: 0.18071 Likelihood ratio test : chisq = 102.56 (p.value = < 2.22e-16)
適合度の指標に関して。
“Log-Likelihood:”は対数尤度。
“McFadden R^2:”は参考文献(B) p.124では、McFadden ρ^2 と書かれている。
“Likelihood ratio test :”は切片モデルに対する検定。
“McFadden R^2:”を計算してみる。
( ll <- summary( FitTM )$logLik )
'log Lik.' -232.4804 (df=11)
切片モデル。
FitTM0 <- mlogit(choice ~ 1, data=TM)
summary( FitTM0 )
Call: mlogit(formula = choice ~ 1, data = TM, method = "nr", print.level = 0) Frequencies of alternatives: air train bus car 0.27619 0.30000 0.14286 0.28095 nr method 4 iterations, 0h:0m:0s g'(-H)^-1g = 0.000622 successive function values within tolerance limits Coefficients : Estimate Std. Error t-value Pr(>|t|) train:(intercept) 0.082692 0.181974 0.4544 0.649529 bus:(intercept) -0.659237 0.224888 -2.9314 0.003374 ** car:(intercept) 0.017094 0.184907 0.0924 0.926341 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log-Likelihood: -283.76 McFadden R^2: 0 Likelihood ratio test : chisq = 0 (p.value = 1)
( ll0 <- summary( FitTM0 )$logLik )
'log Lik.' -283.7588 (df=3)
as.numeric( 1 - (ll / ll0) )
[1] 0.1807111
となる。
対数尤度を利用した計算からも分かるとおり、モデルに変数を加えれば加えるほど値が高くなっていくことに注意。
この指標の読み方であるが、参考文献(B)が(C)を参照している。
(参考文献(C) p.124より)
横軸が“McFadden R^2:”で、縦軸が通常のR^2である。
どうやって導き出したかは不明だが、経験則的にこのように判断できるとなっている。0.5程度のモデルが作れれば、十分合格点ということになるだろう。
各個人に対する、各選択肢の選択予測値は以下の操作で出力を行う。
( fitRes <- FitTM$probabilities )
# fitted(FitTM, outcome=FALSE) # これでも同じ
air train bus car [1,] 0.119817364 0.23053702 0.090576734 0.559068885 [2,] 0.317031379 0.28838285 0.102175448 0.292410320 [3,] 0.651515719 0.04527835 0.049321451 0.253884481 [4,] 0.492368062 0.03999761 0.029395135 0.438239197 [5,] 0.117289099 0.49291692 0.237541845 0.152252133 [6,] 0.087258878 0.53618734 0.131474081 0.245079699 [7,] 0.829788104 0.02232870 0.027580156 0.120303036 [8,] 0.251498768 0.34596333 0.176161979 0.226375922 [9,] 0.154158841 0.21992794 0.094726675 0.531186546 [10,] 0.526192120 0.04275839 0.031663300 0.399386190 [11,] 0.179378500 0.42617592 0.112891226 0.281554353 [12,] 0.114347018 0.29971102 0.117604796 0.468337167 [13,] 0.144166110 0.17389664 0.089848946 0.592088309 [14,] 0.223915354 0.27859546 0.119081403 0.378407786 [15,] 0.498162161 0.25977304 0.200121678 0.041943118 ..... skip ..... [207,] 0.218547728 0.24899471 0.211234675 0.321222889 [208,] 0.056502916 0.56022179 0.276306236 0.106969061 [209,] 0.170163004 0.37815127 0.183175726 0.268510000 [210,] 0.405440673 0.04138912 0.057691327 0.495478876
各個人の実際の選択に対する予測値のみを返したい場合。
FitTM$fitted.values
# fitted( FitTM, outcome=TRUE ) # これでも同じ
1.air 2.air 3.air 4.air 5.air 6.air 0.559068885 0.292410320 0.253884481 0.438239197 0.152252133 0.536187342 7.air 8.air 9.air 10.air 11.air 12.air 0.829788104 0.226375922 0.531186546 0.399386190 0.281554353 0.468337167 13.air 14.air 15.air 16.air 17.air 18.air 0.592088309 0.378407786 0.041943118 0.625623550 0.558350920 0.738326313 ..... skip ..... 199.air 200.air 201.air 202.air 203.air 204.air 0.510837490 0.361912483 0.258362415 0.315766795 0.202403945 0.267393122 205.air 206.air 207.air 208.air 209.air 210.air 0.166773948 0.128082868 0.218547728 0.276306236 0.268510000 0.495478876
混同行列(Confusion Matrix)を作成してみる。
標準出力にあるのかもしれないが、分からないので、選択予測確率(fitResに格納している)から作成してみる。
このfitResの各行において、最大の確率に対応する列名を取り出す関数を作成。
## fitResの各行において、maxとなる列名をまとめたvectorを返す関数
fn_findPred1 <- function(df0) {
df1 <- as.data.frame(df0)
maxLevel <- apply(df1, 1, max)
nc <- ncol(df1)
temp_return <- NULL
for(i in 1:nc) {
#print(i)
temp_col <- ifelse( df1[i]==maxLevel, names(df1[i]), "" )
#print(head(temp_col))
if(i==1){
temp_return <- temp_col
}else{
temp_return <- paste(temp_return, temp_col, sep="")
}
#print(head(temp_return))
}
return( temp_return )
}
fitResを代入して関数を実行。
( predChoiceCol <- fn_findPred1(fitRes) )
[1] "car" "air" "air" "air" "train" "train" "air" "train" "car" [10] "air" "train" "car" "car" "car" "air" "train" "train" "train" [19] "train" "air" "train" "train" "air" "air" "train" "air" "air" [28] "train" "train" "train" "train" "train" "air" "train" "train" "train" ..... skip ..... [181] "train" "train" "train" "train" "train" "air" "train" "train" "car" [190] "air" "train" "train" "car" "bus" "train" "train" "train" "air" [199] "air" "air" "train" "air" "train" "train" "air" "air" "car" [208] "train" "train" "car"
各個人の予測選択肢が出力される。
この予測選択による選択肢の割合を計算するには以下。
table(predChoiceCol)/210
predChoiceCol air bus car train 0.30000000 0.04761905 0.25238095 0.40000000
各個人の実際の選択をデータ(TM)から作成。
( actChoiceCol <- as.character(TM$mode[TM$choice==TRUE]) )
[1] "car" "car" "car" "car" "car" "train" "air" "car" "car" [10] "car" "car" "car" "car" "car" "car" "train" "train" "train" [19] "train" "train" "train" "train" "air" "air" "air" "air" "air" [28] "air" "train" "train" "train" "train" "train" "train" "train" "train" ..... skip ..... [181] "car" "car" "train" "train" "train" "bus" "air" "air" "car" [190] "car" "car" "air" "air" "bus" "train" "train" "car" "air" [199] "air" "air" "bus" "car" "bus" "bus" "bus" "car" "air" [208] "bus" "car" "car"
predChoiceColとactChoiceColでクロス集計。
( ConfMatrix <- table(actChoiceCol, predChoiceCol) )
predChoiceCol actChoiceCol air bus car train air 25 0 20 13 bus 9 10 3 8 car 17 0 30 12 train 12 0 0 51
と混同行列(Confusion Matrix)が作成できる。
誤分類率を求める場合は、混同行列の対角要素を利用し、以下のように計算できる。
1 - sum(diag(ConfMatrix))/210 # 210 = #individuals
[1] 0.447619
では、変数を増やしたモデルも作成しておく。
FitTM2 <- mlogit(choice ~ wait+vcost | income+size | travel, data=TM)
summary( FitTM2 )
Call: mlogit(formula = choice ~ wait + vcost | income + size | travel, data = TM, method = "nr", print.level = 0) Frequencies of alternatives: air train bus car 0.27619 0.30000 0.14286 0.28095 nr method 5 iterations, 0h:0m:0s g'(-H)^-1g = 7.32E-07 gradient close to zero Coefficients : Estimate Std. Error t-value Pr(>|t|) train:(intercept) -2.1240906 1.2159544 -1.7469 0.0806633 . bus:(intercept) -3.3311281 1.3317591 -2.5013 0.0123739 * car:(intercept) -7.5393851 1.2585583 -5.9905 2.092e-09 *** wait -0.0973882 0.0113005 -8.6180 < 2.2e-16 *** vcost -0.0058938 0.0086924 -0.6780 0.4977476 train:income -0.0727496 0.0172146 -4.2260 2.378e-05 *** bus:income -0.0333524 0.0176288 -1.8919 0.0585002 . car:income -0.0163247 0.0141088 -1.1571 0.2472481 train:size 1.1771968 0.3094744 3.8039 0.0001425 *** bus:size 0.8055721 0.3962361 2.0331 0.0420464 * car:size 1.0019958 0.2754798 3.6373 0.0002755 *** air:travel -0.0329873 0.0072772 -4.5330 5.816e-06 *** train:travel -0.0069413 0.0013858 -5.0090 5.472e-07 *** bus:travel -0.0067683 0.0016560 -4.0872 4.367e-05 *** car:travel -0.0070640 0.0012872 -5.4879 4.068e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log-Likelihood: -162.79 McFadden R^2: 0.42631 Likelihood ratio test : chisq = 241.94 (p.value = < 2.22e-16)
McFrdden R^2など大きく上昇している。
各個人における、選択肢の予測確率。
fitRes2 <- fitted(FitTM2, outcome=FALSE)
上で定義した関数を適用。
predChoiceCol2 <- fn_findPred1(fitRes2)
混同行列。
( ConfMatrix2 <- table(actChoiceCol, predChoiceCol2) )
predChoiceCol2 actChoiceCol air bus car train air 40 1 14 3 bus 5 22 2 1 car 5 0 46 8 train 5 1 9 48
誤分類率。
1 - sum(diag(ConfMatrix2))/210
[1] 0.2571429
実データで扱う場合は、それはそれで大変なのだが、今回は使い方メモということで…。
---------- 参考文献 ----------
(A) Estimation of multinomial logit models in R : The mlogit Packages
Yves Croissant
http://cran.r-project.org/web/packages/mlogit/vignettes/mlogit.pdf
(B) マーケティングのデータ分析―分析手法と適用事例
岡太彬訓, 守口剛
http://www.asakura.co.jp/books/isbn/978-4-254-12822-2/
(C) Urban Travel Demand: A Behavioral Analysis
Tom Domencich, Daniel L. McFadden
http://eml.berkeley.edu/~mcfadden/travel.html
Ch.5
http://eml.berkeley.edu/~mcfadden/travel/ch5.pdf