I try to add membership attribute in Hybrid_LC_with_OL model example.
I just add two parameters and one variable to in membership.
Here is the code
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Initialise
rm(list = ls())
library(apollo)
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "Hybrid_LC_with_OL",
modelDescr = "Hybrid latent class choice model on drug choice data, using ordered measurement model for indicators",
indivID = "ID",
nCores = 5,
outputDirectory = "output"
)
# ################################################################# #
#### 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_swissRouteChoiceData
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(# Choice parameters
b_brand_Artemis_A = 0, b_brand_Artemis_B = 0,
b_brand_Novum_A = 0, b_brand_Novum_B = 0.1,
b_brand_BestValue_A = 0, b_brand_BestValue_B = 0.1,
b_brand_Supermarket_A = 0, b_brand_Supermarket_B = 0.1,
b_brand_PainAway_A = 0, b_brand_PainAway_B = 0.1,
b_country_CH_A = 0, b_country_CH_B = 0,
b_country_DK_A = 0, b_country_DK_B = 0,
b_country_USA_A = 0, b_country_USA_B = 0,
b_country_IND_A = 0, b_country_IND_B = 0,
b_country_RUS_A = 0, b_country_RUS_B = 0,
b_country_BRA_A = 0, b_country_BRA_B = 0,
b_char_standard_A = 0, b_char_standard_B = 0,
b_char_fast_A = 0, b_char_fast_B = 0,
b_char_double_A = 0, b_char_double_B = 0,
b_risk_A = 0, b_risk_B = 0,
b_price_A = 0, b_price_B = 0,
lambda_A = 1, lambda_B = 1,
# Class allocation parameters
delta_A = 0, delta_B = 0,
b_second_pref_a = 0,b_second_pref_b = 0,
# Measurement equations parameters
zeta_quality = 1, zeta_ingredient = 1, zeta_patent = 1, zeta_dominance = 1,
tau_quality_1 =-2, tau_quality_2 =-1, tau_quality_3 = 1, tau_quality_4 = 2,
tau_ingredients_1 =-2, tau_ingredients_2 =-1, tau_ingredients_3 = 1, tau_ingredients_4 = 2,
tau_patent_1 =-2, tau_patent_2 =-1, tau_patent_3 = 1, tau_patent_4 = 2,
tau_dominance_1 =-2, tau_dominance_2 =-1, tau_dominance_3 = 1, tau_dominance_4 = 2)
### 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_A", "b_country_USA_A", "b_char_standard_A", "delta_A","b_second_pref_a",
"b_brand_Artemis_B", "b_country_USA_B", "b_char_standard_B")
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 100,
interNormDraws = c("eta")
)
### Create random parameters
apollo_randCoeff=function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["LV"]] = eta
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_brand_Artemis" ]] = list(b_brand_Artemis_A , b_brand_Artemis_B )
lcpars[["b_brand_Novum" ]] = list(b_brand_Novum_A , b_brand_Novum_B )
lcpars[["b_brand_BestValue" ]] = list(b_brand_BestValue_A , b_brand_BestValue_B )
lcpars[["b_brand_Supermarket"]] = list(b_brand_Supermarket_A, b_brand_Supermarket_B)
lcpars[["b_brand_PainAway" ]] = list(b_brand_PainAway_A , b_brand_PainAway_B )
lcpars[["b_country_CH" ]] = list(b_country_CH_A , b_country_CH_B )
lcpars[["b_country_DK" ]] = list(b_country_DK_A , b_country_DK_B )
lcpars[["b_country_USA" ]] = list(b_country_USA_A , b_country_USA_B )
lcpars[["b_country_IND" ]] = list(b_country_IND_A , b_country_IND_B )
lcpars[["b_country_RUS" ]] = list(b_country_RUS_A , b_country_RUS_B )
lcpars[["b_country_BRA" ]] = list(b_country_BRA_A , b_country_BRA_B )
lcpars[["b_char_standard" ]] = list(b_char_standard_A , b_char_standard_B )
lcpars[["b_char_fast" ]] = list(b_char_fast_A , b_char_fast_B )
lcpars[["b_char_double" ]] = list(b_char_double_A , b_char_double_B )
lcpars[["b_risk" ]] = list(b_risk_A , b_risk_B )
lcpars[["b_price" ]] = list(b_price_A , b_price_B )
lcpars[["lambda" ]] = list(lambda_A , lambda_B )
##### Changed #######
V=list()
V[["class_a"]] = delta_A + b_second_pref_a * second_pref
V[["class_b"]] = delta_B + b_second_pref_b * second_pref
### Settings for class allocation models
classAlloc_settings = list(
classes = c(class_a=1, class_b=2),
utilities = V
)
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
return(lcpars)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
### Initialise
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()
### Likelihood of indicators
ol_settings1 = list(outcomeOrdered = attitude_quality,
V = zeta_quality*LV,
tau = list(tau_quality_1, tau_quality_2, tau_quality_3, tau_quality_4),
rows = (task==1))
ol_settings2 = list(outcomeOrdered = attitude_ingredients,
V = zeta_ingredient*LV,
tau = list(tau_ingredients_1, tau_ingredients_2, tau_ingredients_3, tau_ingredients_4),
rows = (task==1))
ol_settings3 = list(outcomeOrdered = attitude_patent,
V = zeta_patent*LV,
tau = list(tau_patent_1, tau_patent_2, tau_patent_3, tau_patent_4),
rows = (task==1))
ol_settings4 = list(outcomeOrdered = attitude_dominance,
V = zeta_dominance*LV,
tau = list(tau_dominance_1, tau_dominance_2, tau_dominance_3, tau_dominance_4),
rows = (task==1))
P[["indic_quality"]] = apollo_ol(ol_settings1, functionality)
P[["indic_ingredients"]] = apollo_ol(ol_settings2, functionality)
P[["indic_patent"]] = apollo_ol(ol_settings3, functionality)
P[["indic_dominance"]] = apollo_ol(ol_settings4, functionality)
P[["indic_quality"]] = apollo_panelProd(P[["indic_quality"]] , apollo_inputs, functionality)
P[["indic_ingredients"]] = apollo_panelProd(P[["indic_ingredients"]], apollo_inputs, functionality)
P[["indic_patent"]] = apollo_panelProd(P[["indic_patent"]] , apollo_inputs, functionality)
P[["indic_dominance"]] = apollo_panelProd(P[["indic_dominance"]] , apollo_inputs, functionality)
### Likelihood of choices inside each class
S <- 2
for(s in 1:S){
### Utilities for alternatives
V = list()
V[["alt1"]] = b_brand_Artemis[[s]] *(brand_1=="Artemis") + b_brand_Novum[[s]] *(brand_1=="Novum") + b_country_CH[[s]]* (country_1=="Switzerland") + b_country_DK[[s]] *(country_1=="Denmark") + b_country_USA[[s]]*(country_1=="USA") + b_char_standard[[s]]*(char_1=="standard") + b_char_fast[[s]]*(char_1=="fast acting") + b_char_double[[s]]*(char_1=="double strength") + b_risk[[s]]*side_effects_1 + b_price[[s]]*price_1 + lambda[[s]]*LV
V[["alt2"]] = b_brand_Artemis[[s]] *(brand_2=="Artemis") + b_brand_Novum[[s]] *(brand_2=="Novum") + b_country_CH[[s]]* (country_2=="Switzerland") + b_country_DK[[s]] *(country_2=="Denmark") + b_country_USA[[s]]*(country_2=="USA") + b_char_standard[[s]]*(char_2=="standard") + b_char_fast[[s]]*(char_2=="fast acting") + b_char_double[[s]]*(char_2=="double strength") + b_risk[[s]]*side_effects_2 + b_price[[s]]*price_2 + lambda[[s]]*LV
V[["alt3"]] = b_brand_BestValue[[s]]*(brand_3=="BestValue") + b_brand_Supermarket[[s]]*(brand_3=="Supermarket") + b_brand_PainAway[[s]]*(brand_3=="PainAway") + b_country_USA[[s]]*(country_3=="USA") + b_country_IND[[s]]*(country_3=="India") + b_country_RUS[[s]]*(country_3=="Russia") + b_country_BRA[[s]]*(country_3=="Brazil") + b_char_standard[[s]]*(char_3=="standard") + b_char_fast[[s]]*(char_3=="fast acting") + b_risk[[s]]*side_effects_3 + b_price[[s]]*price_3
V[["alt4"]] = b_brand_BestValue[[s]]*(brand_4=="BestValue") + b_brand_Supermarket[[s]]*(brand_4=="Supermarket") + b_brand_PainAway[[s]]*(brand_4=="PainAway") + b_country_USA[[s]]*(country_4=="USA") + b_country_IND[[s]]*(country_4=="India") + b_country_RUS[[s]]*(country_4=="Russia") + b_country_BRA[[s]]*(country_4=="Brazil") + b_char_standard[[s]]*(char_4=="standard") + b_char_fast[[s]]*(char_4=="fast acting") + b_risk[[s]]*side_effects_4 + b_price[[s]]*price_4
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3, alt4=4),
choiceVar = best,
utilities = V
)
### Compute within-class choice probabilities using MNL model
P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs ,functionality)
}
### Compute latent class model probabilities
lc_settings = list(inClassProb = P[paste0("Class_", 1:S)], classProb=pi_values)
P[["choice"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
### Comment out as necessary
P = apollo_combineModels(P, apollo_inputs, functionality)
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### CALCULATE LL AT STARTING VALUES ####
# ################################################################# #
apollo_llCalc(apollo_beta, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Best regards,
Xin Dou