From 1bb5187c3b8972e6617f3627dbc753e63a7d3f47 Mon Sep 17 00:00:00 2001 From: Cavallaro <nino.cavallaro@idiv.de> Date: Wed, 6 Mar 2024 09:14:12 +0100 Subject: [PATCH] Case C NR Rent INT --- Scripts/compare_split_samples.R | 16 +- .../Split_samples/mxl_wtp_space_not_tr_A.R | 8 +- .../mxl/Split_samples/mxl_wtp_space_tr_A.R | 8 +- .../Split_samples/mxl_wtp_space_vol_tr_A.R | 8 +- Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R | 193 ++++++++++++++++++ 5 files changed, 216 insertions(+), 17 deletions(-) create mode 100644 Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R diff --git a/Scripts/compare_split_samples.R b/Scripts/compare_split_samples.R index fa39922..011aa55 100644 --- a/Scripts/compare_split_samples.R +++ b/Scripts/compare_split_samples.R @@ -23,8 +23,8 @@ colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", mxl_melt_info <- melt(mxl_info_compare[1:4], id = "Coefficent") mxl_melt_info$ME <- mxl_info_compare$Margin_of_error_TR -mxl_melt_info$ME[8:14] <- mxl_info_compare$Margin_of_error_VOL_TR -mxl_melt_info$ME[15:21] <- mxl_info_compare$Margin_of_error_NOT_TR +mxl_melt_info$ME[9:16] <- mxl_info_compare$Margin_of_error_VOL_TR +mxl_melt_info$ME[17:24] <- mxl_info_compare$Margin_of_error_NOT_TR # Figure paper compare info treatments # @@ -60,16 +60,16 @@ colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR") -mxl_melt_info <- melt(mxl_info_compare[1:4], id = "Coefficent") -mxl_melt_info$ME <- mxl_info_compare$Margin_of_error_TR -mxl_melt_info$ME[8:14] <- mxl_info_compare$Margin_of_error_VOL_TR -mxl_melt_info$ME[15:21] <- mxl_info_compare$Margin_of_error_NOT_TR +mxl_melt_info_log <- melt(mxl_info_compare[1:4], id = "Coefficent") +mxl_melt_info_log$ME <- mxl_info_compare$Margin_of_error_TR +mxl_melt_info_log$ME[8:14] <- mxl_info_compare$Margin_of_error_VOL_TR +mxl_melt_info_log$ME[15:21] <- mxl_info_compare$Margin_of_error_NOT_TR -mxl_melt_info <- mxl_melt_info %>% mutate(ME = case_when(abs(value) < 1 ~ ME*100, TRUE~ME), +mxl_melt_info_log <- mxl_melt_info_log %>% mutate(ME = case_when(abs(value) < 1 ~ ME*100, TRUE~ME), value = case_when(abs(value) < 1 ~ value*100, TRUE~value)) # Figure paper compare info treatments # -ggplot(data=mxl_melt_info, aes(x=Coefficent, y=abs(value), fill=variable)) + +ggplot(data=mxl_melt_info_log, aes(x=Coefficent, y=abs(value), fill=variable)) + geom_bar(stat="identity", position='dodge', width = 0.9) + geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) + ylab("Absolute Value") + diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R index 546f7c8..efe4601 100644 --- a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R +++ b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R @@ -39,7 +39,8 @@ database <- database_full %>% ASC_sq = 0, sig_natural = 15, sig_walking = 2, - sig_rent = 2) + sig_rent = 2, + sig_ASC_sq = 0) ### specify parameters that should be kept fixed, here = none apollo_fixed = c() @@ -49,7 +50,7 @@ database <- database_full %>% interDrawsType = "sobol", interNDraws = n_draws, interUnifDraws = c(), - interNormDraws = c("draws_natural", "draws_walking", "draws_rent"), + interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"), intraDrawsType = "halton", intraNDraws = 0, intraUnifDraws = c(), @@ -63,6 +64,7 @@ database <- database_full %>% 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) } @@ -90,7 +92,7 @@ database <- database_full %>% V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 - Rent_2) - V[['alt3']] = -b_mu_rent*(ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - + V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - Rent_3) diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R index 62564b8..d238b0d 100644 --- a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R +++ b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R @@ -39,7 +39,8 @@ database <- database_full %>% ASC_sq = 0, sig_natural = 15, sig_walking = 2, - sig_rent = 2) + sig_rent = 2, + sig_ASC_sq = 2) ### specify parameters that should be kept fixed, here = none apollo_fixed = c() @@ -49,7 +50,7 @@ database <- database_full %>% interDrawsType = "sobol", interNDraws = n_draws, interUnifDraws = c(), - interNormDraws = c("draws_natural", "draws_walking", "draws_rent"), + interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"), intraDrawsType = "halton", intraNDraws = 0, intraUnifDraws = c(), @@ -63,6 +64,7 @@ database <- database_full %>% 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) } @@ -90,7 +92,7 @@ database <- database_full %>% V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 - Rent_2) - V[['alt3']] = -b_mu_rent*(ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - + V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - Rent_3) diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R index 6fc7212..52d1642 100644 --- a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R +++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R @@ -39,7 +39,8 @@ database <- database_full %>% ASC_sq = 0, sig_natural = 15, sig_walking = 2, - sig_rent = 2) + sig_rent = 2, + sig_ASC_sq = 0) ### specify parameters that should be kept fixed, here = none apollo_fixed = c() @@ -49,7 +50,7 @@ database <- database_full %>% interDrawsType = "sobol", interNDraws = n_draws, interUnifDraws = c(), - interNormDraws = c("draws_natural", "draws_walking", "draws_rent"), + interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"), intraDrawsType = "halton", intraNDraws = 0, intraUnifDraws = c(), @@ -63,6 +64,7 @@ database <- database_full %>% 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) } @@ -90,7 +92,7 @@ database <- database_full %>% V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 - Rent_2) - V[['alt3']] = -b_mu_rent*(ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - + V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 - Rent_3) diff --git a/Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R b/Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R new file mode 100644 index 0000000..b7cd46a --- /dev/null +++ b/Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R @@ -0,0 +1,193 @@ +#### Apollo standard script ##### + +library(apollo) # Load apollo package + + + +# Test treatment effect + +database <- database_full %>% + filter(!is.na(Treatment_new)) %>% + mutate(Dummy_Video_1 = case_when(Treatment_new == 1 ~ 1, TRUE ~ 0), + Dummy_Video_2 = case_when(Treatment_new == 5 ~ 1, TRUE ~ 0), + Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0), + Dummy_Info_nv1 = case_when(Treatment_new == 2 ~1, TRUE~0), + Dummy_Info_nv2 = case_when(Treatment_new == 4 ~1 , TRUE~0)) + +#initialize model + +apollo_initialise() + + +### Set core controls +apollo_control = list( + modelName = "MXL_wtp_NR_Case_C Rent INT", + modelDescr = "MXL wtp space NR Case C Rent INT", + indivID ="id", + mixing = TRUE, + HB= FALSE, + nCores = n_cores, + outputDirectory = "Estimation_results/mxl" +) + +##### 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, + mu_nat_NR = 0, + mu_wd_NR = 0, + mu_asc_NR = 0, + mu_ASC_sq_vid1 = 0, + mu_ASC_sq_vid2 = 0, + mu_ASC_sq_no_info = 0, + mu_ASC_sq_info_nv1 = 0, + mu_ASC_sq_info_nv2 = 0, + mu_rent_vid1 = 0, + mu_rent_vid2 = 0, + mu_rent_no_info = 0, + mu_rent_info_nv1 = 0, + mu_rent_info_nv2 = 0, + mu_nat_vid1 =0, + mu_nat_vid2 = 0, + mu_nat_no_info = 0, + mu_nat_info_nv1 = 0, + mu_nat_info_nv2 = 0, + mu_walking_vid1 =0, + mu_walking_vid2 = 0, + mu_walking_no_info = 0, + mu_walking_info_nv1 = 0, + mu_walking_info_nv2 = 0, + sig_natural = 15, + sig_walking = 2, + sig_rent = 2, + sig_ASC_sq = 2) + +### 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 + mu_rent_vid1 *Dummy_Video_1 + mu_rent_vid2 * Dummy_Video_2 + mu_rent_no_info * Dummy_no_info + + mu_rent_info_nv1 * Dummy_Info_nv1 + mu_rent_info_nv2 * Dummy_Info_nv2) * + (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 + + mu_nat_NR * Naturalness_1 * Z_Mean_NR + mu_wd_NR * WalkingDistance_1 * Z_Mean_NR + + mu_nat_vid1 * Naturalness_1 *Dummy_Video_1 + mu_nat_no_info * Naturalness_1 * Dummy_no_info + + mu_nat_info_nv1 * Naturalness_1 *Dummy_Info_nv1 + mu_nat_vid2 * Naturalness_1 * Dummy_Video_2 + + mu_nat_info_nv2 * Naturalness_1 *Dummy_Info_nv2 + + mu_walking_vid1 * WalkingDistance_1 *Dummy_Video_1 + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info + + mu_walking_info_nv1 * WalkingDistance_1 *Dummy_Info_nv1 + mu_walking_vid2 * WalkingDistance_1 * Dummy_Video_2 + + mu_walking_info_nv2 * WalkingDistance_1 *Dummy_Info_nv2- Rent_1) + + V[['alt2']] = -(b_mu_rent + mu_rent_vid1 *Dummy_Video_1 + mu_rent_vid2 * Dummy_Video_2 + mu_rent_no_info * Dummy_no_info + + mu_rent_info_nv1 * Dummy_Info_nv1 + mu_rent_info_nv2 * Dummy_Info_nv2)* + (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 + + mu_nat_NR * Naturalness_2 * Z_Mean_NR + mu_wd_NR * WalkingDistance_2 * Z_Mean_NR + + mu_nat_vid1 * Naturalness_2 *Dummy_Video_1 + mu_nat_no_info * Naturalness_2 * Dummy_no_info + + mu_nat_info_nv1 * Naturalness_2 *Dummy_Info_nv1 + mu_nat_vid2 * Naturalness_2 * Dummy_Video_2 + + mu_nat_info_nv2 * Naturalness_2 *Dummy_Info_nv2+ + mu_walking_vid1 * WalkingDistance_2 *Dummy_Video_1 + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info + + mu_walking_info_nv1 * WalkingDistance_2 *Dummy_Info_nv1 + mu_walking_vid2 * WalkingDistance_2 * Dummy_Video_2 + + mu_walking_info_nv2 * WalkingDistance_2 *Dummy_Info_nv2 - Rent_2) + + V[['alt3']] = -(b_mu_rent + mu_rent_vid1 *Dummy_Video_1 + mu_rent_vid2 * Dummy_Video_2 + mu_rent_no_info * Dummy_no_info + + mu_rent_info_nv1 * Dummy_Info_nv1 + mu_rent_info_nv2 * Dummy_Info_nv2) * + (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + + mu_asc_NR * Z_Mean_NR + mu_nat_NR * Naturalness_3 * Z_Mean_NR + + mu_wd_NR * WalkingDistance_3 * Z_Mean_NR + + mu_nat_vid1 * Naturalness_3 *Dummy_Video_1 + mu_nat_no_info * Naturalness_3 * Dummy_no_info + + mu_nat_info_nv1 * Naturalness_3 *Dummy_Info_nv1 + mu_nat_vid2 * Naturalness_3 * Dummy_Video_2 + + mu_nat_info_nv2 * Naturalness_3 *Dummy_Info_nv2+ + mu_walking_vid1 * WalkingDistance_3 *Dummy_Video_1 + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info + + mu_walking_info_nv1 * WalkingDistance_3 *Dummy_Info_nv1 + mu_walking_vid2 * WalkingDistance_3 * Dummy_Video_2 + + mu_walking_info_nv2 * WalkingDistance_3 *Dummy_Info_nv2 + + mu_ASC_sq_vid1 * Dummy_Video_1 + mu_ASC_sq_vid2 * Dummy_Video_2 + + mu_ASC_sq_no_info * Dummy_no_info + mu_ASC_sq_info_nv1 * Dummy_Info_nv1 + + mu_ASC_sq_info_nv2 * Dummy_Info_nv2 - 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_NR_case_c_rentINT = apollo_estimate(apollo_beta, apollo_fixed, + apollo_probabilities, apollo_inputs, + estimate_settings=list(maxIterations=400, + estimationRoutine="bfgs", + hessianRoutine="analytic")) + + + +# ################################################################# # +#### MODEL OUTPUTS ## +# ################################################################# # +apollo_saveOutput(mxl_wtp_NR_case_c_rentINT) + + -- GitLab