MRS derivation from MMNL with Box-Cox transformation
Posted: 28 Oct 2024, 18:33
Dear prof. dr. Hess and dr. Palma,
In the MMNL model that I’ve estimated, I’ve specified six of the random parameters to follow a normal distribution (i.e., the policy attributes), two a positive lognormal distribution (i.e., the effectiveness attributes), and one a negative lognormal distribution (i.e., the cost attribute). Additionally, to adjust for non-linearities in the utilities derived from different levels of the effectiveness attributes, I’ve applied a Box-Cox transformation to both of these attributes. Finally, to mitigate exploding implicit prices and ensure finite moments for the WTP distribution, I’ve applied a “mu-shift” to the cost attribute following Crastes dit Sourd (2024) . See the code and output below.
After estimation, I’d now like to derive the MRS between each of the other attributes and the cost attribute. I’ve got two questions on this:
1) I’ve already tried to derive the MRS between each of the (normally distributed) policy attributes and the cost attribute. This seems to work for most attributes, but not for all. In the code provided below, this applies to ‘Price sunscreen 30% lower’. When I run the code for this attribute, the values of all the draws for ‘b_lowerprice_sunscreen’ are ‘NaN’ and, as a result, also those for ‘wtp_lowerprice_sunscreen’.
I suspect this has to do with this attribute having a negative coefficient for the sigma parameter (I’ve run three country-specific models, and I’m having the same issue for those (different) attributes that have a negative coefficient for the sigma parameter in the other countries). I believe that the sign of the coefficient of the sigma parameter does not matter for the interpretation, but perhaps it may matter for the derivation of the MRS in the approach used here? Do you think this could have caused the ‘NaN’ results? And, if so, how should I adapt my code accordingly?
2) For the derivation of the MRS between the two (positive lognormally distributed) effectiveness attributes and the cost attribute, I’m wondering how to adapt the code to reflect the Box-Cox transformation of these attributes. Preferably, I’d obtain MRS estimates for each of the four levels of each effectiveness attribute (as included in the DCE) and/or plot the non-linear WTP curve for different levels of the attribute. I’ve looked at the suggestions by Daly et al. (2012), but I’m not sure how to implement this in Apollo.
Thank you very much in advance for your reply!
Kind regards,
Sander
In the MMNL model that I’ve estimated, I’ve specified six of the random parameters to follow a normal distribution (i.e., the policy attributes), two a positive lognormal distribution (i.e., the effectiveness attributes), and one a negative lognormal distribution (i.e., the cost attribute). Additionally, to adjust for non-linearities in the utilities derived from different levels of the effectiveness attributes, I’ve applied a Box-Cox transformation to both of these attributes. Finally, to mitigate exploding implicit prices and ensure finite moments for the WTP distribution, I’ve applied a “mu-shift” to the cost attribute following Crastes dit Sourd (2024) . See the code and output below.
Code: Select all
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_left = 0,
asc_right_mu = 0,
asc_right_sig = 0,
b_informationcampaign_no = 0,
b_informationcampaign_yes_mu = 0.3248,
b_informationcampaign_yes_sig = 0.001,
b_solarbed_ban_no = 0,
b_solarbed_ban_yes_mu = -0.0193,
b_solarbed_ban_yes_sig = 0.001,
b_solarium_ban_no = 0,
b_solarium_ban_yes_mu = 0.0706,
b_solarium_ban_yes_sig = 0.001,
b_lowerprice_sunscreen_no = 0,
b_lowerprice_sunscreen_yes_mu =0.2836,
b_lowerprice_sunscreen_yes_sig = 0.001,
b_free_sunscreen_no = 0,
b_free_sunscreen_yes_mu = 0.0919,
b_free_sunscreen_yes_sig = 0.001,
b_free_app_no = 0,
b_free_app_yes_mu = 0.1863,
b_free_app_yes_sig = 0.001,
b_log_newcases_mu = log(0.0384),
b_log_newcases_sig = 0.001,
lambda_newcases = 2.3123,
b_log_deaths_mu = log(0.0502),
b_log_deaths_sig = 0.001,
lambda_deaths = 1.4772,
b_log_cost_mu = log(0.4149/2),
b_log_cost_sig = 0.001)
### 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("asc_left", "b_informationcampaign_no", "b_solarbed_ban_no",
"b_solarium_ban_no", "b_lowerprice_sunscreen_no", "b_free_sunscreen_no",
"b_free_app_no")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "sobol",
interNDraws = 5000,
interNormDraws = c("draws_asc", "draws_informationcampaign_yes","draws_solarbed_ban_yes","draws_solarium_ban_yes","draws_lowerprice_sunscreen_yes", "draws_free_sunscreen_yes", "draws_free_app_yes", "draws_newcases", "draws_deaths", "draws_cost")
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["asc_right"]] = asc_right_mu + asc_right_sig * draws_asc
randcoeff[["b_informationcampaign_yes"]] = b_informationcampaign_yes_mu + b_informationcampaign_yes_sig * draws_informationcampaign_yes
randcoeff[["b_solarbed_ban_yes"]] = b_solarbed_ban_yes_mu + b_solarbed_ban_yes_sig * draws_solarbed_ban_yes
randcoeff[["b_solarium_ban_yes"]] = b_solarium_ban_yes_mu + b_solarium_ban_yes_sig * draws_solarium_ban_yes
randcoeff[["b_lowerprice_sunscreen_yes"]] = b_lowerprice_sunscreen_yes_mu + b_lowerprice_sunscreen_yes_sig * draws_lowerprice_sunscreen_yes
randcoeff[["b_free_sunscreen_yes"]] = b_free_sunscreen_yes_mu + b_free_sunscreen_yes_sig * draws_free_sunscreen_yes
randcoeff[["b_free_app_yes"]] = b_free_app_yes_mu + b_free_app_yes_sig * draws_free_app_yes
randcoeff[["b_newcases"]] = exp(b_log_newcases_mu + b_log_newcases_sig * draws_newcases)
randcoeff[["b_deaths"]] = exp(b_log_deaths_mu + b_log_deaths_sig * draws_deaths)
randcoeff[["b_cost"]] = -exp(b_log_cost_mu)-exp(b_log_cost_mu + b_log_cost_sig * draws_cost)
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))
### Create list of probabilities P
P = list()
### Likelihood of choices
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[["PackageA"]] = asc_left + b_informationcampaign_no*(att1information_campaigns1==1) + b_informationcampaign_yes*(att1information_campaigns1==2) +
b_solarbed_ban_no*(att2prohibiting_solarbeds1==1) + b_solarbed_ban_yes*(att2prohibiting_solarbeds1==2) +
b_solarium_ban_no*(att3prohibiting_solarstudios1==1) + b_solarium_ban_yes*(att3prohibiting_solarstudios1==2) +
b_lowerprice_sunscreen_no*(att4pricereduction_sunscreen1==1) + b_lowerprice_sunscreen_yes*(att4pricereduction_sunscreen1==2) +
b_free_sunscreen_no*(att5free_sunscreen1==1) + b_free_sunscreen_yes*(att5free_sunscreen1==2) +
b_free_app_no*(att6app1==1) + b_free_app_yes*(att6app1==2) +
b_newcases*(att7new_cases1^lambda_newcases-1)/lambda_newcases + b_deaths*(att8deaths1^lambda_deaths-1)/lambda_deaths + b_cost*att9costs1
V[["PackageB"]] = asc_right + b_informationcampaign_no*(att1information_campaigns2==1) + b_informationcampaign_yes*(att1information_campaigns2==2) +
b_solarbed_ban_no*(att2prohibiting_solarbeds2==1) + b_solarbed_ban_yes*(att2prohibiting_solarbeds2==2) +
b_solarium_ban_no*(att3prohibiting_solarstudios2==1) + b_solarium_ban_yes*(att3prohibiting_solarstudios2==2) +
b_lowerprice_sunscreen_no*(att4pricereduction_sunscreen2==1) + b_lowerprice_sunscreen_yes*(att4pricereduction_sunscreen2==2) +
b_free_sunscreen_no*(att5free_sunscreen2==1) + b_free_sunscreen_yes*(att5free_sunscreen2==2) +
b_free_app_no*(att6app2==1) + b_free_app_yes*(att6app2==2) +
b_newcases*(att7new_cases2^lambda_newcases-1)/lambda_newcases + b_deaths*(att8deaths2^lambda_deaths-1)/lambda_deaths + b_cost*att9costs2
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(PackageA=0, PackageB=1),
avail = 1,
choiceVar = Y2,
utilities = V
)
### Compute probabilities for MNL model component
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 model
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
###Output
Estimation method : bgw
Model diagnosis : Relative function convergence
Optimisation diagnosis : Maximum found
hessian properties : Negative definite
maximum eigenvalue : -2.775176
reciprocal of condition number : 0.00396026
Number of individuals : 793
Number of rows in database : 9516
Number of modelled outcomes : 9516
Number of cores used : 4
Number of inter-individual draws : 5000 (sobol)
LL(start) : -6092.37
LL at equal shares, LL(0) : -6595.99
LL at observed shares, LL(C) : -6593.54
LL(final) : -5508.67
Rho-squared vs equal shares : 0.1648
Adj.Rho-squared vs equal shares : 0.1615
Rho-squared vs observed shares : 0.1645
Adj.Rho-squared vs observed shares : 0.1613
AIC : 11061.34
BIC : 11218.88
Estimated parameters : 22
Time taken (hh:mm:ss) : 13:35:52.88
pre-estimation : 02:09:2.47
estimation : 01:44:47.57
initial estimation : 01:32:17.58
estimation after rescaling : 00:12:29.99
post-estimation : 09:42:2.83
Iterations : 20
initial estimation : 19
estimation after rescaling : 1
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) p(2-sided)
asc_left 0.00000 NA NA NA
asc_right_mu -0.13915 0.04234 -3.2865 0.001014
asc_right_sig 0.60248 0.05426 11.1046 0.000000
b_informationcampaign_no 0.00000 NA NA NA
b_informationcampaign_yes_mu 0.66858 0.06306 10.6028 0.000000
b_informationcampaign_yes_sig 0.89721 0.08185 10.9615 0.000000
b_solarbed_ban_no 0.00000 NA NA NA
b_solarbed_ban_yes_mu 0.04378 0.06225 0.7033 0.481847
b_solarbed_ban_yes_sig 0.81315 0.08498 9.5690 0.000000
b_solarium_ban_no 0.00000 NA NA NA
b_solarium_ban_yes_mu 0.03001 0.06719 0.4466 0.655167
b_solarium_ban_yes_sig 1.01701 0.08436 12.0550 0.000000
b_lowerprice_sunscreen_no 0.00000 NA NA NA
b_lowerprice_sunscreen_yes_mu 0.50803 0.05722 8.8785 0.000000
b_lowerprice_sunscreen_yes_sig -0.60792 0.09560 -6.3592 2.028e-10
b_free_sunscreen_no 0.00000 NA NA NA
b_free_sunscreen_yes_mu 0.17983 0.05930 3.0324 0.002426
b_free_sunscreen_yes_sig 0.61790 0.09914 6.2326 4.587e-10
b_free_app_no 0.00000 NA NA NA
b_free_app_yes_mu 0.43040 0.05358 8.0323 8.882e-16
b_free_app_yes_sig 0.67710 0.08010 8.4529 0.000000
b_log_newcases_mu -3.04217 0.46284 -6.5728 4.937e-11
b_log_newcases_sig 1.58243 0.16345 9.6814 0.000000
lambda_newcases 1.71702 0.40339 4.2565 2.076e-05
b_log_deaths_mu -3.58676 0.51303 -6.9913 2.724e-12
b_log_deaths_sig 2.53755 0.24172 10.4977 0.000000
lambda_deaths 0.72174 0.38281 1.8854 0.059378
b_log_cost_mu -2.01922 0.11949 -16.8984 0.000000
b_log_cost_sig 2.74803 0.14313 19.1993 0.000000
Rob.s.e. Rob.t.rat.(0) p(2-sided)
asc_left NA NA NA
asc_right_mu 0.04270 -3.2587 0.001119
asc_right_sig 0.06628 9.0905 0.000000
b_informationcampaign_no NA NA NA
b_informationcampaign_yes_mu 0.06599 10.1317 0.000000
b_informationcampaign_yes_sig 0.08583 10.4537 0.000000
b_solarbed_ban_no NA NA NA
b_solarbed_ban_yes_mu 0.06019 0.7274 0.466996
b_solarbed_ban_yes_sig 0.08630 9.4218 0.000000
b_solarium_ban_no NA NA NA
b_solarium_ban_yes_mu 0.06528 0.4596 0.645786
b_solarium_ban_yes_sig 0.09033 11.2592 0.000000
b_lowerprice_sunscreen_no NA NA NA
b_lowerprice_sunscreen_yes_mu 0.06165 8.2401 2.220e-16
b_lowerprice_sunscreen_yes_sig 0.10441 -5.8226 5.794e-09
b_free_sunscreen_no NA NA NA
b_free_sunscreen_yes_mu 0.06093 2.9516 0.003161
b_free_sunscreen_yes_sig 0.10315 5.9903 2.095e-09
b_free_app_no NA NA NA
b_free_app_yes_mu 0.05757 7.4756 7.683e-14
b_free_app_yes_sig 0.08543 7.9258 2.220e-15
b_log_newcases_mu 0.46651 -6.5211 6.981e-11
b_log_newcases_sig 0.16047 9.8611 0.000000
lambda_newcases 0.41196 4.1679 3.074e-05
b_log_deaths_mu 0.50448 -7.1098 1.162e-12
b_log_deaths_sig 0.23884 10.6247 0.000000
lambda_deaths 0.43479 1.6600 0.096918
b_log_cost_mu 0.13258 -15.2306 0.000000
b_log_cost_sig 0.14798 18.5704 0.000000
Overview of choices for MNL model component :
PackageA PackageB
Times available 9516.00 9516.00
Times chosen 4866.00 4650.00
Percentage chosen overall 51.13 48.87
Percentage chosen when available 51.13 48.87
1) I’ve already tried to derive the MRS between each of the (normally distributed) policy attributes and the cost attribute. This seems to work for most attributes, but not for all. In the code provided below, this applies to ‘Price sunscreen 30% lower’. When I run the code for this attribute, the values of all the draws for ‘b_lowerprice_sunscreen’ are ‘NaN’ and, as a result, also those for ‘wtp_lowerprice_sunscreen’.
Code: Select all
### WTP estimates policies
model = apollo_loadModel("MMNL_skincancer_AT_PolicyNormal_EffectsMuShiftedLogNorm_MNLstartingvalues_linearTax_5000draws_BoxCox_Effectiveness")
apollo_modelOutput(model, modelOutput_settings=list(printPVal=2))
#information campaigns
N = 1000
b_informationcampaign = rnorm(N, mean=model$estimate["b_informationcampaign_yes_mu"], sd=model$estimate["b_informationcampaign_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_informationcampaign = b_informationcampaign/b_cost
mean(wtp_informationcampaign); sd(wtp_informationcampaign)
#prohibition of tanning beds
N = 1000
b_solarbed_ban = rnorm(N, mean=model$estimate["b_solarbed_ban_yes_mu"], sd=model$estimate["b_solarbed_ban_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_solarbed_ban = b_solarbed_ban/b_cost
mean(wtp_solarbed_ban); sd(wtp_solarbed_ban)
#prohibition of solaria
N = 1000
b_solarium_ban = rnorm(N, mean=model$estimate["b_solarium_ban_yes_mu"], sd=model$estimate["b_solarium_ban_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_solarium_ban = b_solarium_ban/b_cost
mean(wtp_solarium_ban); sd(wtp_solarium_ban)
#price sunscreen 30% lower
N = 1000
b_lowerprice_sunscreen = rnorm(N, mean=model$estimate["b_lowerprice_sunscreen_yes_mu"], sd=model$estimate["b_lowerprice_sunscreen_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_lowerprice_sunscreen = b_lowerprice_sunscreen/b_cost
mean(wtp_lowerprice_sunscreen); sd(wtp_lowerprice_sunscreen)
#free sunscreen in public areas
N = 1000
b_free_sunscreen = rnorm(N, mean=model$estimate["b_free_sunscreen_yes_mu"], sd=model$estimate["b_free_sunscreen_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_free_sunscreen = b_free_sunscreen/b_cost
mean(wtp_free_sunscreen); sd(wtp_free_sunscreen)
#free skin cancer detection app
N = 1000
b_free_app = rnorm(N, mean=model$estimate["b_free_app_yes_mu"], sd=model$estimate["b_free_app_yes_sig"])
b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
wtp_free_app = b_free_app/b_cost
mean(wtp_free_app); sd(wtp_free_app)
#output - Warning message
N = 1000
> b_lowerprice_sunscreen = rnorm(N, mean=model$estimate["b_lowerprice_sunscreen_yes_mu"], sd=model$estimate["b_lowerprice_sunscreen_yes_sig"])
Warning message:
In rnorm(N, mean = model$estimate["b_lowerprice_sunscreen_yes_mu"], :
NAs produced
> b_cost = -exp(model$estimate["b_log_cost_mu"])-exp(rnorm(N, mean=model$estimate["b_log_cost_mu"], sd=model$estimate["b_log_cost_sig"]))
> wtp_lowerprice_sunscreen = b_lowerprice_sunscreen/b_cost
> mean(wtp_lowerprice_sunscreen); sd(wtp_lowerprice_sunscreen)
[1] NaN
[1] NA
2) For the derivation of the MRS between the two (positive lognormally distributed) effectiveness attributes and the cost attribute, I’m wondering how to adapt the code to reflect the Box-Cox transformation of these attributes. Preferably, I’d obtain MRS estimates for each of the four levels of each effectiveness attribute (as included in the DCE) and/or plot the non-linear WTP curve for different levels of the attribute. I’ve looked at the suggestions by Daly et al. (2012), but I’m not sure how to implement this in Apollo.
Thank you very much in advance for your reply!
Kind regards,
Sander