多項ロジットモデル Last modified: Jul 27, 2011
目的
牧厚志ら「経済・経営のための統計学」有斐閣アルマの第 8 章「マーケティング ブランド選択行動の分析」濱岡豊
本文に示されている R プログラムを汎用化して使いやすくした。
使用法
mlm(d, x=NULL, z=NULL, constants=TRUE)
引数
d 本文中で $d_{ij}$ と記されているもの
x 本文中で $x_{ijm}$ と記されているもの（変数の順序（列）に注意）
z 本文中で $z_{i}$ と記されているもの
constants 本文中で「ブランド定数」と記されているものを加えるか否か
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/mlm.R", encoding="euc-jp")
# 牧厚志ら「経済・経営のための統計学」有斐閣アルマの
# 第8章「マーケティング ブランド選択行動の分析」濱岡豊
# 本文に示されている R プログラムを一般化して使いやすくした
mlm <- function( d, # 本文中で $d_{ij}$ と記されているもの
x=NULL, # 本文中で $x_{ijm}$ と記されているもの（変数の順序（列）に注意）
z=NULL, # 本文中で $z_{i}$ と記されているもの
constants=TRUE) # 本文中で「ブランド定数」と記されているものを加えるか否か
{ # x, z, constant は 1 つ以上選択可。途中の引数が省略される場合には引数名付きで指定すること
getV <- function(par) # 効用 $V_{ij}$ の計算
{
V <- matrix(0, n, N)
if (include.x) { # $x_{ijm}$ について
beta <- par[1:p]
V <- V+t(apply(x3, 3, function(z) z%*%beta))
}
if (include.z) { # $z_{i}$ について
beta <- cbind(matrix(par[include.x*p+1:(ncol.z*(n-1))], ncol.z, n-1), rep(0, ncol.z))
V <- V+t(z%*%beta)
}
if (constants) { # ブランド定数について
V <- V+c(par[(n.parameters-n+2):n.parameters], 0)
}
return(V)
}
mlogit <- function(par) # 対数尤度を求める関数
{
V <- getV(par)
return(sum(V[cbind(d, 1:N)]-log(colSums(exp(V)))))
}
include.x <- !is.null(x) # x が指定されると TRUE
include.z <- !is.null(z) # z が指定されると TRUE
stopifnot(include.x || include.z || constants) # どれか 1 つは指定されていないといけない
if (is.data.frame(d)) { # データフレームなら行列にしておく
d <- data.matrix(d)
}
n <- length(unique(d)) # 選択対象の種類
N <- length(d) # サンプルサイズ
if (include.x && is.data.frame(x)) { # データフレームなら行列にしておく
x <- data.matrix(x)
p <- ncol(x)/n # 選択対象の属性の種類
if (p != floor( p )) {
stop("属性の総項目数が選択対象の種類の整数倍になっていない")
}
x3 <- array(unlist(x), dim=c(N, p, n)) # 選択対象別に計算するために 3 次元配列にする
}
if (include.z && is.data.frame(z)) { # データフレームなら行列にしておく
z <- data.matrix(z)
}
n.parameters <- 0 # 推定すべきパラメータ数
if (include.x) {
n.parameters <- n.parameters+p # + 選択対象の属性の種類
}
if (include.z) {
ncol.z <- ncol(z) # + z として使う変数の個数*(選択対象の種類-1)
n.parameters <- n.parameters+ncol.z*(n-1)
}
if (constants) {
n.parameters <- n.parameters+n-1 # + 選択対象の種類-1
}
par <- numeric(n.parameters) # パラメータベクトル（初期値は 0 でよい）
res <- optim(par, mlogit, method="BFGS", hessian=TRUE, control=list(fnscale=-1))
LL <- res$value # 対数尤度
beta.estimated <- res$par # 推定されたパラメータ
AIC <- -2*(LL-n.parameters) # AIC
ses <- sqrt(diag(solve(-res$hessian))) # 標準誤差
t <- abs(beta.estimated)/ses # t 値
P <- pt(t, N-n.parameters, lower.tail=FALSE)*2 # P 値
df <- data.frame(beta.estimated, ses, t=round(t, 3), P=round(P, 3))
if(include.x) { # 名前の編集
rownames(df)[1:p] <- c(sub("[0-9]*$", "", colnames(x)[1:p]))
}
if (include.z) {
for (i in 1:ncol.z) {
rownames(df)[include.x*p++(i-1)*(n-1)+1:(n-1)] <- paste(colnames(z)[i], 1:(n-1), sep="-")
}
}
if (constants) {
rownames(df)[(n.parameters-n+2):n.parameters] <- paste("Constant", 1:(n-1), sep="-")
}
V <- t(getV(beta.estimated)) # 効用
v <- exp(V)
P <- v/rowSums(v) # 選択確率
predict <- apply(P, 1, which.max) # 選択予測（P が一番大きいもの）
res <- list(df=df, LL=LL, AIC=AIC, V=V, P=P, predict=predict)
class(res) <- "mlm"
return(res)
}
summary.mlm <- function(object, digits=5) # sumamry メソッド
{
print(object$df, digits=digits)
cat("LL =", object$LL, "\nAIC =", object$AIC, "\n")
}
print.mlm <- function(object, digits=5) # print メソッド
{
print.default(object, digits=digits)
}
predict.mlm <- function(object) # predict メソッド
{
object$predict
}
使用例
データ
#
# コーヒーデータ(回答者の属性など)の読み込み
COFFEE <- read.delim("COFFEE.txt", na.strings = ".")
#コーヒーデータ(広告への反応など)の読み込み
COFFEEcm <- read.delim("COFFEE-cm.txt", na.strings = ".")
#これらを結合して新しいデータフレームに
COFFEE <- data.frame(COFFEE, COFFEEcm)
attach(COFFEE)
d <- blikbnr
x <- data.frame( btaste1, bfraver1, bdesign1,
btaste2, bfraver2, bdesign2,
btaste3, bfraver3, bdesign3,
btaste4, bfraver4, bdesign4,
btaste5, bfraver5, bdesign5,
btaste6, bfraver6, bdesign6,
btaste7, bfraver7, bdesign7 )
ans8.2 <- mlm(d, x, constants=TRUE)
summary(ans8.2)
table(d, predict(ans8.2))
x <- data.frame( btaste1, bfraver1, bdesign1, aplesur1,
btaste2, bfraver2, bdesign2, aplesur2,
btaste3, bfraver3, bdesign3, aplesur3,
btaste4, bfraver4, bdesign4, aplesur4,
btaste5, bfraver5, bdesign5, aplesur5,
btaste6, bfraver6, bdesign6, aplesur6,
btaste7, bfraver7, bdesign7, aplesur7 )
ans8.4 <- mlm(d, x, constants=TRUE)
summary(ans8.4)
table(d, predict(ans8.4))
z <- data.frame(hindo) # 変数が 1 個のときでも data.frame を使うこと
ans8.6 <- mlm(d, x, z, constants=TRUE)
summary(ans8.6)
table(d, predict(ans8.6))
detach()
出力結果例
> # コーヒーデータ(回答者の属性など)の読み込み
> COFFEE <- read.delim("COFFEE.txt", na.strings = ".")
> #コーヒーデータ(広告への反応など)の読み込み
> COFFEEcm <- read.delim("COFFEE-cm.txt", na.strings = ".")
> #これらを結合して新しいデータフレームに
> COFFEE <- data.frame(COFFEE, COFFEEcm)
> attach(COFFEE)
> d <- blikbnr
> x <- data.frame( btaste1, bfraver1, bdesign1,
+ btaste2, bfraver2, bdesign2,
+ btaste3, bfraver3, bdesign3,
+ btaste4, bfraver4, bdesign4,
+ btaste5, bfraver5, bdesign5,
+ btaste6, bfraver6, bdesign6,
+ btaste7, bfraver7, bdesign7 )
> ans8.2 <- mlm(d, x, constants=TRUE)
> summary(ans8.2)
beta.estimated ses t P
btaste 0.60309 0.31009 1.945 0.055
bfraver 0.16744 0.28954 0.578 0.565
bdesign 0.61000 0.18954 3.218 0.002
Constant-1 -0.82574 0.40731 2.027 0.046
Constant-2 -0.39088 0.34058 1.148 0.254
Constant-3 -0.78882 0.52331 1.507 0.136
Constant-4 -0.26770 0.35458 0.755 0.452
Constant-5 -1.72052 0.62964 2.733 0.008
Constant-6 -1.97588 0.75332 2.623 0.010
LL = -126.1016
AIC = 270.2032
> table(d, predict(ans8.2))
d 1 2 3 4 7
1 2 5 0 1 1
2 0 21 0 2 6
3 0 3 0 2 0
4 0 3 1 5 7
5 0 3 0 0 0
6 0 1 0 0 1
7 0 4 0 1 21
>
> x <- data.frame( btaste1, bfraver1, bdesign1, aplesur1,
+ btaste2, bfraver2, bdesign2, aplesur2,
+ btaste3, bfraver3, bdesign3, aplesur3,
+ btaste4, bfraver4, bdesign4, aplesur4,
+ btaste5, bfraver5, bdesign5, aplesur5,
+ btaste6, bfraver6, bdesign6, aplesur6,
+ btaste7, bfraver7, bdesign7, aplesur7 )
> ans8.4 <- mlm(d, x, constants=TRUE)
> summary(ans8.4)
beta.estimated ses t P
btaste 0.57406 0.31043 1.849 0.068
bfraver 0.20588 0.28590 0.720 0.474
bdesign 0.56703 0.19257 2.945 0.004
aplesur 0.41502 0.17601 2.358 0.021
Constant-1 -0.78553 0.41262 1.904 0.061
Constant-2 0.20921 0.42214 0.496 0.622
Constant-3 -0.34342 0.55895 0.614 0.541
Constant-4 -0.33678 0.35880 0.939 0.351
Constant-5 -1.61925 0.63338 2.557 0.012
Constant-6 -1.50040 0.77897 1.926 0.058
LL = -123.2123
AIC = 266.4246
> table(d, predict(ans8.4))
d 1 2 3 4 7
1 3 4 0 1 1
2 0 19 0 3 7
3 0 3 0 2 0
4 0 5 1 4 6
5 1 2 0 0 0
6 0 1 0 0 1
7 0 3 0 2 21
>
> z <- data.frame(hindo) # 変数が 1 個のときでも data.frame を使うこと
> ans8.6 <- mlm(d, x, z, constants=TRUE)
> summary(ans8.6)
beta.estimated ses t P
btaste 0.492064 0.31908 1.542 0.127
bfraver 0.288856 0.29674 0.973 0.334
bdesign 0.626502 0.20240 3.095 0.003
aplesur 0.422931 0.17928 2.359 0.021
hindo-1 -0.354418 0.42709 0.830 0.409
hindo-2 -0.083622 0.29511 0.283 0.778
hindo-3 0.512617 0.50968 1.006 0.318
hindo-4 -0.822852 0.36821 2.235 0.028
hindo-5 0.323660 0.60486 0.535 0.594
hindo-6 0.893519 0.86130 1.037 0.303
Constant-1 0.248476 1.29336 0.192 0.848
Constant-2 0.439327 1.01733 0.432 0.667
Constant-3 -2.125203 1.93440 1.099 0.275
Constant-4 1.981968 1.06288 1.865 0.066
Constant-5 -2.645172 2.08739 1.267 0.209
Constant-6 -4.728288 3.47013 1.363 0.177
LL = -117.3403
AIC = 266.6807
> table(d, predict(ans8.6))
d 1 2 3 4 5 7
1 3 5 0 0 0 1
2 0 16 0 5 1 7
3 0 3 1 1 0 0
4 0 6 1 6 0 3
5 1 2 0 0 0 0
6 0 1 0 0 0 1
7 0 4 1 2 0 19
> detach()
直前のページへ戻る E-mail to Shigenobu AOKI
