-
Notifications
You must be signed in to change notification settings - Fork 0
/
examples.R
358 lines (296 loc) · 8.39 KB
/
examples.R
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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
set.seed(5678)
library(cbcTools)
library(idefix)
library(logitr)
library(jph)
library(tidyverse)
library(fastDummies)
# Define the attributes and levels
levels <- list(
brand = c("GM", "BMW", "Ferrari"),
price = c("20k", "40k", "100k")
)
profiles <- cbc_profiles(
brand = c("GM", "BMW", "Ferrari"),
price = c("20k", "40k", "100k")
)
design <- cbc_design(
profiles = profiles,
n_resp = 3, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 3 # Number of questions per respondent
)
cbc_balance(design)
cbc_overlap(design)
design <- cbc_design(
profiles = profiles,
n_resp = 30, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 3 # Number of questions per respondent
)
cbc_balance(design)
cbc_overlap(design)
# Make randomized design ----
design <- cbc_design(
profiles = profiles,
n_resp = 1200, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 3 # Number of questions per respondent
)
# Make D-efficient design with priors using {idefix} ----
cs <- Profiles(lvls = c(3, 3), coding = c("D", "D"))
mu <- c(1, 2, -1, -4)
sigma <- diag(length(mu))
M <- MASS::mvrnorm(n = 500, mu = mu, Sigma = sigma)
D <- Modfed(
cand.set = cs,
n.sets = 24,
n.alts = 3,
alt.cte = c(0, 0, 0),
par.draws = M,
max.iter = 25,
n.start = 5
)
DD <- Decode(des = D$design, lvl.names = levels, coding = c("D", "D"), n.alts = 3)
design_deff <- DD$design
# Repeat rows to match size of randomized design
n <- nrow(design) / nrow(design_deff)
design_deff <- design_deff[rep(seq_len(nrow(design_deff)), n), ]
row.names(design_deff) <- NULL
names(design_deff) <- c("brand", "price")
design_deff$brand <- factor(design_deff$brand, levels$brand)
design_deff$price <- factor(design_deff$price, levels$price)
design_deff$respID <- design$respID
design_deff$qID <- design$qID
design_deff$altID <- design$altID
design_deff$obsID <- design$obsID
# Compare balance
cbc_balance(design)
cbc_balance(design_deff)
# Sim data without interaction ----
data <- cbc_choices(
design = design,
obsID = "obsID",
priors = list(
price = c(-1, -4),
brand = c(1, 2)
)
)
data_deff <- cbc_choices(
design = design_deff,
obsID = "obsID",
priors = list(
price = c(-1, -4),
brand = c(1, 2)
)
)
# Sim data with interaction ----
data_int <- cbc_choices(
design = design,
obsID = "obsID",
priors = list(
price = c(-1, -4),
brand = c(1, 2),
`price*brand` = c(0.25, 0.5, 0.5, 1)
)
)
data_int_deff <- cbc_choices(
design = design_deff,
obsID = "obsID",
priors = list(
price = c(-1, -4),
brand = c(1, 2),
`price*brand` = c(0.25, 0.5, 0.5, 1)
)
)
# Power analysis on a design ----
# Estimate models with different sample sizes
results <- cbc_power(
nbreaks = 10,
n_q = 3,
data = data,
pars = c("price", "brand"),
outcome = "choice",
obsID = "obsID"
)
results_deff <- cbc_power(
nbreaks = 10,
n_q = 3,
data = data_deff,
pars = c("price", "brand"),
outcome = "choice",
obsID = "obsID"
)
results_int <- cbc_power(
nbreaks = 10,
n_q = 3,
data = data_int,
pars = c("price", "brand", "price*brand"),
outcome = "choice",
obsID = "obsID"
)
results_int_deff <- cbc_power(
nbreaks = 10,
n_q = 3,
data = data_int_deff,
pars = c("price", "brand", "price*brand"),
outcome = "choice",
obsID = "obsID"
)
# Visualize
results %>%
mutate(type = "Random") %>%
rbind(results_deff %>% mutate(type = "D-efficient")) %>%
ggplot() +
geom_hline(yintercept = 0.05, color = "red", linetype = 2) +
geom_point(
aes(x = sampleSize, y = se, color = type),
size = 1.8
) +
facet_wrap(vars(coef), nrow = 1) +
scale_x_continuous(limits = c(0, 1250)) +
scale_y_continuous(limits = c(0, 0.25)) +
expand_limits(y = 0) +
theme_bw(base_size = 14) +
theme(panel.grid.minor = element_blank()) +
labs(
color = "Design",
x = "Sample size",
y = "Standard error"
)
ggsave(file.path("images", "design_compare.png"), width = 12, height = 3.5)
results_int %>%
mutate(type = "Random") %>%
rbind(results_deff %>% mutate(type = "D-efficient")) %>%
mutate(
coef = fct_relevel(coef, c(
"price40k", "price100k", "brandBMW", "brandFerrari",
"price40k:brandBMW", "price100k:brandBMW",
"price40k:brandFerrari", "price100k:brandFerrari"))) %>%
ggplot() +
geom_hline(yintercept = 0.05, color = "red", linetype = 2) +
geom_point(
aes(x = sampleSize, y = se, color = type),
size = 1.8
) +
facet_wrap(vars(coef), nrow = 2) +
scale_x_continuous(limits = c(0, 1250)) +
scale_y_continuous(limits = c(0, 1)) +
expand_limits(y = 0) +
theme_bw(base_size = 14) +
theme(panel.grid.minor = element_blank()) +
labs(
color = "Design",
x = "Sample size",
y = "Standard error"
)
ggsave(file.path("images", "design_compare_int.png"), width = 12, height = 6)
# Apple conjoint example ----
# Define the attributes and levels
levels <- list(
price = seq(1, 4, 0.5), # $ per pound
type = c("Fuji", "Gala", "Honeycrisp"),
freshness = c("Excellent", "Average", "Poor")
)
# Generate all all possible profiles
profiles <- cbc_profiles(levels)
# Make a randomized survey design
design_rand <- cbc_design(
profiles = profiles,
n_resp = 300, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 6 # Number of questions per respondent
)
# Make a randomized survey design with a "no choice" option
design_rand_nochoice <- cbc_design(
profiles = profiles,
n_resp = 300, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 6, # Number of questions per respondent
no_choice = TRUE
)
# Make randomized labeled survey design with each "type" appearing in each
# choice question
design_rand_labeled <- cbc_design(
profiles = profiles,
n_resp = 300, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 6, # Number of questions per respondent
label = "type"
)
# Make randomized labeled survey design with each "type" appearing in each
# choice question and with a "no choice" option
design_rand_labeled_nochoice <- cbc_design(
profiles = profiles,
n_resp = 300, # Number of respondents
n_alts = 3, # Number of alternatives per question
n_q = 6, # Number of questions per respondent
no_choice = TRUE,
label = "type"
)
# Inspect survey design ----
# View the attribute balance and pairwise attribute balance
cbc_balance(design_rand)
# View the amount of overlap in a design
cbc_overlap(design_rand)
# Simulate choices ----
# Simulate random choices for a survey design
data_rand <- cbc_choices(
design = design_rand,
obsID = "obsID"
)
# Simulate choices based on a utility model with the following parameters:
# - 1 continuous "price" parameter
# - 2 categorical parameters for "type" (first level is reference)
# - 2 categorical parameters for "freshness" (first level is reference)
data_prior <- cbc_choices(
design = design_rand,
obsID = "obsID",
priors = list(
price = -0.1,
type = c(0.1, 0.2),
freshness = c(0.1, -0.2)
)
)
# Simulate choices based on a utility model with the following parameters:
# - 1 continuous "price" parameter
# - 2 categorical parameters for "type" (first level is reference)
# - 2 categorical parameters for "freshness" (first level is reference)
# - 2 interaction parameters between "price" and "type"
data_prior_int <- cbc_choices(
design = design_rand,
obsID = "obsID",
priors = list(
price = -0.1,
type = c(0.1, 0.2),
freshness = c(0.1, -0.2),
`price*type` = c(0.1, 0.5)
)
)
# Simulate choices based on a utility model with the following parameters:
# - 1 continuous "price" parameter
# - 2 random normal discrete parameters for "type" (first level is reference)
# - 2 categorical parameters for "freshness" (first level is reference)
data_prior_mixed <- cbc_choices(
design = design_rand,
obsID = "obsID",
priors = list(
price = -0.1,
type = randN(mu = c(0.1, 0.2), sigma = c(0.5, 1)),
freshness = c(0.1, -0.2)
)
)
# Estimate models with different sample sizes
results <- cbc_power(
nbreaks = 10,
n_q = 6,
data = data_rand,
pars = c("price", "type", "freshness"),
outcome = "choice",
obsID = "obsID"
)
# Preview
head(results)
tail(results)
# Visualize
plot(results)