DGPをいじりながらコンジョイントをやってみよう

パッケージの読み込み

library(dplyr)  # 使っても使わなくても読み込むパッケージ1

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)  # 使っても使わなくても読み込むパッケージ2
library(ade4)   # Factor型をダミー型に変換
library(cjoint)

## Loading required package: sandwich

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: ggplot2

## Loading required package: survey

## Loading required package: grid

## 
## Attaching package: 'survey'

## The following object is masked from 'package:graphics':
## 
##     dotchart

## cjoint: AMCE Estimator for Conjoint Experiments
## Version: 2.0.4

DGP

実験の設定は

回答者数: 10,000名
プロファイル数: 2
タスク数: 4
属性数: 3(A, B, C)
水準数: 2, 3, 4(属性1の場合、非復元抽出、他は復元抽出)

データ生成過程は

\[Y = \beta_0 + \beta_1 A_2 + \beta_2 B_2 + \beta_3 B_3 + \beta_4 C_2 + \beta_5 C_3 + \beta_6 C_4 + \varepsilon\]

だとする。各係数は以下のように設定する。

母数	係数
\(\beta_0\)	0
\(\beta_1\)	0.04651905
\(\beta_2\)	0.09523810
\(\beta_3\)	0.14285714
\(\beta_4\)	0.19047619
\(\beta_5\)	0.23809524
\(\beta_6\)	0.28571429
\(\sigma_{\varepsilon}\)	1

早速データセットを作ってみる。

n.Respondents <- 10000
n.Task        <- 4
n.Profile     <- 2
n.Obs         <- n.Respondents * n.Task * n.Profile

dgp.df <- data.frame(matrix(rep(NA, n.Obs * 7), 
                            nrow = n.Obs))
names(dgp.df) <- c("ID", "Task", "Profile", 
                   "Attr1", "Attr2", "Attr3", "Selected")

dgp.df$ID      <- rep(1:n.Respondents, each = n.Task * n.Profile)
dgp.df$Task    <- rep(rep(1:n.Task, each = n.Profile), n.Profile)
dgp.df$Profile <- rep(1:n.Profile, n.Respondents * n.Task)
dgp.df$Attr2   <- sample(c("B1", "B2", "B3"), n.Obs, replace = TRUE)
dgp.df$Attr3   <- sample(c("C1", "C2", "C3", "C4"), n.Obs, replace = TRUE)

Attr1 <- rep(NA, n.Obs)
for(i in seq(1, (n.Obs), by = 2)){
    Attr1[i:(i+1)] <- sample(c("A1", "A2"), 2, replace = FALSE)
}

dgp.df$Attr1   <- Attr1
dgp.df$Attr1   <- factor(dgp.df$Attr1)
dgp.df$Attr2   <- factor(dgp.df$Attr2)
dgp.df$Attr3   <- factor(dgp.df$Attr3)

# ここからselectedを計算
beta_vec <- (1:6)/sum(1:6)
temp.df  <- acm.disjonctif(dgp.df[, 4:6])[-c(1, 3, 6)]
e        <- rnorm(n.Obs, 0, 0.1)
Selected <- rep(NA, n.Obs)
for(i in 1:n.Obs){
    Selected[i] <- sum(temp.df[i,] * beta_vec)
}

Selected <- Selected + e

# 今のSelectedは選択確率。
# 別にこのまま計算してもいいが、リアリティーのために
# 選択確率が高い方に1をつける
for(i in seq(1, n.Obs, by = 2)){
    if(Selected[i] == Selected[i+1]){
        temp.result <- sample(0:1, 1, replace = FALSE)
    }else if(Selected[i] > Selected[i+1]){
        temp.result <- c(1, 0)
    }else{
        temp.result <- c(0, 1)
    }
    Selected[i:(i+1)] <- temp.result
}

dgp.df$Selected <- Selected

# これでデータセット作りは完了。
# 中身を見てみる
head(dgp.df)

##   ID Task Profile Attr1 Attr2 Attr3 Selected
## 1  1    1       1    A2    B1    C1        0
## 2  1    1       2    A1    B2    C3        1
## 3  1    2       1    A2    B1    C3        0
## 4  1    2       2    A1    B1    C3        1
## 5  1    3       1    A2    B1    C2        0
## 6  1    3       2    A1    B1    C3        1

tail(dgp.df)

##          ID Task Profile Attr1 Attr2 Attr3 Selected
## 79995 10000    2       1    A2    B3    C2        1
## 79996 10000    2       2    A1    B2    C1        0
## 79997 10000    3       1    A1    B1    C3        0
## 79998 10000    3       2    A2    B3    C4        1
## 79999 10000    4       1    A2    B2    C4        1
## 80000 10000    4       2    A1    B2    C1        0

コンジョイント分析

実際にcjointパッケージを用いて推定してみる。

Est1 <- amce(Selected ~ Attr1 + Attr2 + Attr3, data = dgp.df,
             respondent.id = "ID")
Est2 <- amce(Selected ~ Attr1 + Attr2 + Attr3, 
             data = filter(dgp.df, Profile == 1),
             respondent.id = "ID")

summary(Est1)$amce

##   Attribute Level  Estimate    Std. Err   z value Pr(>|z|)    
## 1     Attr1    A2 0.1671127 0.004166654  40.10716        0 ***
## 2     Attr2    B2 0.1717942 0.003819910  44.97337        0 ***
## 3     Attr2    B3 0.2550115 0.003758975  67.84069        0 ***
## 4     Attr3    C2 0.3234665 0.004212308  76.79080        0 ***
## 5     Attr3    C3 0.4083121 0.004109691  99.35349        0 ***
## 6     Attr3    C4 0.5032019 0.003952208 127.32169        0 ***

summary(Est2)$amce

##   Attribute Level  Estimate    Std. Err  z value      Pr(>|z|)    
## 1     Attr1    A2 0.1676327 0.004423632 37.89482  0.000000e+00 ***
## 2     Attr2    B2 0.1770276 0.005443581 32.52043 5.484219e-232 ***
## 3     Attr2    B3 0.2536654 0.005395813 47.01152  0.000000e+00 ***
## 4     Attr3    C2 0.3234017 0.006156375 52.53120  0.000000e+00 ***
## 5     Attr3    C3 0.4092543 0.005983901 68.39256  0.000000e+00 ***
## 6     Attr3    C4 0.5030460 0.005918553 84.99475  0.000000e+00 ***

結果の比較

ケースが半減するし、標準誤差を見ても意味がないので、係数だけみる。

属性	水準	プロファイル両方使用	プロファイル1のみ使用
A
	A2	0.167	0.168
B
	B2	0.172	0.177
	B3	0.255	0.254
C
	C2	0.323	0.323
	C3	0.408	0.409
	C4	0.503	0.503

どっちも係数はほぼ変わらないことが確認できた。

DGPと係数が異なるのはSelectedの値を選択確率でなく、バイナリー変数へ変換したため。（係数の比較という観点では問題ない）

DGPをいじりながらコンジョイントをやってみよう

Jaehyun Song

June 13, 2016

パッケージの読み込み

DGP

コンジョイント分析

結果の比較