WTP for random coefficients in mixed logit

[rory_m42]

Hello,

I have a question about estimating WTP in the preference space.

I am estimating a mixed logit model. I have assumed all variables to be normally distributed, except for the cost variable which I have assumed as a negative lognormal, and the reference variables which I have fixed. The code for model estimation is as follows:

Code: Select all

choiceAnalysis_settings <- list(
  alternatives = c(alt1=1, alt2=2, alt3=3),
  avail        = list(alt1=database$alt1.a1, alt2=database$alt2.a1, alt3=database$alt3.a1),
  choiceVar    = database$Preference)



# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #


### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(mu_community_engagement      = 0,
              mu_basic_community_services     = 0,
              mu_community_holidays      = 0,
              mu_music_community     = 0,
              mu_no_maintenance     = 0,
              mu_essential_repair     = 0,
              mu_Comprehensive_repair     = 0,
              mu_open_6days    = 0,
              mu_weekly_service     = 0,
              mu_daily_reflection       = 0,
              mu_minimal_tech  = 0,
              mu_basic_tech       = 0,
              mu_high_tech      = 0,
              mu_price         = -3,
              
              
            
              sigma_basic_community_services     = 0,
              sigma_community_holidays      = 0,
              sigma_music_community     = 0,
       
              sigma_essential_repair     = 0,
              sigma_Comprehensive_repair     = 0,
       
              sigma_weekly_service     = 0,
              sigma_daily_reflection       = 0,
      
              sigma_basic_tech       = 0,
              sigma_high_tech      = 0,
              sigma_price         = 0
              
              
              
)


### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c("mu_community_engagement", "mu_no_maintenance", "mu_open_6days", "mu_minimal_tech")

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "mlhs",
  interNDraws    = 200,
  interUnifDraws = c(),
  interNormDraws = c("draws_basic_community","draws_community_holidays","draws_music_community", "draws_essential_repair", "draws_comprehensive_repair", "draws_weekly_service", "draws_daily_reflection",  "draws_basic_tech", "draws_high_tech", "draws_price"), 
  
  intraDrawsType = "pmc",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)


apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  

  randcoeff[["basic_community_services"]] = ( mu_basic_community_services + sigma_basic_community_services * draws_basic_community)
  randcoeff[["community_holidays"]] = (mu_community_holidays + sigma_community_holidays * draws_community_holidays)
  randcoeff[["music_community"]] = ( mu_music_community + sigma_music_community * draws_music_community)

  randcoeff[["essential_repair"]] = ( mu_essential_repair + sigma_essential_repair * draws_essential_repair)
  randcoeff[["Comprehensive_repair"]] = ( mu_Comprehensive_repair + sigma_Comprehensive_repair * draws_comprehensive_repair)

  randcoeff[["weekly_service"]] = ( mu_weekly_service + sigma_weekly_service * draws_weekly_service)
  randcoeff[["daily_reflection"]] = ( mu_daily_reflection + sigma_daily_reflection* draws_daily_reflection)

  randcoeff[["basic_tech"]] = ( mu_basic_tech + sigma_basic_tech * draws_basic_tech)
  randcoeff[["high_tech"]] = ( mu_high_tech + sigma_high_tech * draws_high_tech)
  
  randcoeff[["price"]] = -exp(mu_price + sigma_price * draws_price)
  
  
  return(randcoeff)
  
}




apollo_inputs = apollo_validateInputs()

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))
  
  
  
  ### Create list of probabilities P
  P = list()
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[['A']]  = ( 
    mu_community_engagement     * ( alt1.a1 == 1)    
    + basic_community_services  * ( alt1.a1 == 2)
    + community_holidays       * ( alt1.a1 == 3)        
    + music_community       * ( alt1.a1 == 4)       
    + mu_no_maintenance   * ( alt1.b1 == 1)   
    + essential_repair      * ( alt1.b1 == 2)      
    + Comprehensive_repair      * ( alt1.b1 == 3)
    + mu_open_6days        * ( alt1.c1 == 1)
    + weekly_service       * ( alt1.c1 == 2)
    + daily_reflection    *  (alt1.c1 == 3)
    + mu_minimal_tech        * ( alt1.d1 == 1)
    + basic_tech       * ( alt1.d1 == 2)
    + high_tech    *  (alt1.d1 == 3)
    + price                * ( ( alt1.e1 == 2) * 2
                               +( alt1.e1 == 4) * 4
                               +( alt1.e1 == 6) * 6
                               +( alt1.e1 == 8) * 8
                               +( alt1.e1 ==10)  * 10
                               +( alt1.e1 ==12)  * 12 
                               +( alt1.e1 ==14)  * 14 
                               +( alt1.e1 ==16)  * 16 ))
  
  
  
  
  V[['B']]  = ( 
    mu_community_engagement     * ( alt2.a1 == 1)    
    + basic_community_services  * ( alt2.a1 == 2)
    + community_holidays       * ( alt2.a1 == 3)        
    + music_community       * ( alt2.a1 == 4)       
    + mu_no_maintenance   * ( alt2.b1 == 1)   
    + essential_repair      * ( alt2.b1 == 2)      
    + Comprehensive_repair      * ( alt2.b1 == 3)
    + mu_open_6days        * ( alt2.c1 == 1)
    + weekly_service       * ( alt2.c1 == 2)
    + daily_reflection    *  (alt2.c1 == 3)
    + mu_minimal_tech        * ( alt2.d1 == 1)
    + basic_tech       * ( alt2.d1 == 2)
    + high_tech    *  (alt2.d1 == 3)
    + price                * ( ( alt2.e1 == 2) * 2
                               +( alt2.e1 == 4) * 4
                               +( alt2.e1 == 6) * 6
                               +( alt2.e1 == 8) * 8
                               +( alt2.e1 ==10)  * 10
                               +( alt2.e1 ==12)  * 12 
                               +( alt2.e1 ==14)  * 14 
                               +( alt2.e1 ==16)  * 16 ))
  
  V[['C']]  = (  
    mu_community_engagement     *   ( alt3.a1 == 1)    
    + mu_no_maintenance   * ( alt3.b1 == 1)      
    + mu_open_6days         * ( alt3.c1 == 1)
    + mu_minimal_tech     *  (alt3.d1  ==1)
    + price          *  ( alt3.e1_status == 0) *0)
  
  
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(A=1,B=2, C=3),
    choiceVar     = Preference,
    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 = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings=list(maxIterations=1000))

apollo_modelOutput(model)

Based on previous posts in this forum and some online literature, I think I have a number of options for calculating WTP. I am hoping someone can help me decide which is correct/the best one to use.

1) Simply take the ratio of the marginal utilities such that WTP = mu_attribute/mu_price

2) If the moments of the WTP distribution are finite, then one can find the mean of the WTP by sampling from each of the distributions of the coefficients in the WTP calculation. I think for my case this is as follows:

Code: Select all

N=10000

x = rnorm(N, mean= model$estimate['mu_essential_repair'], sd = abs(model$estimate['sigma_essential_repair']))
y = -exp(rlnorm(N, mean= model$estimate['mu_price'], sd = abs(model$estimate['sigma_price'])))
wtp = x/y
mean(wtp)

3) I have seen use of the unconditionals function in the apollo package, of which I am unsure how to proceed.

Any help is really appreciated!

[stephanehess]

Hi

option 1 would incorrect, as that would be a ratio of means, while you want a mean of ratios.

option 2 and 3 are the same, but option 3 is much easier. You can use apollo_unconditionals to get the distributed parameters and then take the ratio of them.

E.g.

unconditionals = apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
wtp_basic_tech = as.vector(unconditionals[["basic_tech"]])/as.vector(unconditionals[["price"]])

[rory_m42]

Hi Stephane,

Thank you for your reply. If its okay, I have a few follow up questions.

We have used the code you provided to get WTP estimates in preference space, but found untennably high mean and sd values. This encouraged us to move to WTP space, where estimates became much more reasonable. In WTP space, we have notied that the standard errors are very small, and the t-values very very high. Is this cause for concern or a common bi-product of working in WTP space? We have seen low standard errors appear in a few papers e.g. https://link.springer.com/article/10.10 ... 21-09406-7.

Furthermore, we have a sample that asks identical questions but for different sites. To analyse the impact of the site on WTP, we want to interact a binary site variable with the other coefficients in the model. This would mean interacting with every coefficient in the model, but most notably with the price coefficient. This complicates the WTP calculation in WTP space. Is it possible to interact the price variable, or is this something that WTP space does not allow?

Thank you!

[dpalma]

Hi,

About the first question, having very small s.e., or equivalently, very large t-tests, is not a problem on itself, but it could be indicative of another issue, notably poor convergence of the model, leading to an ill-defined Hessian and covariance matrix. I would double check that your model did not converge to a saddle point. To do this, check the "Optimisation diagnosis" section in the output. There is an example below. You would expect the hessian to be "Negative definitive", and the maximum eigenvalue to be negative. Again, the output below is just an example, your model does not need to have similar values to these.

Code: Select all

Optimisation diagnosis                      : Maximum found
     hessian properties                     : Negative definite
     maximum eigenvalue                     : -0.79691
     reciprocal of condition number         : 3.93679e-08

Concerning your second question, if you interact all coefficients of your model with a dummy variable, for example by doing something like:

Code: Select all

V[["alt1"]] = (b1_a*dummy_a + b1_b*(1 - dummy_a))*x1 + ...

Then that is the same as estimating two separate models: one only for people in group a (dummy_a=1) and another for people not in that group (dummy_a = 0). So the way you calculate the WTP does not change, you just use b1_a/bCost_a for group "a", and b1_b/bCost_b for group "b".

This can be applied to the case of WTP space as well. Your utility function would look something like the following.

Code: Select all

V[["alt1"]] = (bCost_a*dummy_a + bCost_b*(1 - dummy_a)) * ( (b1_a*dummy_a + b1_b*(1 - dummy_a))*x1 + ... )

I hope this clarifies the situation.

Best wishes
David

[rory_m42]

Hi David,

Thank you for your response. In regard to the second question, the attempted specification is slightly different to the one you provided. We have one sample of respondents who are asked identical questions. However, the questions for each respondent refer to a different site and we want to test for the impact/significance of the site on the respondents WTP. In this case, suppose there are two sites, A & B, and we want to see the impact of being asked about site A. The specification we have attempted is the following (in preference space):

Code: Select all

V[["alt1"]] = bCost*price + b_1*x + bCost*price*dummy_A + b_1*x*dummy_A

The confusion is then whether this specification even makes sense, and then how it would be converted to WTP space given the interaction with price.

Thank you for your help.

[dpalma]

Hi,

If you want to measure differences in the effect of attributes due to the location, you will need additional coefficients for each attribute. So your specification in preference space would look as follows:

Code: Select all

V[["alt1"]] = bCost*price + b_1*x_1 + bCost_A*price*dummy_A + b_1A*x_1*dummy_A

Which you can rewrite as follows:

Code: Select all

V[["alt1"]] = (bCost + bCost_A*dummy_A)*price + (b_1 + b_1A*dummy_A)*x_1

Note that we have two coefficients for each attribute: bCost and bCost_A for price, and b_1 and b_1A for attribute x_1. Assuming that dummy_A takes value 1 when the location is “A”, and zero when location is “B”, then the effect of price for location “A” will be (bCost + bCost_A), and the effect for location “B” will be (bCost). So if bCost_A is not significantly different from zero, the effect of price does not change due to the location, i.e. the effect is the same in both locations.

Based on the second formulation, you can move to willingness to pay space as follows:

Code: Select all

V[["alt1"]] = (bCost + bCost_A*dummy_A)*( (w1 + w1_A)*x_1 + price )

In this formulation, the WTP for attribute x_1 is (w1 + w1_A) in location “A”, and (w1) in location “B”.

I hope this helps clarifies the confusion.

Best wishes,
David