Page 1 of 1
WTP Space Mixed Logit with piecewise price parameter
Posted: 30 Oct 2021, 17:16
by sac111
Dear Stephane,
Thank you for creating this package and the forum it has been really helpful so far.
I am trying to estimate a mixed logit model in the WTP space for a discrete choice experiment.
The price attribute is an expected change in bill in percentage points(this can be either positive negative or zero).
I would like to be able to capture the asymmetry between WTP and WTA so have used a piecewise linear transformation to split the bill attribute into a fee (bill change >=0) or discount (bill change <=0) following the method here:
S. Hess, J. M. Rose, and D. A. Hensher, “Asymmetric preference formation in willingness to pay estimates in discrete choice models,” Transp. Res. Part E Logist. Transp. Rev., vol. 44, no. 5, pp. 847–863, Sep. 2008, doi: 10.1016/j.tre.2007.06.002.
https://linkinghub.elsevier.com/retriev ... 4507000786
This works fine for a mnl specification and there is a difference in the fee and discount parameters (mnl_models.png).
For the mixed logit specification I am using log-normal distributions for the fee and discount, and normal distributions for the non-monetary attributes. In the preference space the mixed logit results look fine. However the main aim of the study is to look at WTP/WTA and associated confidence intervals so would like to estimate the model directly in the WTP space.
I have tried two versions of specifying the utility functions in the WTP:
WTP_V1 - including the other price attribute within the WTP space (within the brackets). I end up with NaNs for the standard errors of some of the variables (WTP_V1_FEE.txt & WTP_V1_DISC.txt).
WTP_V2 - keeping the other price attribute outside the WTP space (outside the brackets). Works for the WTP estimation with plausible results but does not converge for the WTA estimation (WTP_V2_FEE.txt & WTP_V2_DISC.txt).
My question is first whether specification in the WTP space is the best way to go for this is analysis?
Secondly does WTP_V2 make sense?
Lastly are the convergence issues relating to the WTA estimates a problem with starting values?
I have included the model specifications and outputs for the preference space and the two WTP versions.
Thank you,
Saurab
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 01 Nov 2021, 13:01
by stephanehess
Hi
personally, I would avoid WTP space in your case. As soon as you have multiple cost parameters, you would need to make a decision in relation to which cost component you express wtp, and then you'd still need post-estimation calculations to transform for other valuations (WTP vs WTA).
Also be mindful that if you use Normals for the non-price terms, then WTP space and preference space would be different anyway, as the WTP in the preference space model would be a ratio of normal over lognormal, which is different from a normal, which is what you'd have in WTP space.
If you're working in preference space, the model will be numerically easier too, so you can then also explore using semi-non-parametric transforms of the distribution
Stephane
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 02 Nov 2021, 10:57
by sac111
Dear Stephane,
Thank you for the quick response!
I think I will stick to the preference space then and will look into the semi-non-parametric distributions.
I am assuming the best way to compute confidence intervals for the mean WTP/WTA would then be to use simulations (Krinsky-Robb method)?
Thank you,
Saurab
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 02 Nov 2021, 13:21
by stephanehess
Saurab
re the CIs, I would sample from the asymptotic variance covariance matrix for all parameters, including the covariances, and then each time simulate across the random distribution, and then calculate a confidence interval for the means
Stephane
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 04 Nov 2021, 17:01
by sac111
Dear Stephane,
Thank you.
I m trying to implement a zero censored normal for some the attributes as a fairly large minority of respondents would be expected to be indifferent to them. I used the pmax() function in the randcoeff specification however I get the following warning when running apollo_estimate() (full script below).
Code: Select all
##############################################################################################################################################
This model could potentially be estimated faster using analytical gradients, yet some issue is preventing it from using them.
You might want to ask for help in the Apollo forum (http://www.apollochoicemodelling.com/forum) on how to solve this issue.
If you do, please post your code and data (if not confidential).
##############################################################################################################################################
The estimation is indeed very slow. Is there a better way to specify the censored normal (using an analytical expression) in apollo?
Thank you,
Saurab
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "INT_MXL_PREF_CN",
modelDescr ="MXL model with interaction, pref, ln-price and cn anon",
indivID ="ID",
panelData = TRUE,
mixing = TRUE,
nCores = 8,
seed = 423543
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = readRDS(here::here("data/formatted", "wide_data.rds"))
database <- subset(database, filt.trade == 1 & filt.manip2==1)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
apollo_beta=c(
mu.fee = -3,
mu.disc = -3,
mu.anon = 0,
mu.hh = 0,
mu.daily = 0,
mu.anon.hh = 0,
mu.anon.daily = 0,
sg.fee = 0,
sg.disc = 0,
sg.anon = 0,
sg.hh = 0,
sg.daily = 0,
sg.anon.hh = 0,
sg.anon.daily = 0
)
apollo_fixed = c() # Reference real-time
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "sobolOwenFaureTezuka",
interNDraws = 1000,
interUnifDraws = c(),
interNormDraws = c(
"draws.fee",
"draws.disc",
"draws.anon",
"draws.hh",
"draws.daily",
"draws.anon.hh",
"draws.anon.daily"
),
intraDrawsType = "sobol",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["a.fee"]] = -exp(mu.fee + sg.fee*draws.fee)
randcoeff[["a.disc"]] = exp(mu.disc + sg.disc*draws.disc)
randcoeff[["b.anon"]] = pmax(mu.anon + sg.anon*draws.anon,0)
randcoeff[["b.hh"]] = pmax(mu.hh + sg.hh*draws.hh,0)
randcoeff[["b.daily"]] = pmax(mu.daily + sg.daily*draws.daily ,0)
randcoeff[["b.anon.hh"]] = mu.anon.hh + sg.anon.hh*draws.anon.hh
randcoeff[["b.anon.daily"]] = mu.anon.daily + sg.anon.daily*draws.anon.daily
return(randcoeff)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
### Attach inputs and detach after function exit
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()
V = list()
V[['1']] = a.fee*(FeeP_1 + a.disc*DiscP_1 + b.anon*(Anon_1==1) +
b.hh*(HH_1==1) + b.daily*(Daily_1==1) +
b.anon.hh*(HH_1==1 & Anon_1==1) + b.anon.daily*(Daily_1==1 & Anon_1==1))
V[['2']] = a.fee*(FeeP_2 + a.disc*DiscP_2 + b.anon*(Anon_2==1) +
b.hh*(HH_2==1) + b.daily*(Daily_2==1) +
b.anon.hh*(HH_2==1 & Anon_2==1) + b.anon.daily*(Daily_2==1 & Anon_2==1))
mnl_settings = list(
alternatives = c("1"=1, "2"=2),
avail = list("1"=av_1, "2"=av_2),
choiceVar = choice,
V = V
)
### Compute probabilities using MNL model
P[['model']] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
### Average across inter-individual draws
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
estimate_settings <- list(
# scaling = c(mu.fee = 20, mu.disc = 20)
)
model = apollo_estimate(
apollo_beta,
apollo_fixed,
apollo_probabilities,
apollo_inputs,
estimate_settings
)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
modelOutput_settings = list(printPVal=2)
apollo_modelOutput(model,modelOutput_settings)
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 05 Nov 2021, 15:28
by dpalma
Hi Saurab,
We are not sure what is preventing the model from using analytical gradients. It could be a number of things. If you can share the data (or part of it) with us, it would help us figuring out what the problem is. If you are ok with sharing the data, you can either post it in the forum, or write an email do D.Palma [at] leeds.ac.uk
Cheers
David
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 05 Nov 2021, 17:40
by sac111
Hi David,
I have emailed the dataset and full R script.
Thank you,
Saurab
Re: WTP Space Mixed Logit with piecewise price parameter
Posted: 12 Nov 2021, 10:53
by dpalma
Hi Saurab,
Thanks for sending the data and script.
I managed to run your code using analytical gradient, but with a few caveats:
- I used Apollo v0.2.6, which has just been released. To install it:
- If you are using MacOS or Linux: Simply type install.packages("apollo")
- If you are using Windows: Download the package from the webpage and the type install.packages("C:/pathToFile/apollo_0.2.6.zip", repos=NULL) to install it, where pathToFile needs to changed to the appropriate path.
- I had to change the draws type to MLHS. For some reason I am investigating, the sobolOwenFaureTezuka draws did not work.
- I changed the way you are generating the truncated normal. By simply taking the max(0, x), where x ~ N(0,1), you were generating a distribution with half its mass concentrated on 0. A truncated normal with only positive values, instead, should scale its mass to the right of zero so that it integrates to one. I changed your code to achieve that.
- You will probably get a warning about database being different in the global environment and inside apollo_inputs. We are investigating why this happens exactly, but it's nothing to worry about.
The updated code is below.
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "INT_MXL_PREF_CN",
modelDescr ="MXL model with interaction, pref, ln-price and cn anon",
indivID ="ID",
panelData = TRUE,
mixing = TRUE,
nCores = 5,
seed = 423543,
debug=TRUE
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = readRDS("wide_data.rds")
database <- subset(database, filt.trade == 1 & filt.manip2==1)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
apollo_beta=c(
mu.fee =-3,
mu.disc =-3,
mu.anon = 0,
mu.hh = 0,
mu.daily = 0,
mu.anon.hh = 0,
mu.anon.daily = 0,
sg.fee = 0,
sg.disc = 0,
sg.anon = 0,
sg.hh = 0,
sg.daily = 0,
sg.anon.hh = 0,
sg.anon.daily = 0
)
apollo_fixed = c() # Reference real-time
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs", #"sobolOwenFaureTezuka",
interNDraws = 1000,
interUnifDraws = c("draws.anon","draws.hh","draws.daily"),
interNormDraws = c("draws.fee","draws.disc","draws.anon.hh",
"draws.anon.daily")
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["a.fee"]] = -exp(mu.fee + sg.fee*draws.fee)
randcoeff[["a.disc"]] = exp(mu.disc + sg.disc*draws.disc)
randcoeff[["b.anon"]] = mu.anon + sg.anon *qnorm(0.5 + 0.5*draws.anon )
randcoeff[["b.hh"]] = mu.hh + sg.hh *qnorm(0.5 + 0.5*draws.hh )
randcoeff[["b.daily"]] = mu.daily + sg.daily*qnorm(0.5 + 0.5*draws.daily)
randcoeff[["b.anon.hh"]] = mu.anon.hh + sg.anon.hh*draws.anon.hh
randcoeff[["b.anon.daily"]] = mu.anon.daily + sg.anon.daily*draws.anon.daily
return(randcoeff)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
### Attach inputs and detach after function exit
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()
V = list()
V[['1']] = a.fee*FeeP_1 + a.disc*DiscP_1 + b.anon*(Anon_1==1) +
b.hh*(HH_1==1) + b.daily*(Daily_1==1) + b.anon.hh*(HH_1==1 & Anon_1==1) +
b.anon.daily*(Daily_1==1 & Anon_1==1)
V[['2']] = a.fee*FeeP_2 + a.disc*DiscP_2 + b.anon*(Anon_2==1) +
b.hh*(HH_2==1) + b.daily*(Daily_2==1) + b.anon.hh*(HH_2==1 & Anon_2==1) +
b.anon.daily*(Daily_2==1 & Anon_2==1)
mnl_settings = list(
alternatives = c("1"=1, "2"=2),
avail = list("1"=av_1, "2"=av_2),
choiceVar = choice,
V = V
)
### Compute probabilities using MNL model
P[['model']] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
### Average across inter-individual draws
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
#estimate_settings <- list(
# scaling = c(mu.fee = 20, mu.disc = 20)
#)
model = apollo_estimate(
apollo_beta,
apollo_fixed,
apollo_probabilities,
apollo_inputs#,
#estimate_settings
)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
modelOutput_settings = list(printPVal=2)
apollo_modelOutput(model,modelOutput_settings)
Cheers
David