目的 牧厚志ら「経済・経営のための統計学」有斐閣アルマの第 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()