Skip to content
Snippets Groups Projects
Select Git revision
  • 05d9bfb49fc92353343ff9ba6a49de79add8b797
  • main default protected
  • test_coef
  • 21-things-to-take-care-of-before-submission
4 results

mxl_wtp_space_caseG.R

Blame
  • Nino Cavallaro's avatar
    nc71qaxa authored
    05d9bfb4
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    mxl_wtp_space_caseG.R 7.80 KiB
    #### Apollo standard script #####
    
    library(apollo) # Load apollo package 
    
    
    
    # Test treatment effect
    
    database_full$probabilities <- logit_sig %>% predict(database_full, type = "response")
    database_full$predicted.classes <- ifelse(database_full$probabilities  >= cut_off, 1, 0)
    
    database <- database_full %>%
      filter(!is.na(Treatment_new)) %>%
      mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
             Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
             Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
             Dummy_pred_Treated = case_when(Treatment_new == 6 & predicted.classes == 1 ~ 1, TRUE~0),
             Dummy_Treated_predicted = case_when(Treatment_new == 1  & predicted.classes == 1|Treatment_new == 2 & predicted.classes == 1 ~ 1, TRUE ~ 0),
             Dummy_Treated_notpredicted = case_when(Treatment_new == 1  & predicted.classes == 0|Treatment_new == 2 & predicted.classes == 0 ~ 1, TRUE ~ 0)) 
    
    table(database$Dummy_Treated)
    table(database$Dummy_Vol_Treated)
    table(database$Dummy_no_info)
    
    #initialize model 
    
    apollo_initialise()
    
    
    ### Set core controls
    apollo_control = list(
      modelName  = "MXL_wtp_Case_G base",
      modelDescr = "MXL wtp space Case G base",
      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_ASC_sq_treated_predicted = 0,
                  mu_ASC_sq_treated_notpredicted = 0,
                  mu_ASC_sq_vol_treated = 0,
                  mu_ASC_sq_no_info = 0,
                  mu_nat_treated_predicted =0,
                  mu_nat_treated_notpredicted = 0,
                  mu_nat_vol_treated = 0,
                  mu_nat_no_info = 0,
                  mu_walking_treated_predicted =0,
                  mu_walking_treated_notpredicted =0,
                  mu_walking_vol_treated = 0,
                  mu_walking_no_info = 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)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
                                    +  mu_nat_treated_predicted * Naturalness_1 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_1 * Dummy_no_info
                                  +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_1 *Dummy_Treated_notpredicted
                                  +  mu_walking_treated_predicted * WalkingDistance_1 *Dummy_Treated_predicted +  mu_walking_treated_notpredicted * WalkingDistance_1 *Dummy_Treated_notpredicted
                                  + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
                                  +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1)
      
      V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
                                  + mu_nat_treated_predicted * Naturalness_2 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_2 * Dummy_no_info
                                  + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_2 *Dummy_Treated_notpredicted
                                  + mu_walking_treated_predicted * WalkingDistance_2 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_2 *Dummy_Treated_notpredicted
                                  + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
                                  + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2)
      
      V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
                                  +  mu_nat_treated_predicted * Naturalness_3 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_3 * Dummy_no_info
                                  +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_3 *Dummy_Treated_notpredicted
                                  +  mu_walking_treated_predicted * WalkingDistance_3 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_3 *Dummy_Treated_notpredicted
                                  +  mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
                                  +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
                                  +  mu_ASC_sq_treated_predicted * Dummy_Treated_predicted + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
                                  +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3 + mu_ASC_sq_treated_notpredicted * Dummy_Treated_notpredicted)
      
      
      ### 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_case_g = 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_case_g)