################################################################### # # Script for Replication of: Consumers' Preference and Willingness to Pay # For Biofortified Garri New Design Pretest # Description: The script to estimate Stated ANA mixed logit model for ## nutrition + method treatment ## Author: Akinwehinmi Titilayo ####################################################### # ################################################################# # #### LOAD LIBRARY AND DEFINE CORE SETTINGS #### # ################################################################# # ### Clear memory #rm(list = ls()) ### Load Apollo library check.packages <- function(pkg){ new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])] if (length(new.pkg)) install.packages(new.pkg, dependencies = TRUE) sapply(pkg, require, character.only = TRUE) } packages <- c("apollo") check.packages(packages) ### Initialise code apollo_initialise() ### Set core controls apollo_control = list( modelName = "Stated_ANA_Method", modelDescr = "Stated ANA Model Method Treatment", indivID = "id", mixing = TRUE, nCores = 7, analyticGrad = TRUE, # Needs to be set explicitly for models with inter-intra draws (requires extra RAM). outputDirectory = "output" ) # ################################################################# # #### LOAD DATA AND APPLY ANY TRANSFORMATIONS #### # ################################################################# # #Subset database for Control and Treatment 1 #apollo_database = read.csv("data_biofortificaton.csv") database <- subset(apollo_database, consumer==3) #database <- subset(database, id!=199&id!=204) database <- subset(database, (is.na(database$stated_method_att)==F)) database <- subset(database, (is.na(database$stated_va_att)==FALSE)) database <- subset(database, (is.na(database$stated_iron_att)==FALSE)) database <- subset(database, (is.na(database$stated_starch_att)==FALSE)) database <- subset(database, (is.na(database$stated_price_att)==FALSE)) # ################################################################# # #### DEFINE MODEL PARAMETERS #### # ################################################################# # ### Vector of parameters, including any that are kept fixed in estimation apollo_beta = c( #Parameters b_price = 0, #sigma_price_inter = 0, asc3 = 0, sigma_asc3_inter = 0, va = 0, sigma_va_inter = 0, iron25 = 0, sigma_iron25_inter = 0, iron50 = 0, sigma_iron50_inter = 0, lowstarch = 0, sigma_lowstarch_inter = 0, gm = 0, sigma_gm_inter = 0, ge = 0, sigma_ge_inter = 0, price_ana = 0, va_ana = 0, iron25_ana = 0, iron50_ana = 0, lowstarch_ana = 0, gm_ana = 0, ge_ana = 0, # # Correlation in preferences within the control group sigma_va_iron25 = 0, sigma_va_iron50 = 0, sigma_va_lowstarch = 0, sigma_va_gm = 0, sigma_va_ge = 0, sigma_iron25_iron50 = 0, sigma_iron25_lowstarch = 0, sigma_iron25_gm = 0, sigma_iron25_ge = 0, sigma_iron50_lowstarch = 0, sigma_iron50_gm = 0, sigma_iron50_ge = 0, sigma_lowstarch_gm = 0, sigma_lowstarch_ge = 0, sigma_gm_ge = 0, sigma_va_asc3 = 0, sigma_iron25_asc3 = 0, sigma_iron50_asc3 = 0, sigma_lowstarch_asc3 = 0, sigma_gm_asc3 = 0, sigma_ge_asc3 = 0 ) ### 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() ### Overwrite beta starting values #apollo_readBeta(apollo_beta, apollo_fixed, "MMNL_Preference_Corr_Ondo_9", overwriteFixed = FALSE) # ################################################################# # #### DEFINE RANDOM COMPONENTS #### # ################################################################# # ### Set parameters for generating draws apollo_draws = list( interDrawsType = "sobol", interNDraws = 500, interUnifDraws = c("draws_price_inter"), interNormDraws = c( "draws_asc3_inter", "draws_va_inter", "draws_iron25_inter", "draws_iron50_inter", "draws_lowstarch_inter", "draws_gm_inter", "draws_ge_inter", #draws for correlation "draws_va_iron25" , "draws_va_iron50" , "draws_va_lowstarch" , "draws_va_gm" , "draws_va_ge" , "draws_iron25_iron50" , "draws_iron25_lowstarch" , "draws_iron25_gm" , "draws_iron25_ge" , "draws_iron50_lowstarch" , "draws_iron50_gm" , "draws_iron50_ge" , "draws_lowstarch_gm" , "draws_lowstarch_ge" , "draws_gm_ge" , "draws_va_asc3" , "draws_iron25_asc3" , "draws_iron50_asc3" , "draws_lowstarch_asc3" , "draws_gm_asc3" , "draws_ge_asc3" ) ) ### Create random parameters apollo_randCoeff = function(apollo_beta, apollo_inputs){ randcoeff = list() #randcoeff[["b_price"]] = -exp(price + sigma_price_inter*draws_price_inter) randcoeff[["b_asc3"]] = (asc3 + sigma_asc3_inter*draws_asc3_inter) randcoeff[["b_va"]] = (va + sigma_va_inter*draws_va_inter + sigma_va_asc3*draws_va_asc3) randcoeff[["b_iron25"]] = (iron25 + sigma_iron25_inter*draws_iron25_inter + sigma_va_iron25*draws_va_iron25 + sigma_iron25_asc3*draws_iron25_asc3) randcoeff[["b_iron50"]] = (iron50 + sigma_iron50_inter*draws_iron50_inter + sigma_iron50_asc3*draws_iron50_asc3 + sigma_va_iron50*draws_va_iron50 + sigma_iron25_iron50*draws_iron25_iron50) randcoeff[["b_lowstarch"]] = (lowstarch + sigma_lowstarch_inter*draws_lowstarch_inter + sigma_va_lowstarch*draws_va_lowstarch + sigma_iron25_lowstarch*draws_iron25_lowstarch + sigma_lowstarch_asc3*draws_lowstarch_asc3 + sigma_iron50_lowstarch*draws_iron50_lowstarch) randcoeff[["b_gm"]] = (gm + sigma_gm_inter*draws_gm_inter + sigma_va_gm*draws_va_gm + sigma_gm_asc3*draws_gm_asc3 + sigma_iron25_gm*draws_iron25_gm + sigma_iron50_gm*draws_iron50_gm + sigma_lowstarch_gm*draws_lowstarch_gm) randcoeff[["b_ge"]] = (ge + sigma_ge_inter*draws_ge_inter + sigma_ge_asc3*draws_ge_asc3 + sigma_va_ge*draws_va_ge + sigma_iron25_ge*draws_iron25_ge + sigma_iron50_ge*draws_iron50_ge + sigma_lowstarch_ge*draws_lowstarch_ge + sigma_gm_ge*draws_gm_ge) 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"]] = (b_price*stated_price_att1 + b_va*stated_vab_att1 + b_iron25*stated_iron25_att1 + b_iron50*stated_iron50_att1 + b_lowstarch*stated_lowstarch_att1 + b_gm*stated_gm_att1 + b_ge*stated_ge_att1 + price_ana*stated_price_natt1 + va_ana*stated_vab_natt1 + iron25_ana*stated_iron25_natt1 + iron50_ana*stated_iron50_natt1 + lowstarch_ana*stated_lowstarch_natt1 + gm_ana*stated_gm_natt1 + ge_ana*stated_ge_natt1) V[["alt2"]] = (b_price*stated_price_att2 + b_va*stated_vab_att2 + b_iron25*stated_iron25_att2 + b_iron50*stated_iron50_att2 + b_lowstarch*stated_lowstarch_att2 + b_gm*stated_gm_att2 + b_ge*stated_ge_att2 + price_ana*stated_price_natt2 + va_ana*stated_vab_natt2 + iron25_ana*stated_iron25_natt2 + iron50_ana*stated_iron50_natt2 + lowstarch_ana*stated_lowstarch_natt2 + gm_ana*stated_gm_natt2 + ge_ana*stated_ge_natt2) V[["alt3"]] = (b_price*stated_price_att3 + price_ana*stated_price_natt3 + b_asc3*asc33) ### 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) ### Average across intra-individual draws # P = apollo_avgIntraDraws(P, apollo_inputs, 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 #### # ################################################################# # ### optionally load pre-estimated model from memory # model = apollo_loadModel("MMNL_Preference_Batched_Corr_Ondo_2") # apollo_addCovariance(model, apollo_inputs) ### Optional speedTest # speedTest_settings=list( # nDrawsTry = c(100, 200, 500, 1000), # nCoresTry = 8, # nRep = 10) # apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings) ### Optional starting values search model = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs, estimate_settings=list(maxIterations=400)) # ################################################################# # #### MODEL OUTPUTS #### # ################################################################# # # ----------------------------------------------------------------- # #---- FORMATTED OUTPUT (TO SCREEN) ---- # ----------------------------------------------------------------- # modelOutput_settings = list( printOutliers = 40, printPVal = 2, printClassical = FALSE ) apollo_modelOutput(model) # ----------------------------------------------------------------- # #---- FORMATTED OUTPUT (TO FILE, using model name) ---- # ----------------------------------------------------------------- # apollo_saveOutput(model, modelOutput_settings) # ################################################################# # ##### POST-PROCESSING #### # ################################################################# # ###Compute WTP from the Preference space estimates deltaMethod_settings=list(expression=c(wtp_va="va/(-b_price)", wtp_iron25="iron25/(-b_price)", wtp_iron50="iron50/(-b_price)", wtp_lowstarch="lowstarch/(-b_price)", wtp_gm="gm/(-b_price)", wtp_ge="ge/(-b_price)", wtp_sd_va="sigma_va_inter/(-b_price)", wtp_sd_iron25="sigma_iron25_inter/(-b_price)", wtp_sd_iron50="sigma_iron50_inter/(-b_price)", wtp_sd_lowstarch="sigma_lowstarch_inter/(-b_price)", wtp_sd_gm="sigma_gm_inter/(-b_price)", wtp_sd_ge="sigma_ge_inter/(-b_price)")) wtp_method <- apollo_deltaMethod(model, deltaMethod_settings) write.csv(wtp_method, "wtp_stated_method.csv")