diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_not_treat.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_not_treat.R new file mode 100644 index 0000000000000000000000000000000000000000..6c0a05e148f09d9902c16f7ac74f8fe7bf80f495 --- /dev/null +++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_not_treat.R @@ -0,0 +1,136 @@ +#### 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) + + diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R index 8177fbe1d15e9860278bf60372ff976de482cb67..aa19ccd7a08773b25c5a95804bad0f76f03685d0 100644 --- a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R +++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R @@ -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")) diff --git a/Scripts/visualize_distr_spli.R b/Scripts/visualize_distr_spli.R index ad73703cd1df9ba1dc1b0c3750653105156637a3..ecce6d20436f3637ecb0fbbec1d07ca585869934 100644 --- a/Scripts/visualize_distr_spli.R +++ b/Scripts/visualize_distr_spli.R @@ -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