Hi,
I would recommend using an exploded logit but adding error components to introduce correlation between alternatives, approximate a nesting structure. The code below is a modification of example "EL" in the
examples page, where I introduced correlation between alternatives 1 and 2. You could so something similar with your data.
Best wishes
David
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Initialise
rm(list = ls())
library(apollo)
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "EL",
modelDescr = "Exploded logit model on drug choice data",
indivID = "ID"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = apollo_drugChoiceData
### for data dictionary, use ?apollo_drugChoiceData
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(b_brand_Artemis = 0,
b_brand_Novum = 0,
b_brand_BestValue = 0,
b_brand_Supermarket = 0,
b_brand_PainAway = 0,
b_country_CH = 0,
b_country_DK = 0,
b_country_USA = 0,
b_country_IND = 0,
b_country_RUS = 0,
b_country_BRA = 0,
b_char_standard = 0,
b_char_fast = 0,
b_char_double = 0,
b_risk = 0,
b_price = 0,
scale_1 = 1,
scale_2 = 1,
scale_3 = 1,
sigma_EC1 = 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("b_brand_Artemis", "b_country_USA", "b_char_standard", "scale_1")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs",
interNDraws = 500,
interNormDraws = c("eta1")
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list(
EC_nest1 = sigma_EC1*eta1
)
return(randcoeff)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
### 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 el_settings, order is irrelevant
V = list()
V[["alt1"]] = ( b_brand_Artemis*(brand_1=="Artemis") + b_brand_Novum*(brand_1=="Novum")
+ b_country_CH*(country_1=="Switzerland") + b_country_DK*(country_1=="Denmark") + b_country_USA*(country_1=="USA")
+ b_char_standard*(char_1=="standard") + b_char_fast*(char_1=="fast acting") + b_char_double*(char_1=="double strength")
+ b_risk*side_effects_1
+ b_price*price_1 + EC_nest1)
V[["alt2"]] = ( b_brand_Artemis*(brand_2=="Artemis") + b_brand_Novum*(brand_2=="Novum")
+ b_country_CH*(country_2=="Switzerland") + b_country_DK*(country_2=="Denmark") + b_country_USA*(country_2=="USA")
+ b_char_standard*(char_2=="standard") + b_char_fast*(char_2=="fast acting") + b_char_double*(char_2=="double strength")
+ b_risk*side_effects_2
+ b_price*price_2 + EC_nest1)
V[["alt3"]] = ( b_brand_BestValue*(brand_3=="BestValue") + b_brand_Supermarket*(brand_3=="Supermarket") + b_brand_PainAway*(brand_3=="PainAway")
+ b_country_USA*(country_3=="USA") + b_country_IND*(country_3=="India") + b_country_RUS*(country_3=="Russia") + b_country_BRA*(country_3=="Brazil")
+ b_char_standard*(char_3=="standard") + b_char_fast*(char_3=="fast acting")
+ b_risk*side_effects_3
+ b_price*price_3)
V[["alt4"]] = ( b_brand_BestValue*(brand_4=="BestValue") + b_brand_Supermarket*(brand_4=="Supermarket") + b_brand_PainAway*(brand_4=="PainAway")
+ b_country_USA*(country_4=="USA") + b_country_IND*(country_4=="India") + b_country_RUS*(country_4=="Russia") + b_country_BRA*(country_4=="Brazil")
+ b_char_standard*(char_4=="standard") + b_char_fast*(char_4=="fast acting")
+ b_risk*side_effects_4
+ b_price*price_4)
### Define settings for exploded logit
el_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3, alt4=4),
avail = list(alt1=1, alt2=1, alt3=1, alt4=1),
choiceVars = list(best, second_pref, third_pref),
utilities = V,
scales = list(scale_1,scale_2,scale_3)
)
### Compute exploded logit probabilities
P[["model"]]=apollo_el(el_settings, functionality)
### Comment out as necessary
P = apollo_panelProd(P, apollo_inputs, functionality)
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION AND OUTPUT ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
apollo_modelOutput(model)
apollo_saveOutput(model)