dpalma wrote: ↑09 Aug 2021, 18:17
Hi David,
No, you cannot use apollo_deltaMethod to calculate WTP when you have random coefficients. apollo_deltaMethod works only for fixed (as in non-random parameters). Instead, you have a couple of alternatives to calculate WTP when using random coefficients.
- Do simulation: That is what you did in your code. For this to be valid, you need to make sure that the ratio between the coefficients has finite moments (i.e. finite mean and s.d.). This is true if your numerator (b_attribute) follows a normal distribution and your denominator (b_price) is a fixed parameters; or both are log-normals. But it is not defined if both numerator and denominator follow a normal distribution; and in several other cases. I recommend you looking at Daily, Hess and Train (2012) Assuring finite moments for willingness to pay in random coefficient models
Best
David
I have a related but different question.
I have estimated mixed multinomial logistic regressions on a stated preference DCE with inter-individual draws for all attributes including time, the attribute we planned to use for willingness to spend time. Time was estimated using uniform draws, and the other attributes were estimated using normal draws. Dummy-coded models clearly show a non-linear specification for time use, so we have specified the time attribute as entering the utility function in quadratic form based on the following two publications:
van der Pol M, Ryan M. Specification of the Utility Function in Discrete Choice Experiments. Value Heal. 2014;17(2):297-301.
Torres C, Hanley N, Riera A. How wrong can you be? Implications of incorrect utility function specification for welfare measurement in choice experiments. J Environ Econ Manage. 2011;62:111-121.
I have tried conducting simulation for willingness to spend time as recommended in your post above, but it produces NaN for all simulated attribute values. I cannot tell if 1) I've made an error, 2) it is invalid to estimate WTP using an attribute whose estimated mean is centered around zero, 3) something else about quadratic specification needs to be revised, 4) something else I don't know. Code and output provided below.
I am also wondering if it is even possible to directly estimate a quadratic WTP attribute in WTP space in Apollo because of the plus/minus functionality in finding quadratic solutions. Do you have any examples of this or can you advise?
Code: Select all
#### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
>
> #### Check cores
> parallel::detectCores()
[1] 4
>
> #### Set core controls
> apollo_control = list(
+ modelName ="MXL_3",
+ modelDescr ="Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
+ reporting time continous and quadratic specification,
+ main effects only, all attributes + ASCs random, MLHS draws",
+ indivID ="RID",
+ mixing = TRUE,
+ nCores = 2
+ #analyticGrad = TRUE # Needs to be explicit for models with inter-intra draws (requires extra RAM).
+ )
>
> # inter-individual draws are used for variation across individuals
> # intra-individual draws for variation across observations for the same individual
>
> #### Set database
> database = df1
>
> ###################################################################
> #### DEFINE MODEL PARAMETERS ####
> ###################################################################
>
> # Utility U_ij = sum[1 to K](beta_jk * X_jk) + epsilon_ij
> #
> # where U_ij is the ultity of profile j for individual i
> # X_jk where k = 1...K are attributes in the dce
>
> ### Vector of parameters, including any that are kept fixed in estimation
> apollo_beta = c(mu_asc_a = -3,
+ sigma_asc_a = -0.01,
+ mu_asc_b = -3,
+ sigma_asc_b = -0.01,
+ b_violate_low = 0, # Violate_low is kept fixed
+ mu_b_investigate = -3,
+ sigma_b_investigate_inter = -0.01,
+ mu_b_time = -3,
+ sigma_b_time_inter = -0.01,
+ mu_b_time_sq = -3,
+ sigma_b_time_sq_inter = -0.01,
+ mu_b_anonymous = -3,
+ sigma_b_anonymous_inter = -0.01,
+ mu_b_amnesty = -3,
+ sigma_b_amnesty_inter = -0.01,
+ mu_b_violate_med = -3,
+ sigma_b_violate_med_inter = -0.01,
+ mu_b_violate_high = -3,
+ sigma_b_violate_high_inter = -0.01)
>
> ### 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("b_violate_low")
>
> #####################################################################
> #### DEFINE RANDOM COMPONENTS ####
> #####################################################################
>
> ### Set parameters for generating draws
> apollo_draws = list(
+ interDrawsType = "MLHS",
+ interNDraws = 500,
+ interUnifDraws = c("draws_time_inter", "draws_time_sq_inter"),
+ interNormDraws = c("draws_asc_a", "draws_asc_b",
+ "draws_anonymous_inter", "draws_amnesty_inter",
+ "draws_investigate_inter","draws_violate_med_inter",
+ "draws_violate_high_inter"),
+ intraDrawsType = "",
+ intraNDraws = 0,
+ intraNormDraws = c(),
+ intraUnifDraws = c()
+ )
>
> ### Create random parameters
> apollo_randCoeff = function(apollo_beta, apollo_inputs){
+ randcoeff = list()
+
+ randcoeff[["asc_a"]] = mu_asc_a + sigma_asc_a * draws_asc_a
+ randcoeff[["asc_b"]] = mu_asc_b + sigma_asc_b * draws_asc_b
+ randcoeff[["b_time"]] = mu_b_time + sigma_b_time_inter * draws_time_inter
+ randcoeff[["b_time_sq"]] = mu_b_time_sq + sigma_b_time_sq_inter * draws_time_sq_inter
+ randcoeff[["b_anonymous"]] = mu_b_anonymous + sigma_b_anonymous_inter * draws_anonymous_inter
+ randcoeff[["b_amnesty"]] = mu_b_amnesty + sigma_b_amnesty_inter * draws_amnesty_inter
+ randcoeff[["b_investigate"]] = mu_b_investigate + sigma_b_investigate_inter * draws_investigate_inter
+ randcoeff[["b_violate_med"]] = mu_b_violate_med + sigma_b_violate_med_inter * draws_violate_med_inter
+ randcoeff[["b_violate_high"]] = mu_b_violate_high + sigma_b_violate_high_inter * draws_violate_high_inter
+ # randcoeff[[""]] = mu_b_ + sigma_b__inter * draws__inter
+
+ # randcoeff[["b_investigate"]] = mu_b_investigate
+ # + sigma_b_investigate_inter * draws_investigate_inter
+ # + sigma_b_investigate_intra * draws_investigate_intra
+
+ return(randcoeff)
+ }
>
> #####################################################################
> #### GROUP AND VALIDATE INPUTS ####
> #####################################################################
>
> apollo_inputs = apollo_validateInputs()
Several observations per individual detected based on the value of RID. Setting panelData in apollo_control
set to TRUE.
All checks on apollo_control completed.
All checks on database completed.
Generating inter-individual draws ......... Done
>
> #####################################################################
> #### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
> #####################################################################
>
> apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+
+ ### Function initialisation: do not change the following three commands
+ ### 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']] = asc_a + b_anonymous * anonymous_a + b_time * time_a + b_time_sq * time_a * time_a + b_investigate * investigate_a + b_violate_low * (violate_a == 0) + b_violate_med * (violate_a == 1) + b_violate_high * (violate_a == 2) + b_amnesty * amnesty_a
+ V[['b']] = asc_b + b_anonymous * anonymous_b + b_time * time_b + b_time_sq * time_b * time_b + b_investigate * investigate_b + b_violate_low * (violate_b == 0) + b_violate_med * (violate_b == 1) + b_violate_high * (violate_b == 2) + b_amnesty * amnesty_b
+ V[['c']] = 0
+
+ ### Define settings for MNL model component
+ mnl_settings = list(
+ alternatives = c(a=1, b=2, c=3),
+ avail = list(a=1, b=1, c=1),
+ 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 ####
> #####################################################################
>
> # ### Optional speedTest
> # speedTest_settings=list(
> # nDrawsTry = c(100, 200, 500),
> # nCoresTry = 1:10,
> # nRep = 10
> # )
> #
> # apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
>
> MXL_3 = apollo_estimate(apollo_beta, apollo_fixed,
+ apollo_probabilities, apollo_inputs,
+ estimate_settings=list(hessianRoutine="maxLik"))
Testing likelihood function...
Overview of choices for MNL model component :
a b c
Times available 8032.00 8032.00 8032.00
Times chosen 3412.00 3314.00 1306.00
Percentage chosen overall 42.48 41.26 16.26
Percentage chosen when available 42.48 41.26 16.26
Pre-processing likelihood function...
Preparing workers for multithreading...
Testing influence of parameters..................
Starting main estimation
Initial function value: -294565.4
Initial gradient value:
mu_asc_a sigma_asc_a mu_asc_b sigma_asc_b
3411.98161 -150.96857 3313.98355 -143.72177
mu_b_investigate sigma_b_investigate_inter mu_b_time sigma_b_time_inter
3876.99834 -173.52127 17608.96516 8514.50448
mu_b_time_sq sigma_b_time_sq_inter mu_b_anonymous sigma_b_anonymous_inter
56644.96516 25272.75713 3962.99842 -178.90226
mu_b_amnesty sigma_b_amnesty_inter mu_b_violate_med sigma_b_violate_med_inter
4031.99838 -186.53884 1717.00000 -41.88887
mu_b_violate_high sigma_b_violate_high_inter
3502.99699 -146.35327
initial value 294565.439042
iter 2 value 45985.610138
iter 3 value 9908.723343
iter 4 value 9124.078980
iter 5 value 7984.555492
iter 6 value 7856.865565
iter 7 value 7736.338888
iter 8 value 7667.604144
iter 9 value 7650.549835
iter 10 value 7637.570218
iter 11 value 7566.365249
iter 12 value 7536.284985
iter 13 value 7458.709208
iter 14 value 7395.036999
iter 15 value 7271.401961
iter 16 value 7229.744813
iter 17 value 7201.647174
iter 18 value 7107.481110
iter 19 value 7077.874576
iter 20 value 6887.624995
iter 21 value 6877.596519
iter 22 value 6753.768691
iter 23 value 6724.593713
iter 24 value 6646.528796
iter 25 value 6541.358854
iter 26 value 6516.179882
iter 27 value 6501.455471
iter 28 value 6472.533695
iter 29 value 6466.063839
iter 30 value 6462.003358
iter 31 value 6458.907129
iter 32 value 6457.497947
iter 33 value 6456.925265
iter 34 value 6456.645028
iter 35 value 6456.498729
iter 36 value 6456.457763
iter 37 value 6456.398473
iter 38 value 6456.380287
iter 39 value 6456.379749
iter 40 value 6456.373976
iter 41 value 6456.370618
iter 42 value 6456.369985
iter 43 value 6456.369099
iter 44 value 6456.368444
iter 45 value 6456.367087
iter 46 value 6456.362775
iter 47 value 6456.361955
iter 48 value 6456.361509
iter 49 value 6456.361195
iter 50 value 6456.356222
iter 51 value 6456.353742
iter 52 value 6456.353581
iter 53 value 6456.353452
iter 53 value 6456.353389
iter 53 value 6456.353366
final value 6456.353366
converged
Additional convergence test using scaled estimation. Parameters will be scaled by their current estimates
and additional iterations will be performed.
initial value 6456.353366
iter 2 value 6456.352860
iter 3 value 6456.352588
iter 3 value 6456.352566
iter 3 value 6456.352566
final value 6456.352566
converged
Estimated parameters:
Estimate
mu_asc_a -1.50323
sigma_asc_a -0.95509
mu_asc_b -1.59154
sigma_asc_b -0.97108
b_violate_low 0.00000
mu_b_investigate 1.05741
sigma_b_investigate_inter -1.16282
mu_b_time -0.17974
sigma_b_time_inter 1.45248
mu_b_time_sq -0.12585
sigma_b_time_sq_inter -0.01072
mu_b_anonymous 1.16539
sigma_b_anonymous_inter -1.19584
mu_b_amnesty 1.29390
sigma_b_amnesty_inter -1.27754
mu_b_violate_med 1.44376
sigma_b_violate_med_inter -0.93255
mu_b_violate_high 1.74279
sigma_b_violate_high_inter -1.16227
Computing covariance matrix using numerical methods (maxLik). This may take a while, no progress bar
displayed.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a saddle point!
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
Calculating LL of each model component...
Warning message:
In sqrt(diag(varcov)) : NaNs produced
>
> # ----------------------------------------------------------------- #
> #---- FORMATTED OUTPUT (TO SCREEN) ----
> # ----------------------------------------------------------------- #
>
> apollo_modelOutput(MXL_3)
Model run using Apollo for R, version 0.2.6 on Darwin by aliceellyson
www.ApolloChoiceModelling.com
Model name : MXL_3
Model description : Mixed logit (MXL) on DCE SP data, dummy-coded attributes,
reporting time continous and quadratic specification,
main effects only, all attributes + ASCs random, MLHS draws
Model run at : 2022-01-12 15:29:50
Estimation method : bfgs
Model diagnosis : successful convergence
Number of individuals : 1004
Number of rows in database : 8032
Number of modelled outcomes : 8032
Number of cores used : 2
Number of inter-individual draws : 500 (MLHS)
LL(start) : -294565.4
LL(0) : -8824.054
LL(C) : -8227.245
LL(final) : -6456.353
Rho-square (0) : 0.2683
Adj.Rho-square (0) : 0.2663
Rho-square (C) : 0.2152
Adj.Rho-square (C) : 0.2131
AIC : 12948.71
BIC : 13074.55
Estimated parameters : 18
Time taken (hh:mm:ss) : 00:31:21.28
pre-estimation : 00:03:19.52
estimation : 00:07:7.04
post-estimation : 00:20:54.71
Iterations : 59
Min abs eigenvalue of Hessian : 8455.181
Some eigenvalues of Hessian are positive, indicating potential problems!
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
mu_asc_a -1.50323 0.26855 -5.5975 0.297467 -5.0534
sigma_asc_a -0.95509 0.07988 -11.9573 0.097385 -9.8074
mu_asc_b -1.59154 0.27041 -5.8857 0.300556 -5.2953
sigma_asc_b -0.97108 0.08062 -12.0449 0.093964 -10.3346
b_violate_low 0.00000 NA NA NA NA
mu_b_investigate 1.05741 0.06812 15.5221 0.070229 15.0566
sigma_b_investigate_inter -1.16282 0.08624 -13.4837 0.093665 -12.4147
mu_b_time -0.17974 0.26788 -0.6710 0.285514 -0.6295
sigma_b_time_inter 1.45248 0.08103 17.9251 0.084762 17.1361
mu_b_time_sq -0.12585 0.04955 -2.5395 0.054294 -2.3178
sigma_b_time_sq_inter -0.01072 NaN NaN 0.008842 -1.2122
mu_b_anonymous 1.16539 0.06967 16.7264 0.075518 15.4320
sigma_b_anonymous_inter -1.19584 0.08475 -14.1095 0.093386 -12.8053
mu_b_amnesty 1.29390 0.07235 17.8831 0.078648 16.4518
sigma_b_amnesty_inter -1.27754 0.08509 -15.0131 0.091237 -14.0024
mu_b_violate_med 1.44376 0.10844 13.3137 0.120448 11.9866
sigma_b_violate_med_inter -0.93255 0.13866 -6.7253 0.161731 -5.7660
mu_b_violate_high 1.74279 0.08884 19.6165 0.096355 18.0872
sigma_b_violate_high_inter -1.16227 0.08928 -13.0179 0.095641 -12.1524
> # Simulate values
> N = 10000
> # Attributes
> b_attr_anon = rnorm(N, mean=MXL_3$estimate["mu_b_anonymous"], sd=MXL_3$estimate["sigma_b_anonymous_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_anonymous"], sd = MXL_3$estimate["sigma_b_anonymous_inter"]) :
NAs produced
> b_attr_immu = rnorm(N, mean=MXL_3$estimate["mu_b_amnesty"], sd=MXL_3$estimate["sigma_b_amnesty_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_amnesty"], sd = MXL_3$estimate["sigma_b_amnesty_inter"]) :
NAs produced
> b_attr_inve = rnorm(N, mean=MXL_3$estimate["mu_b_investigate"], sd=MXL_3$estimate["sigma_b_investigate_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_investigate"], sd = MXL_3$estimate["sigma_b_investigate_inter"]) :
NAs produced
> b_attr_MED = rnorm(N, mean=MXL_3$estimate["mu_b_violate_med"], sd=MXL_3$estimate["sigma_b_violate_med_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_violate_med"], sd = MXL_3$estimate["sigma_b_violate_med_inter"]) :
NAs produced
> b_attr_HIGH = rnorm(N, mean=MXL_3$estimate["mu_b_violate_high"], sd=MXL_3$estimate["sigma_b_violate_high_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_violate_high"], sd = MXL_3$estimate["sigma_b_violate_high_inter"]) :
NAs produced