濱岡 豊｢8章 マーケティング ブランド選択行動｣
『経済・経営のための統計学』有斐閣
演習用プログラム
実行したい部分をコピーして､Rのコンソールにペーストすると実行されます｡
p.219
options(width=120,digits=4)
#あらかじめ これらのファイルがある場所をワーキングディレクトリworking directory指定すること｡
#コーヒーデータ(回答者の属性など)の読み込み
COFFEE<-read.delim("COFFEE.txt",na.strings = ".")
#コーヒーデータ(広告への反応など)の読み込み
COFFEEcm<-read.delim("COFFEE-cm.txt",na.strings = ".")
#これらを結合して新しいデータフレームに
COFFEE<-data.frame(COFFEE,COFFEEcm)
save(COFFEE,file="COFFEE.rda")
attach(COFFEE)
summary(COFFEE)
p.220
#もっとも好きな広告
table(blikcmr)
#もっとも好きなブランド
table(blikbnr)
#もっとも買いたいブランド
table(binpbnr)
#ブランド選択ダミーの定義(もっとも好きなブランドとして選ばれたものについては1)
COFFEE$d1<-as.numeric(blikbnr==1)
COFFEE$d2<-as.numeric(blikbnr==2)
COFFEE$d3<-as.numeric(blikbnr==3)
COFFEE$d4<-as.numeric(blikbnr==4)
COFFEE$d5<-as.numeric(blikbnr==5)
COFFEE$d6<-as.numeric(blikbnr==6)
COFFEE$d7<-as.numeric(blikbnr==7)
attach(COFFEE)
#----------------------------Rの出力8-1 各ブランドの3つの属性の平均値と､選択されたシェア
mean(data.frame(btaste1,bfraver1,bdesign1,d1))
mean(data.frame(btaste2,bfraver2,bdesign2,d2))
mean(data.frame(btaste3,bfraver3,bdesign3,d3))
mean(data.frame(btaste4,bfraver4,bdesign4,d4))
mean(data.frame(btaste5,bfraver5,bdesign5,d5))
mean(data.frame(btaste6,bfraver6,bdesign6,d6))
mean(data.frame(btaste7,bfraver7,bdesign7,d7))
#----------------------------Rの出力8-2-ブランドへの属性評価とブランド定数で説明
#大阪大学経済学研究科 里村助教授のページを参考にロジット分析ルーチン作成
#http://www2.econ.osaka-u.ac.jp/~satomura/
#R では添え字を使うと処理が遅くなるので添え字無しに変更&ここでのデータに適用できるように定式化しただけ
mlogit <- function(beta) {
#3つのβについて初期値
beta1 <- beta[1]
beta2 <- beta[2]
beta3 <- beta[3]
#ブランド定数の初期値
bc01<- beta[4]
bc02<- beta[5]
bc03<- beta[6]
bc04<- beta[7]
bc05<- beta[8]
bc06<- beta[9]
#確定的効用
V1<- beta1*btaste1+beta2*bfraver1+beta3*bdesign1+bc01
V2<- beta1*btaste2+beta2*bfraver2+beta3*bdesign2+bc02
V3<- beta1*btaste3+beta2*bfraver3+beta3*bdesign3+bc03
V4<- beta1* btaste4+beta2*bfraver4+beta3*bdesign4+bc04
V5<- beta1* btaste5+beta2*bfraver5+beta3*bdesign5+bc05
V6<- beta1* btaste6+beta2*bfraver6+beta3*bdesign6+bc06
V7<- beta1* btaste7+beta2*bfraver7+beta3*bdesign7
#対数尤度関数
denom<-exp(V1)+exp(V2)+exp(V3)+exp(V4)+exp(V5)+exp(V6)+exp(V7)
LLL<-d1*V1 +d2*V2+d3*V3+d4*V4+d5*V5+d6*V6+d7*V7 -log(denom)
LL<- sum(LLL)
return(LL)
}
#-------------------------------Rの出力8-3
#初期値の設定 多項ロジットは初期値に関係なく推定できるので0にしてしまう｡
beta<-c(0,0,0,0,0,0,0,0,0)
#optimルーチンによって計算｡結果をresに
res<-optim(beta,mlogit, method = "BFGS", hessian = TRUE, control=list(fnscale=-1))
res
#最大対数尤度
LL<-res$value
#推定値
beta.estimated<-res$par
#赤池の情報量基準 AIC=-2(最大対数尤度- 推定したパラメータ数)
AIC<- -2*(LL-length(beta.estimated))
#標準誤差=-1*ヘシアン行列の逆行列 の対角項
ses<-sqrt(diag(solve(-res$hessian) ))
#t値
t<-beta.estimated/ses
#これらを一覧
data.frame(beta.estimated,ses,t)
LL
AIC
#--------------------------------------------------Rの出力8.4 感情的反応を入れたロジットモデル
mlogit <- function(beta) {
beta1 <- beta[1]
beta2 <- beta[2]
beta3 <- beta[3]
beta4 <- beta[4]
#Brand constant
bc01<- beta[5]
bc02<- beta[6]
bc03<- beta[7]
bc04<- beta[8]
bc05<- beta[9]
bc06<- beta[10]
#Vij
V1<- beta1*btaste1+beta2*bfraver1+beta3*bdesign1+beta4*aplesur1+bc01
V2<- beta1*btaste2+beta2*bfraver2+beta3*bdesign2+beta4*aplesur2+bc02
V3<- beta1*btaste3+beta2*bfraver3+beta3*bdesign3+beta4*aplesur3+bc03
V4<- beta1*btaste4+beta2*bfraver4+beta3*bdesign4+beta4*aplesur4+bc04
V5<- beta1*btaste5+beta2*bfraver5+beta3*bdesign5+beta4*aplesur5+bc05
V6<- beta1*btaste6+beta2*bfraver6+beta3*bdesign6+beta4*aplesur6+bc06
V7<- beta1*btaste7+beta2*bfraver7+beta3*bdesign7+beta4*aplesur7
#LL
denom<-exp(V1)+exp(V2)+exp(V3)+exp(V4)+exp(V5)+exp(V6)+exp(V7)
LLL<-d1*V1 +d2*V2+d3*V3+d4*V4+d5*V5+d6*V6+d7*V7 -log(denom)
LL<- sum(LLL)
return(LL)
}
#-------------------------------Rの出力8-5
#初期値の設定
beta<-c(0,0,0,0,0,0,0,0,0,0)
res<-optim(beta,mlogit, method = "BFGS", hessian = TRUE, control=list(fnscale=-1))
res
LL<-res$value
beta.estimated<-res$par
AIC<- -2*(LL-length(beta.estimated)) #AIC=-2(LL- # of parameters)
ses<-sqrt(diag(solve(-res$hessian)) )
t<-beta.estimated/ses
data.frame(beta.estimated,ses,t)
LL
AIC
#--------------------------------------------------Rの出力8-6 飲用頻度を入れたモデル。
mlogit <- function(beta) {
beta1 <- beta[1]
beta2 <- beta[2]
beta3 <- beta[3]
beta4 <- beta[4]
#ブランド乗数
bc01<- beta[5]
bc02<- beta[6]
bc03<- beta[7]
bc04<- beta[8]
bc05<- beta[9]
bc06<- beta[10]
#ブランド×飲用頻度(の少なさ)
sbc01<- beta[11]
sbc02<- beta[12]
sbc03<- beta[13]
sbc04<- beta[14]
sbc05<- beta[15]
sbc06<- beta[16]
#確定的効用
V1<- beta1*btaste1+beta2*bfraver1+beta3*bdesign1+beta4*aplesur1+bc01+sbc01*hindo
V2<- beta1*btaste2+beta2*bfraver2+beta3*bdesign2+beta4*aplesur2+bc02+sbc02*hindo
V3<- beta1*btaste3+beta2*bfraver3+beta3*bdesign3+beta4*aplesur3+bc03+sbc03*hindo
V4<- beta1*btaste4+beta2*bfraver4+beta3*bdesign4+beta4*aplesur4+bc04+sbc04*hindo
V5<- beta1*btaste5+beta2*bfraver5+beta3*bdesign5+beta4*aplesur5+bc05+sbc05*hindo
V6<- beta1*btaste6+beta2*bfraver6+beta3*bdesign6+beta4*aplesur6+bc06+sbc06*hindo
V7<- beta1*btaste7+beta2*bfraver7+beta3*bdesign7+beta4*aplesur7
#対数尤度関数
denom<-exp(V1)+exp(V2)+exp(V3)+exp(V4)+exp(V5)+exp(V6)+exp(V7)
LLL<-d1*V1 +d2*V2+d3*V3+d4*V4+d5*V5+d6*V6+d7*V7 -log(denom)
LL<- sum(LLL)
return(LL)
}
#-------------------------------Rの出力8-5
#初期値の設定
beta<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
res<-optim(beta,mlogit, method = "BFGS", hessian = TRUE, control=list(fnscale=-1))
res
LL<-res$value
beta.estimated<-res$par
AIC<- -2*(LL-length(beta.estimated))
ses<-sqrt(diag(solve(-res$hessian)) )
t<-beta.estimated/ses
data.frame(beta.estimated,ses,t)
LL
AIC