Data
Choice data should be provided in a `long’ format tibble or
data.frame, i.e. one row per alternative, choice task and
respondent.
It should contain the following columns
id (integer; respondent identifier)task (integer; task number)alt (integer; alternative #no within task)x (double; quantity purchased)p (double; price)- attributes defining the choice alternatives (factor, ordered)
and/or
- continuous attributes defining the choice alternatives
(numeric)
By default, it is assumed that respondents were able to choose the
‘outside good’, i.e. there is a ‘no-choice’ option. The no-choice option
should not be explicitly stated, i.e. there is
no “no-choice-level” in any of the attributes. This differs common
practice in discrete choice modeling, however, the
implicit outside good is consistent with both
volumetric and discrete choice models based on
economic assumptions.
Dummy-coding is performed automatically to ensure consistency between
estimation and prediction, and to set up screening models in a
consistent manner. Categorical attributes should therefore be provided
as factor or ordered, and not be
dummy-coded.
Let’s consider an example dataset from a volumetric conjoint about
frozen pizza:
data(pizza)
head(pizza)
#> # A tibble: 6 × 11
#> id task alt x p brand size crust topping coverage cheese
#> <int> <int> <int> <dbl> <dbl> <fct> <fct> <fct> <fct> <fct> <fct>
#> 1 2 1 1 0 4.5 Tony forOne StufCr Pepperoni densetop realcheese
#> 2 2 1 2 0 4.5 DiGi ForTwo Thin Surp densetop NoInfo
#> 3 2 1 3 0 3 Fresc ForTwo TrCr Veg ModCover NoInfo
#> 4 2 1 4 0 1.75 Priv forOne RisCr HI ModCover realcheese
#> 5 2 1 5 2 3.25 Tomb forOne TrCr Cheese ModCover NoInfo
#> 6 2 1 6 0 3.75 RedBa ForTwo StufCr PepSauHam densetop realcheese
There are n_distinct(pizza$id) unique despondents who
respond to pizza$task %>% n_distinct() choice tasks with
n_distinct(pizza$alt) alternatives each. The convenience
functionec_summarize_attrlvls provides a quick glance at
the attributes and levels (for categorical attributes). The function
uses standard ‘tidyverse’ functions to generate that summary.
pizza %>% ec_summarize_attrlvls
#> # A tibble: 6 × 2
#> attribute levels
#> <chr> <chr>
#> 1 brand DiGi, Fresc, Priv, RedBa, Tomb, Tony
#> 2 size forOne, ForTwo
#> 3 crust Thin, RisCr, StufCr, TrCr
#> 4 topping Pepperoni, Cheese, HI, PepSauHam, Surp, Veg
#> 5 coverage densetop, ModCover
#> 6 cheese NoInfo, realcheese
Looking across choice tasks, we can see that sometimes respondents
choose 0 quantity of all alternatives in a choice task. In this case,
their entire budget is allocated towards the ‘outside good’, which is
not explicitly stated in the data. Even when they choose any positive
quantity of alternatives, they would still consume some amount of
outside good. Volumetric demand models regularly assume that respondents
never consume their entire category budget, there will always some
amount left (and if it is only a fraction of a cent).
pizza %>%
group_by(id, task) %>%
summarise(x_total = sum(x), .groups = "keep") %>%
group_by(x_total) %>% count(x_total)
#> # A tibble: 22 × 2
#> # Groups: x_total [22]
#> x_total n
#> <dbl> <int>
#> 1 0 282
#> 2 1 316
#> 3 2 330
#> 4 3 207
#> 5 4 233
#> 6 5 102
#> 7 6 406
#> 8 7 36
#> 9 8 52
#> 10 9 15
#> # … with 12 more rows
The long “tidy” data structure makes it easy to work with choice
data. Using just a few lines of code, we can visualize the distribution
of quantity demanded per choice task, on average. From this, we can
clearly see that respondents frequently buy more than one unit of frozen
pizza, which suggest that indeed a volumetric choice experiment is more
realistic than a discrete choice one.
pizza %>%
group_by(id, task) %>%
summarise(sum_x=sum(x), .groups = "keep") %>%
group_by(id) %>%
summarise(mean_x=mean(sum_x), .groups = "keep") %>%
ggplot(aes(x=mean_x)) +
geom_density(fill="lightgrey", linewidth=1) +
geom_vline(xintercept=1, linewidth=2, color="darkred") +
theme_minimal() +
xlab("Average units per task")
Holdout
For hold-out validation, we keep 1 task per respondent. In v-fold
cross-validation, this is done several times. However, each re-run of
the model may take a while. For this example, we only use 1 set of
holdout tasks. Hold-out evaluations results may vary slightly between
publications that discuss this dataset.
#randomly assign hold-out group, use seed for reproducible plot
set.seed(1.2335252)
pizza_ho_tasks=
pizza %>%
distinct(id,task) %>%
mutate(id=as.integer(id))%>%
group_by(id) %>%
summarise(task=sample(task,1), .groups = "keep")
set.seed(NULL)
#calibration data
pizza_cal= pizza %>% mutate(id=as.integer(id)) %>%
anti_join(pizza_ho_tasks, by=c('id','task'))
#'hold-out' data
pizza_ho= pizza %>% mutate(id=as.integer(id)) %>%
semi_join(pizza_ho_tasks, by=c('id','task'))
Estimation
Estimate both models using at least 200,000 draws. Saving
each 50th or 100th draw is sufficient. The vd_est_vdm fits
the compensatory volumetric demand model, while
vd_est_vdm_screen fits the model with attribute-based
conjunctive screening. Using the error_dist argument, the
type of error distribution can be specified. While KHKA 2022
assume Normal-distributed errors, here we assume Extreme Value Type 1
errors.
#compensatory
out_pizza_cal = pizza_cal %>% vd_est_vdm(R=200000, keep=50, error_dist = "EV1")
dir.create("draws")
save(out_pizza_cal,file='draws/out_pizza_cal.rdata')
#conjunctive screening
out_pizza_screening_cal = pizza_cal %>% vd_est_vdm_screen(R=200000, keep=50, error_dist = "EV1")
save(out_pizza_screening_cal,file='draws/out_pizza_screening_cal.rdata')
Diagnostics
Quick check of convergence, or stationarity of the traceplots.
Compensatory:
out_pizza_cal %>% ec_trace_MU(burnin = 100)
Conjunctive Screening:
out_pizza_screening_cal %>% ec_trace_MU(burnin = 100)
Distribution of Log-likehoods to check for outlier respondents:
vd_LL_vdm(out_pizza_cal ,pizza_cal, fromdraw = 3000) %>%
apply(1,mean) %>% tibble(LL=.) %>%
ggplot(aes(x=LL)) +
geom_density(fill="lightgrey", linewidth=1) + theme_minimal() +
xlab("Individual average Log-Likelihood")
Here the left tail looks fine, and there do not appear to be
extremely “bad” responders.
Demand Curves
Elasticities in volumetric demand models are not constant. To better
understand implications of pricing decisions, demand curves can be a
helpful tool.
The ec_demcurve function can be applied to all demand
simulators. Based of an initial market scenario, it generates a series
of scenarios where the price of a focal product is changed over an
interval. It then runs the demand simulator several times, and we can
use the output to draw a demand curve.
Pre-simulating error terms using ec_gen_err_ev1 helps to
smooth these demand curves. It generates one error term per draw for
each of the alternatives and tasks in the corresponding design
dataset.
Define base case
We have to start by defining a base case, i.e. a
scenario that we want to start with.
#define a set of offerings
testm_pizza =
tibble(
id=1L, task=1L, alt=1:6,
brand = c("DiGi", "Fresc", "Priv", "RedBa", "Tomb", "Tony"),
size = "forOne",
crust = "Thin",
topping = "Veg",
coverage= "ModCover",
cheese = "NoInfo",
p = c(3.5,3,2,2,2,1.5)) %>%
prep_newprediction(pizza)
#for this experiment, offer this assortment to each of the respondents
testmarket=
tibble(
id = rep(seq_len(n_distinct(pizza$id)),each=nrow(testm_pizza)),
task = 1,
alt = rep(1:nrow(testm_pizza),n_distinct(pizza$id))) %>%
bind_cols(
testm_pizza[rep(1:nrow(testm_pizza),n_distinct(pizza$id)),-(1:3)])
Demand curves
The output of ec_demcurve is a list, containing demand
summaries for each product under the different price scenarios. The list
can be stacked into a data.frame, which can be used to generate
plots.
#Define focal brand for analysis
focal_alternatives <-
testmarket %>% transmute(focal=brand=='Priv') %>% pull(focal)
#pre-sim error terms
eps_not <- testmarket %>% ec_gen_err_ev1(out_pizza_cal, 55667)
#demand curve compensatory
vd_demc_comp <-
testmarket %>%
ec_demcurve(focal_alternatives,
seq(0.5,1.5,,9),
vd_dem_vdm,
out_pizza_cal,
eps_not)
#demand curve conjunctive screening
vd_demc_screen <-
testmarket %>%
ec_demcurve(focal_alternatives,
seq(0.5,1.5,,9),
vd_dem_vdm_screen,
out_pizza_screening_cal,
eps_not)
#combine demand curves from both models
vd_outputs=rbind(
vd_demc_comp %>% do.call('rbind',.) %>% bind_cols(model='comp') %>% bind_cols(demand='volumetric'),
vd_demc_screen %>% do.call('rbind',.) %>% bind_cols(model='screenpr') %>% bind_cols(demand='volumetric'))
#plot
vd_outputs%>%
filter(brand=="Priv") %>%
ggplot(aes(x=scenario, y=`E(demand)`, color=model)) + geom_line(size=2, alpha=.7) +
xlab("Price (as % of original)") +
scale_x_continuous(labels = scales::percent_format(), n.breaks = 5) +
theme_minimal()
Incidence curves
While demand curves look similar, incidence curves reveal that
drastic price decreases lead to a smaller increase in people buying when
accounting for screening:
#demand curve compensatory
vd_demc_comp_inci =
testmarket %>%
ec_demcurve_inci(focal_alternatives,
seq(0.25,1.5,,9),
vd_dem_vdm,
out_pizza_cal,
eps_not)
#demand curve conjunctive screening
vd_demc_screen_inci =
testmarket %>%
ec_demcurve_inci(focal_alternatives,
seq(0.25,1.5,,9),
vd_dem_vdm_screen,
out_pizza_screening_cal,
eps_not)
#combine demand curves from both models
vd_outputs_inci=rbind(
vd_demc_comp_inci %>% do.call('rbind',.) %>% bind_cols(model='comp') %>% bind_cols(demand='volumetric'),
vd_demc_screen_inci %>% do.call('rbind',.) %>% bind_cols(model='screenpr') %>% bind_cols(demand='volumetric')) %>%
mutate(`E(purchase)`=`E(demand)`)
#plot
vd_outputs_inci%>%
filter(brand=="Priv") %>%
ggplot(aes(x=scenario, y=`E(purchase)`, color=model)) + geom_line(size=2, alpha=.7) +
xlab("Price (as % of original)") +
scale_x_continuous(labels = scales::percent_format(), n.breaks = 5) +
theme_minimal()
Incidence curves for all brands
Generating the same plot for all major brands, we can see differences
in incidence between the different models.
allbrands=names(table(testmarket$brand))
demcout=list()
for(kk in 1:6){
demcout[[kk]]=
ec_demcurve_inci(testmarket,
testmarket$brand==allbrands[kk],
c(.5,.6,.7,.8,.9,1),
vd_dem_vdm_screen,
out_pizza_screening_cal,
ec_gen_err_ev1(testmarket,out_pizza_screening_cal,seed = 34543563))
}
demcout_vd=list()
for(kk in 1:6){
demcout_vd[[kk]]=
ec_demcurve_inci(testmarket,
testmarket$brand==allbrands[kk],
c(.5,.6,.7,.8,.9,1),
vd_dem_vdm,
out_pizza_cal,
ec_gen_err_ev1(testmarket,out_pizza_cal,seed = 34543563))
}
names(demcout)=allbrands
names(demcout_vd)=allbrands
demcurves_pizza=
demcout %>% purrr::map_dfr(. %>% do.call('rbind',.),.id='focal') %>%
filter(focal==brand) %>% bind_cols(model='with screening') %>%
bind_rows(
demcout_vd %>% purrr::map_dfr(. %>% do.call('rbind',.),.id='focal') %>%
filter(focal==brand) %>% bind_cols(model='without screening'))
demcurves_pizza %>%
mutate(model=factor(model,levels=c('without screening',"with screening"))) %>%
ggplot(aes(x=scenario,y=`E(demand)`, linetype=model)) +
geom_line() +
facet_wrap(~focal, scales = 'free_y') +
xlab("Price (as % of original)") + ylab('Number of people buying') +
scale_x_continuous(labels = scales::percent_format(), n.breaks = 3) + theme_minimal()
