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

ネスティッドロジットモデル（Nested Logit）をoptimを使って解きました。
階層構造（入れ子）になっている選択を通常のロジットで解くとIIA問題があり、上手く推定できません。

対処方法としては、プロビットモデルで解くか、ネスティッドロジットモデル（Nested Logit）を使って解くかのアプローチになるわけですが、ネスティッドロジットモデルに対応しているパッケージなどが少なく、自分で尤度関数を作成して解くのが一番分かりやすそうです。
Rだとmlogitというパッケージがありますね。

【ポイント 1】
log sum で上がってくる部分のλは0≦λ≦1である必要があるため、
exp(bx) / (1 + exp(bx))
という形式で与えることで、0≦λ≦1を保証しています。

【ポイント 2】
optim 関数はデフォルトで最小化を行うため、control 引数に fnscale = -1 を代入して、最大化を行うようにする

【Nested LogitのRコード】
hh <- nrow(dat_ss) # レコード数
b0 <- c(0, 0, 0, 0, 0) # パラメータの初期化

# ロジットモデルの尤度関数
fr <- function(x)
{
b1 <- x[1] # 階層1の説明変数 1
b2 <- x[2] # 階層1の説明変数 2
b3 <- x[3] # 階層2の説明変数 1
b4 <- x[4] # 階層2の説明変数 2
bx <- x[7] # 0 <= bx <= 1

LL = 0
for(i in 1:hh) {
# 効用の計算
　 U_K2 <- b3 * dat_ss[i,3] + b4*dat_ss[i,4]
U_K1 <- b1 * dat_ss[i,1] + b2*dat_ss[i,2] + (exp(bx)/(1+exp(bx)))*log(1+exp(U_K2))

# 選択確率の計算
PP1 <- exp(U_K1)/(exp(U_K1)+1)
PP2 <- exp(U_K2)/(exp(U_K2)+1)

# 対数尤度の計算
LLL <- dat_ss[i,5]*log(PP1)+(1-dat_ss[i,5])*log(1-PP1)+dat_ss[i,5]*(dat_ss[i,6]*log(PP2)+(1-dat_ss[i,6])*log(1-PP2))
LL <- LL + LLL
}
return(LL)
}

# 対数尤度関数の最大化
res_nl <- optim(b0, fr, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))