In Apollo manual (p20), you indicate that correlations between MMNL random parameters can be calculated empirically from the draws produced by apollo_unconditionals. You you please provide the code to get that ? I could not find into the manual.
Thanks
Marianne
My code for the MMNL with 1 lognormal and 6 normal coefficients is:
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
#library(xtable)
### Initialise code
apollo_initialise()
### set working directory
setwd("C:/Users/marianne.lefebvre/Nextcloud/RECHERCHE/BEHAVE/BEHAVEPartagePauline/EAU/DCE_Eau/Data_analysis/Rapollo_RevisionERAE")
### Set core controls
apollo_control = list(
modelName = "MMNL_DCEEau_WTPspace",
modelDescr = "MMNL on Behave water DCE data in WTP space",
indivID = "idLS",
nCores = 3,
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
database = read.csv("Final_Result_Wide.csv",header=TRUE)
settings<-list(
alternative_column="Quchoice",
alternative_specific_attributes=c("DeclaFarm","DeclaPlotSoil","DeclaAuto","ManualWeekly","Smart","Vol"),
choice_column="Choice",
ID_column="idLS",
observation_column="Choicenumber"
)
#database = apollo_longToWide(databaseLong2,longToWide_settings=settings)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
#Estimates MNL
# Estimate s.e. t.rat.(0)
#ASC 0.96904 0.242249 4.0002
#b_Vol 0.01871 0.002577 7.2585
#b_DeclaFarm 0.66822 0.138640 4.8198
#b_DeclaPlotSoil 0.85329 0.170235 5.0124
#b_DeclaAuto 0.31873 0.168379 1.8929
#b_ManualWeekly -0.12685 0.230557 -0.5502
#b_Smart -0.20943 0.220074 -0.9516
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(ASC_mu = 13.4403,
ASC_sig = 123.2308,
v_Vol_mu= -2.6916 ,
v_Vol_sig=0.7512,
v_DeclaFarm_mu = 47.3078,
v_DeclaFarm_sig = 29.6476,
v_DeclaPlotSoil_mu = 47.7545,
v_DeclaPlotSoil_sig = 10.1956,
v_DeclaAuto_mu = 26.4988,
v_DeclaAuto_sig = -30.1430,
v_ManualWeekly_mu = -20.3407,
v_ManualWeekly_sig = -53.2512,
v_Smart_mu = -15.6283,
v_Smart_sig = 31.1956)
### 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()
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws: inter because we have panel data with several choices per ind. we want heterogeneity at ind level (not choice level)
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 2000,
interNormDraws = c("draws_DeclaFarm","draws_DeclaPlotSoil","draws_DeclaAuto","draws_ManualWeekly","draws_Smart","draws_Vol","draws_ASC" )
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_Vol"]] = exp(v_Vol_mu + v_Vol_sig * draws_Vol )
randcoeff[["v_DeclaFarm"]] = v_DeclaFarm_mu + v_DeclaFarm_sig * draws_DeclaFarm
randcoeff[["v_DeclaPlotSoil"]] = v_DeclaPlotSoil_mu + v_DeclaPlotSoil_sig * draws_DeclaPlotSoil
randcoeff[["v_DeclaAuto"]] = v_DeclaAuto_mu + v_DeclaAuto_sig * draws_DeclaAuto
randcoeff[["v_ManualWeekly"]] = v_ManualWeekly_mu + v_ManualWeekly_sig * draws_ManualWeekly
randcoeff[["v_Smart"]] = v_Smart_mu + v_Smart_sig * draws_Smart
randcoeff[["ASC"]] = ASC_mu + ASC_sig * draws_ASC
return(randcoeff)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### 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[["alt1"]] = (v_DeclaFarm * DeclaFarm1 + v_DeclaPlotSoil * DeclaPlotSoil1 + v_DeclaAuto * DeclaAuto1 + v_ManualWeekly * ManualWeekly1 + v_Smart * Smart1 + Vol1)*b_Vol
V[["alt2"]] = (v_DeclaFarm * DeclaFarm2 + v_DeclaPlotSoil * DeclaPlotSoil2 + v_DeclaAuto * DeclaAuto2 + v_ManualWeekly * ManualWeekly2 + v_Smart * Smart2 + Vol2)*b_Vol
V[["alt3"]] = (ASC + v_DeclaFarm * DeclaFarm3 + v_DeclaPlotSoil * DeclaPlotSoil3 + v_DeclaAuto * DeclaAuto3 + v_ManualWeekly * ManualWeekly3 + v_Smart * Smart3 + Vol3)*b_Vol
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = Choice,
utilities = 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 = c(1,3,5,7),
# nRep = 30
#)
#apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model,modelOutput_settings = list(printPVal=1)) # for one sided t-test
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
apollo_saveOutput(model,saveOutput_settings = list(printPVal=1)) # for one sided t-test
#As expected, the estimated coefficients for data sharing are positive (interpreted as a decrease in WTA: They ask for less compensation in terms of prioritized water volumes)
#The effect size is largely higher for water data reporting (farmers ask for more/less ? compensation in terms of )
#EXPORT TABLE TEX: commande à vérifier
#print(xtable(model,digits = 2,caption="RPL in WTA space",include.rownames = F)
# ################################################################# #
##### POST-PROCESSING ####
# ################################################################# #
### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
apollo_sink()
beta= apollo_unconditionals(model,apollo_probabilities,apollo_inputs)