多項ロジットモデル     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

Made with Macintosh