パッケージ {mlogit} の利用法
パッケージ {mlogit} を用いた離散選択モデルの推定法について解説します．パッケージのマニュアル（英語）はこちらで閲覧できます．パッケージのチュートリアル（英語）はこちら（中段の Vignetts）で閲覧できます．
１．多項ロジットモデル
以下のコードを Rスクリプトファイルに記述していきます．
（１）ライブラリの読み込み
library(mlogit)（２）データの読み込み
今回は data( ) 関数を用いてパッケージで用意されているサンプルデータを利用します．
# {mlogit} パッケージの "Heating" というサンプルデータの読み込み
data("Heating", package = "mlogit")Environment タブに Heating という Data を確認できます．クリックするとデータの中身を確認できます．表形式のデータになっていますが R ではデータフレーム型と呼ばれます．整数型・実数型・文字列型・因子型などの異なるデータ型の列ベクトルを１つにまとめたものと考えてください．
・他のサンプルデータを確認したい場合には data( ) をコンソールで実行してください．
今回は必要ありませんが Excel で作成したデータを読み込みたい場合は Rスクリプトに下記の命令文を書いてください．
# ワーキングディレクトリの変更
# Rスクリプトとデータを同じフォルダに置いておきます
setwd("C:/Users/fukumoto/ドキュメント/R/講義/都市システム計画演習/R04/")
# ワーキングディレクトリの取得
wd1 <- getwd()
# ファイル名の指定
file1 <- "/データ.xlsx"
# ディレクトリ＋ファイル名
path1 <- paste0(wd1, file1)
# データの読み込み
# {readxl}パッケージの関数 read_excel を使用
# シート1にデータを入れておく．1行目は各列の名前，2行目以降に各サンプルの情報．
df1 <- readxl::read_excel(path1, sheet = 1)先ほど読み込んだデータの概要をコンソールで確認するには以下のようにします．
# データフレームの構造
str(Heating)
# 先頭10行の表示
head(Heating, 10)
# 列ベクトルの表示
# データフレームの特定の列にアクセスする場合に "$" を使います．
Heating$depvar # 列ベクトル
Heating$depvar[1:10] # 列ベクトルの要素へのアクセス（３）データ形式
{mlogit} パッケージでは Wide形式 と Long形式 の異なるデータ形式を扱えます．ここでは Wide形式について説明します．
上のデータは，調査対象の各個人の最大16回の選択行動（機関Aか機関Bの二択）の結果を記録したものです．"id" が個人のID，"choiceid" が各個人の選択機会のIDを表します．"choice" は実際に選択した行動，"変数_A/B" は代替案A/Bの属性情報になります．属性情報は以下の４つです．
価格
時間
乗り換え回数
快適さ（０・１・２の３段階；０が最も快適なクラス）
このデータをmlogitで分析するには以下のようにデータ形式をmlogit.data型に変換します．
# 2929個のサンプルを一意に区別するために choiceid に通し番号を振ります．
Train$choiceid <- 1:nrow(Train)
# 名前の混同を避けるため，choice を choice.alt に名前を変更します．
names(Train)[3] <- "choice.alt"
# mlogit.data型への変更
# choice は選択結果，varying は各代替案の共変量に対応します．
# sep は変数名の区切り文字です．
Tr <- mlogit.data(Train, shape = "wide",
choice = "choice.alt", varying = c(4:11), sep = "_")
# 単位の変更
exchange.rate = 2.20371 #為替レート
Tr$price <- Tr$price / 100 * exchange.rate
Tr$time <- Tr$time / 60
# 中身の確認
head(Tr,10)変換後のデータの中身は以下のとおりです．変換前は選択機会別にサンプルが並んでいましたが変換後は選択機会・選択肢別にサンプルが並んでいます．それに伴って新たに chid列とidx 列が追加されており下図下段にあるように ”選択機会, 選択肢” という情報を格納しています．
（４）ロジットモデルの推定
個人 i が代替案 j を選択した場合の効用が次式で表されるとします．
$$V_{ij} = \alpha_j+\beta x_{ij} +\gamma_jz_i+\delta_jw_{ij}+\epsilon_{ij} $$
αj は代替案 j の定数項です．ロジットモデルでは効用の差によって選択確率が決まるため，N個の代替案がある多項ロジットモデルでは，N-1個の代替案の定数項が推定されます．
xij とwij は個人i と代替案j により異なる変数（機関Aを利用する場合の運賃や時間），zi は個人i により異なる変数（年収や地域ダミー）
xij はすべての代替案の効用に同じように影響するのに対し，wij は代替案ごとに影響の大きさが異なります．例えば，運賃や時間は前者に相当し，乗り換え回数は交通機関（あるいは結節点）によって影響が異なり後者に相当する可能性があります．
モデルの推定と結果の表示は以下の命令文で実行できます．mlogit関数を用いて，データ "Tr" を利用して，"choice.alt" を選択結果，"price" と "time" を上式の xij とする多項ロジットモデルを推定しています．代替案の定数項は自動的に加わります．
分析結果を res0 という変数に格納しており，結果の表示には summary関数を用います．
# モデルの推定
res0 <- mlogit(choice.alt ~ price + time, Tr)
# 結果の表示
summary(res0)
## Call:
## mlogit(formula = choice.alt ~ price + time, data = Tr, method = "nr")
##
## Frequencies of alternatives:choice
## A B
## 0.50324 0.49676
##
## nr method
## 4 iterations, 0h:0m:0s
## g'(-H)^-1g = 1.87E-07
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):B -1.8735e-02 3.9362e-02 -0.4760 0.6341
## price -1.0237e-03 5.9389e-05 -17.2364 < 2.2e-16 ***
## time -1.3968e-02 2.2876e-03 -6.1061 1.021e-09 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Log-Likelihood: -1845.5
## McFadden R^2: 0.090937
## Likelihood ratio test : chisq = 369.23 (p.value = < 2.22e-16)Coefficients 欄の Estimate がパラメータの推定値になります．Std.Error（標準誤差）は推定結果のばらつきの大きさを表し，Pr(>|z|) （p値）はパラメータの真値が 0 とした場合に今回の推定値が得られる確率を表します．p値が小さいほど真値が0でない可能性が大きいことを意味し，逆にp値が大きい場合には統計学的に有意な変数とはいえないのでモデルから落とすことを検討する必要があります．
上の推定結果では定数項が明らかに有意でないので外すことを検討すべきです．その場合は以下のとおり推定します．
# "choice.alt ~ xij | zi | wij" に対応し，0 は該当変数がないことを意味します．
res1 <- mlogit(choice.alt ~ price + time | 0 | 0, Tr)
summary(res1)以下のとおり変数の組み合わせを変更することもできます．
res2 <- mlogit(choice.alt ~ price + time | 0 | change + comfort, Tr)
summary(res2)この記事が気に入ったらサポートをしてみませんか？