R

選択型コンジョイント分析１

survival , support.CEs , conjoint , knitr パッケージ

(使用したサンプルデータ)
選択型コンジョイント分析(pdf) で分析しているデータ

条件付きロジット・モデル

質問表の作成手順
1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(survival)
library(support.CEs)
library(conjoint)
library(knitr)
#
expand<-expand.grid(
	HACCP=c("なし","あり"),
	エコ=c("なし","あり"),
	賞味=c("６日間","７日間"),
	価格=c("145","150","155","160") )
#因子の直交表を作成
design<-caFactorialDesign(data=expand,type="orthogonal")
#因子の直交表を表示
design

HACCP エコ 賞味 価格
2 あり なし ６日間 145
7 なし あり ７日間 145
11 なし あり ６日間 150
14 あり なし ７日間 150
20 あり あり ６日間 155
21 なし なし ７日間 155
25 なし なし ６日間 160
32 あり あり ７日間 160 

1
2
3
4
# コピーを作る
design2<-design
#
rownames(design) == rownames(design2)

TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

乱数を発生させて並び替える
1
2
3
4
5
6
7
8
# 乱数を発生させて並び替える
design <- transform(design,r1=runif(8)) 
design_sort <- design[order(design$r1),] 
# 
design2 <- transform(design2,r1=runif(8)) 
design2_sort <- design2[order(design2$r1),]
# 確認
rownames(design_sort) == rownames(design2_sort)

FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

選択１，２がすべて異なるのでこれをもとに質問を作成する。

TRUE があったら、再度、乱数を発生させて並び替える

1
design_sort ; design2_sort

HACCP エコ 賞味 価格 r1
32 あり あり ７日間 160 0.1268605
20 あり あり ６日間 155 0.1501261
11 なし あり ６日間 150 0.2241221
2 あり なし ６日間 145 0.4913059
21 なし なし ７日間 155 0.5257258
14 あり なし ７日間 150 0.7125287
25 なし なし ６日間 160 0.8284758
7 なし あり ７日間 145 0.9818714

HACCP エコ 賞味 価格 r1
11 なし あり ６日間 150 0.09701797
25 なし なし ６日間 160 0.13295188
14 あり なし ７日間 150 0.16370145
7 なし あり ７日間 145 0.24600342
20 あり あり ６日間 155 0.66435520
21 なし なし ７日間 155 0.67742812
2 あり なし ６日間 145 0.81968328
32 あり あり ７日間 160 0.97952398 

1
2
3
4
5
6
7
#因子の直交表をコード化
code1 <- caEncodedDesign(design_sort[,1:4]) -1
code2 <- caEncodedDesign(design2_sort[,1:4]) -1
#
#コード化された直交表の表示
#
code1　; code2

HACCP エコ 賞味 価格
32 1 1 1 3
20 1 1 0 2
11 0 1 0 1
2 1 0 0 0
21 0 0 1 2
14 1 0 1 1
25 0 0 0 3
7 0 1 1 0

HACCP エコ 賞味 価格
11 0 1 0 1
25 0 0 0 3
14 1 0 1 1
7 0 1 1 0
20 1 1 0 2
21 0 0 1 2
2 1 0 0 0
32 1 1 1 3 

1
2
3
4
5
# 価格を数値に変換する場合 : as.numeric(as.vector()) とすること。
code1 <- data.frame(caEncodedDesign(design_sort[,1:3]) -1,価格=as.numeric(as.vector(design_sort[,4])) )
code2 <- data.frame(caEncodedDesign(design2_sort[,1:3]) -1,価格=as.numeric(as.vector(design2_sort[,4])) )
#
code1 ; code2

HACCP エコ 賞味 価格
32 1 1 1 160
20 1 1 0 155
11 0 1 0 150
2 1 0 0 145
21 0 0 1 155
14 1 0 1 150
25 0 0 0 160
7 0 1 1 145

HACCP エコ 賞味 価格
11 0 1 0 150
25 0 0 0 160
14 1 0 1 150
7 0 1 1 145
20 1 1 0 155
21 0 0 1 155
2 1 0 0 145
32 1 1 1 160

調査を行う。

data.frame : select1 , select2 を作成しているが、上の手順で作成した質問表を使ったならcode1 , code2 に「定数」を加えればよい。
select3 は「どちらも選ばない」が選択肢にある場合に必要。kaito=3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 回答者数
sn<-10
# 質問数
qn<-8
#
select1<-data.frame(
	HACCP=c(1,1,0,0,0,1,0,1),
	エコ=c(0,0,0,1,0,1,1,1),
	賞味=c(1,0,1,1,0,0,0,1),
	価格=c(155,150,160,150,145,160,155,145),
	定数 = rep(0,qn) )
#
select2<-data.frame(
	HACCP=c(0,1,0,0,0,1,1,1),
	エコ=c(1,1,0,0,1,0,1,0),
	賞味=c(1,0,0,1,0,0,1,1),
	価格=c(150,160,145,160,155,150,145,155),
	定数 = rep(0,qn) )
#
select3<-data.frame(
	HACCP = rep(0 ,qn) ,
	エコ = rep(0 ,qn) ,
	賞味 = rep(0 ,qn) ,
	価格 = rep(0 ,qn) ,
	定数 = rep(1 ,qn) )
#
kaito<-c(2,1,2,1,1,2,2,1,1,2,3,3,3,1,2,1,2,2,3,1,2,1,1,1,2,3,1,1,3,3,2,1,1,2,3,3,3,
	1,2,1,2,2,3,1,2,1,1,1,1,2,3,3,3,1,2,1,2,2,3,1,2,1,1,1,2,3,1,1,3,3,2,1,2,2,3,1,2,1,1,1)

分析のための整形

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
select<-NULL
for (i in 1 : sn ){
	for (j in 1 : qn ){
		select<-rbind(select,select1[j,],select2[j,],select3[j,])
		}
	}
#
回答<-NULL
for (i in 1 : length(kaito) ){
	num<-c(0,0,0)
	 switch(kaito[i] ,               # switch(文字列,
	   "1" = (num[1] = 1),   #  "1" のときに実行
	   "2" = (num[2] = 1),   #  "2" のときに実行
	   "3" = (num[3] = 1),   #  "3" のときに実行 (「どちらも選ばない」が選択肢にある場合に必要)
 		)  
	回答 = c(回答,num)
	}
select$回答 <- 回答
#
select$str<-paste0(rep(1:sn,each=qn*3),rep(0,sn*qn),rep(rep(1:qn,each=3),sn))
#
head(select) ; tail(select)

HACCP エコ 賞味 価格 定数 回答 str
1 1 0 1 155 0 0 101
2 0 1 1 150 0 1 101
3 0 0 0 0 1 0 101
21 1 0 0 150 0 1 102
22 1 1 0 160 0 0 102
23 0 0 0 0 1 0 102

HACCP エコ 賞味 価格 定数 回答 str
727 0 1 0 155 0 1 1007
728 1 1 1 145 0 0 1007
729 0 0 0 0 1 0 1007
827 1 1 1 145 0 1 1008
828 1 0 1 155 0 0 1008
829 0 0 0 0 1 0 1008

分析

1
2
out <- clogit( 回答 ~ 定数 + HACCP + エコ + 賞味 + 価格 + strata(str),data=select)
out

coef exp(coef) se(coef) z p
定数 -1.8819 0.1523 5.4111 -0.35 0.728
HACCP 0.9399 2.5598 0.4541 2.07 0.038
エコ 2.2603 9.5863 0.3998 5.65 1.6e-08
賞味 0.5721 1.7720 0.4983 1.15 0.251
価格 -0.0230 0.9772 0.0349 -0.66 0.510

Likelihood ratio test=54.1 on 5 df, p=2.01e-10
n= 240, number of events= 80 

1
2
# 対数尤度
out$loglik

[1] -87.88898 -60.84726

1
2
# McFaddenの決定係数
gofm(out)

Adjusted rho-squared = 0.2507905 > 0.2 —— 自由度調整済み疑似決定係数の値が0.2を超えると適合度が良いと判断される

限界支払意思額
1
2
3
4
5
6
7
8
mwtp1 <- mwtp(
output = out ,
monetary.variables = c("価格"),
nonmonetary.variables =
c("HACCP","エコ","賞味"),
seed = 1)
#
mwtp1

MWTP 2.5% 97.5%
HACCP 40.84 -304.16 347.16
エコ 98.20 -793.43 826.89
賞味 24.86 -237.01 221.51

method = Krinsky and Robb 

1
out$coeff

定数 HACCP エコ 賞味 価格
-1.88185292 0.93992097 2.26033843 0.57208637 -0.02301687

参考にしたテキストの　- VI -40 - は、確定効用に「選択肢固有定数項」をたしているが、
「選択肢固有定数項」は、「どちらも買わない」選択肢のみに設定されているのでたすのは間違いなのでは。
（　- VI -34 - では　牛乳１の確定効用、牛乳２の確定効用の計算では「選択肢固有定数項」をたしていない。）

確定効用
1
2
3
4
5
6
7
8
9
価格<-c(seq(100,200,10))
#
h0h1<-c(0,out$coeff[2])
e0e1<-c(0,out$coeff[3])
s0s1<-c(0,out$coeff[4])
#
exg<-expand.grid(h0h1,e0e1,s0s1) 
colnames(exg)<-c("HACCP","エコ","賞味") 
exg

HACCP エコ 賞味
1 0.000000 0.000000 0.0000000
2 0.939921 0.000000 0.0000000
3 0.000000 2.260338 0.0000000
4 0.939921 2.260338 0.0000000
5 0.000000 0.000000 0.5720864
6 0.939921 0.000000 0.5720864
7 0.000000 2.260338 0.5720864
8 0.939921 2.260338 0.5720864

(例)
HACCPラベル :あり
エコラベル :あり
賞味期限:7日間 

1
2
3
4
確定効用<-sum(exg[8,]) + out$coeff[5] * 価格
# 確定効用<-out$coeff[2] + out$coeff[3] + out$coeff[4] + out$coeff[5] * 価格　#でもよい
#
kable(data.frame(価格,確定効用=round(確定効用 ,3)) )
価格 確定効用
100 1.471
110 1.240
120 1.010
130 0.780
140 0.550
150 0.320
160 0.090
170 -0.141
180 -0.371
190 -0.601
200 -0.831
購入確率 EXP()/(EXP()+EXP(0))
1
2
3
購入確率 <- exp(確定効用)/(exp(確定効用) + exp(0))
#
kable(data.frame(価格,購入確率=round(購入確率 ,3)) )
価格 購入確率
100 0.813
110 0.776
120 0.733
130 0.686
140 0.634
150 0.579
160 0.522
170 0.465
180 0.408
190 0.354
200 0.303 
1
2
3
# png("clogit01.png")
plot(購入確率,価格,type="o",pch=16,las=1,main="購入確率曲線",ylab="価格( 円 )")
# dev.off()

需要量

(仮定)
二人以上世帯における牛乳の年間購入量 93.13 リットル
二人以上世帯数 267,270 

1
2
3
需要量<- 93.13 * 267270 * 購入確率/1000
#
kable(data.frame(価格,需要量=round(需要量 ,0)) )
価格 需要量
100 20240
110 19307
120 18247
130 17068
140 15784
150 14419
160 13003
170 11572
180 10165
190 8815
200 7553 
1
2
3
4
# png("clogit02.png")
par(family="TakaoMincho")
plot(需要量,価格,type="o",pch=16,las=1,main="需要曲線",xlab="需要量( kl )",ylab="価格( 円 )")
# dev.off()