Skip to content
Snippets Groups Projects
Commit 61d9e39b authored by nc71qaxa's avatar nc71qaxa
Browse files

adjust plots

parent 130b8a26
Branches
No related tags found
No related merge requests found
#### Apollo standard script #####
library(apollo) # Load apollo package
# Test treatment effect
database <- database_full %>%
filter(Treatment_D == "No Info 2")
#initialize model
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "MXL_wtp Vol_Not_Treat",
modelDescr = "MXL wtp space Vol_Not_Treat",
indivID ="id",
mixing = TRUE,
HB= FALSE,
nCores = n_cores,
outputDirectory = "Estimation_results/mxl/Split_samples"
)
##### Define model parameters depending on your attributes and model specification! ####
# set values to 0 for conditional logit model
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
ASC_sq = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
sig_ASC_sq = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
### Set parameters for generating draws, use 2000 sobol draws
apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters, define distribution of the parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
return(randcoeff)
}
### validate
apollo_inputs = apollo_validateInputs()
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 (later integrated in mnl_settings below) ####
# Define utility functions here:
V = list()
V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3),
avail = 1, # all alternatives are available in every choice
choiceVar = choice,
V = V#, # tell function to use list vector defined above
)
### 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 - nur bei Mixed Logit!
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 with bfgs algorithm
mxl_wtp_Vol_Not_Treat = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs,
estimate_settings=list(maxIterations=400,
estimationRoutine="bgw",
hessianRoutine="analytic"))
# ################################################################# #
#### MODEL OUTPUTS ##
# ################################################################# #
apollo_saveOutput(mxl_wtp_Vol_Not_Treat)
......@@ -123,7 +123,7 @@ database <- database_full %>%
mxl_wtp_Vol_Treat = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs,
estimate_settings=list(maxIterations=400,
estimationRoutine="bfgs",
estimationRoutine="bgw",
hessianRoutine="analytic"))
......
......@@ -6,6 +6,8 @@ MXL_wtp_Opt_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Sp
MXL_wtp_Treated_Not_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Not_Pred")
MXL_wtp_Treated_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Pred")
MXL_wtp_opt_treated <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Treat")
MXL_wtp_opt_not_treated <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Not_Treat")
mxl_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Treated A")
mxl_vol_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Treated A")
......@@ -18,6 +20,8 @@ models <- list(
"Control Pred" = MXL_wtp_Control_Pred_model,
"Opt Not Pred" = MXL_wtp_Opt_Not_Pred_model,
"Opt Pred" = MXL_wtp_Opt_Pred_model,
"Opt No Info" = MXL_wtp_opt_not_treated,
"Opt Info" = MXL_wtp_opt_treated,
"Treated Pred" = MXL_wtp_Treated_Pred_model,
"Control Split Sample" = mxl_not_tr,
"Opt Split Sample" = mxl_vol_tr,
......@@ -45,6 +49,8 @@ group_mapping <- c(
"Treated Pred" = "Treated",
"Opt Not Pred" = "Optional",
"Opt Pred" = "Optional",
"Opt No Info" = "Optional",
"Opt Info" = "Optional",
"Control Split Sample" = "Control",
"Opt Split Sample" = "Optional",
"Treated Split Sample" = "Treated"
......@@ -60,7 +66,9 @@ prediction_mapping <- c(
"Treated Split Sample" = "Split Sample",
"Opt Not Pred" = "Not Predicted",
"Opt Pred" = "Predicted",
"Opt Split Sample" = "Split Sample"
"Opt Split Sample" = "Split Sample",
"Opt Info" = "Predicted",
"Opt No Info" = "Not Predicted"
)
# Define the x-range for plotting (1000 values)
......@@ -128,7 +136,7 @@ calculate_z_value <- function(model_1, model_2) {
}
# Example of using the function
z_value <- calculate_z_value(MXL_wtp_Treated_Pred_model, MXL_wtp_Control_Pred_model)
z_value <- calculate_z_value(mxl_tr, mxl_vol_tr)
# Print the z-value
z_value
......@@ -198,6 +206,14 @@ for (i in seq_along(models)) {
ci_data$Group <- factor(ci_data$Group, levels = c("Treated", "Optional", "Control"))
ci_data$Model <- factor(ci_data$Model, levels = unique(ci_data$Model))
model_order <- c(
"Treated Not Pred", "Treated Pred", "Treated Split Sample",
"Opt Not Pred", "Opt Pred", "Opt No Info", "Opt Info", "Opt Split Sample",
"Control Not Pred", "Control Pred", "Control Split Sample"
)
ci_data$Model <- factor(ci_data$Model, levels = model_order)
ggplot(ci_data, aes(x = Model, y = Mean, color = Group)) +
geom_point(size = 3, position = position_dodge(width = 0.5)) + # Mean points
geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2, position = position_dodge(width = 0.5)) + # CI bars
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment