I have a simple unlabelled experiment in which respondents choose between two drugs: A or B, depending on 5 variables:
[QOL] - a categorical scenario variable that is included in the model as an interaction (dummy coded)
[LIFE] - expectancy - a continuous scenario variable also included in the model as an interaction
The level of [BENEFIT] the drug provides - a categorical variable (dummy coded)
How [CERT]ain we are the drug works - a categorical variable (dummy coded)
[WAIT] time until available - a continuous variable.
QOL, LIFE and BENEFIT are constant across alternatives (to ease cognitive burden for respondents) but these do vary between choice tasks. Because of this, there is no trade-off data between these attributes to produce main effect beta's. Instead, I have included interactions of these with CERT and WAIT and plan to analyse how they influence sensitivity to of the two other attributes. However, because 3/5 of my attributes are dummy-coded categorical variables, the available interactions has quickly become burdensome as there are so many to choose from, particularly for generating priors for my Bayesian design with only 10% (n=80) of my sample.
My questions:
1)
- Instead of interacting all categorical levels of attributes together which quickly becomes cumbersome (and an over-specification), other than using common sense, what is the best way to prioritise which to include in the model?
- For example, CERT has 4 categorical levels, and QOL has 4 categorical levels. When interacting them, I have interacted the lowest level of CERT as a ref (keeping this lowest level "fixed" in Apollo_fixed) with the highest level of CERT for comparison.
-Using this strategy, you can see the categorical levels I have left out as I have marked them with a # in Apollo_beta, Apollo_fixed and in the Utility functions. To me this makes sense, as in the output, I have at least two levels to compare agains (one being kept as a reference category). Is this use of a reference category for interactions correct?
2)
As I have so many interactions, would it be easier to "Create coefficients" for interactions of scenario variables (QOL, LIFE and BENEFIT). I have noticed many do this before the utility functions for sociodemographics? If you think this would be beneficial, would you be able to provide an example of how to do this for one of my interaction coefficients?
3)
I would value any other feedback / guidance on my code if there is anything that you would do differently before using it to generate priors for a bayesian efficient design.
n.b. the betas here are from bots (simulation) as I don't have any data yet, but this is the model I plan to use to obtain Bayesian Priors.
Code:
Code: Select all
memory
rm(list = ls())
### Load libraries
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "Pilot_v18_SimulationData_Randomised_05_July_2023",
modelDescr = "Simple MNL model with categorical data using example (simulation) (n = ~90)",
indivID = "RID",
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################ #
database = read.csv("pilot_v4_18_SimulationData_Randomised.csv",header=TRUE)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(
#asc:
asc_DrugA = 0,
#main effects:
b_CERT_very_low = 0,
b_CERT_low = 0,
b_CERT_moderate = 0,
b_CERT_high = 0,
b_WAIT = 0,
#Do not include the following QOL,
#LIFE and BENEFIT on their
#own as there is no trade-off
#data so cannot be estimated (constant across alternatives)
# b_QOL_light_work = 0,
# b_QOL_no_work = 0,
# b_QOL_bed_50 = 0,
# b_QOL_comp_disabled = 0,
# b_LIFE = 0,
# b_BENEFIT_small = 0,
# b_BENEFIT_substantial = 0,
# b_BENEFIT_moderate = 0,
#interaction effects
#(with WAIT) for
#scenario variables:
b_LIFExWAIT = 0,
b_QOL_light_workxWAIT = 0,
b_QOL_no_workxWAIT = 0,
b_QOL_bed_50xWAIT = 0,
b_QOL_comp_disabledxWAIT = 0,
b_BENEFIT_smallxWAIT = 0,
b_BENEFIT_moderatexWAIT = 0,
b_BENEFIT_substantialxWAIT= 0,
#interaction effects
#(with CERT) for
#scenario variables:
b_LIFExCERT_very_low = 0,
b_LIFExCERT_low = 0,
b_LIFExCERT_moderate = 0,
b_LIFExCERT_high = 0,
b_QOL_light_workxCERT_very_low = 0,
b_QOL_light_workxCERT_low = 0,
b_QOL_light_workxCERT_moderate = 0,
b_QOL_light_workxCERT_high = 0,
#b_QOL_no_workxCERT_very_low = 0,
#b_QOL_no_workxCERT_low = 0,
#b_QOL_no_workxCERT_moderate = 0,
#b_QOL_no_workxCERT_high = 0,
b_QOL_bed_50xCERT_very_low = 0,
b_QOL_bed_50xCERT_low = 0,
b_QOL_bed_50xCERT_moderate = 0,
b_QOL_bed_50xCERT_high = 0,
b_QOL_comp_disabledxCERT_very_low = 0,
b_QOL_comp_disabledxCERT_low = 0,
b_QOL_comp_disabledxCERT_moderate = 0,
b_QOL_comp_disabledxCERT_high = 0,
b_BENEFIT_smallxCERT_very_low = 0,
b_BENEFIT_smallxCERT_low = 0,
b_BENEFIT_smallxCERT_moderate = 0,
b_BENEFIT_smallxCERT_high = 0,
#b_BENEFIT_moderatexCERT_very_low = 0,
#b_BENEFIT_moderatexCERT_low = 0,
#b_BENEFIT_moderatexCERT_moderate = 0,
#b_BENEFIT_moderatexCERT_high = 0,
b_BENEFIT_substantialxCERT_very_low = 0,
b_BENEFIT_substantialxCERT_low = 0,
b_BENEFIT_substantialxCERT_moderate = 0,
b_BENEFIT_substantialxCERT_high = 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(
"asc_DrugA",
"b_CERT_very_low",
"b_QOL_light_workxWAIT",
"b_BENEFIT_smallxWAIT",
"b_LIFExCERT_very_low",
"b_QOL_light_workxCERT_very_low",
#"b_QOL_no_workxCERT_very_low",
"b_QOL_bed_50xCERT_very_low",
"b_QOL_comp_disabledxCERT_very_low",
"b_BENEFIT_smallxCERT_very_low",
#"b_BENEFIT_moderatexCERT_very_low",
"b_BENEFIT_substantialxCERT_very_low"
)
# ################################################################# #
#### 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()
### Create coefficients for interactions of scenario variables (QOL, LIFE and BENEFIT) ###
#
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[["DrugA"]] =
asc_DrugA +
b_CERT_very_low * (DrugA_CERT == 1) +
b_CERT_low * (DrugA_CERT == 2) +
b_CERT_moderate * (DrugA_CERT == 3) +
b_CERT_high * (DrugA_CERT == 4) +
b_WAIT * DrugA_WAIT +
#Should there be brackets around these interaction effects? e.g.:
# (b_LIFExWAIT * DrugA_LIFE) * DrugA_WAIT +
#Interactions with WAIT:
b_LIFExWAIT * DrugA_LIFE * DrugA_WAIT +
b_QOL_light_workxWAIT * ( DrugA_WAIT * (DrugA_QOL == 1)) +
b_QOL_no_workxWAIT * ( DrugA_WAIT * (DrugA_QOL == 2)) +
b_QOL_bed_50xWAIT * ( DrugA_WAIT * (DrugA_QOL == 3)) +
b_QOL_comp_disabledxWAIT * ( DrugA_WAIT * (DrugA_QOL == 4)) +
b_BENEFIT_smallxWAIT * ( DrugA_WAIT * (DrugA_BENEFIT == 1)) +
b_BENEFIT_moderatexWAIT * ( DrugA_WAIT * (DrugA_BENEFIT == 2)) +
b_BENEFIT_substantialxWAIT * ( DrugA_WAIT * (DrugA_BENEFIT == 3)) +
#Interactions with CERT (all categorical levels):
b_LIFExCERT_very_low * ( DrugA_CERT == 1) * DrugA_LIFE +
b_LIFExCERT_low * ( DrugA_CERT == 2) * DrugA_LIFE +
b_LIFExCERT_moderate * ( DrugA_CERT == 3) * DrugA_LIFE +
b_LIFExCERT_high * ( DrugA_CERT == 4) * DrugA_LIFE +
b_QOL_light_workxCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_QOL == 1)) +
b_QOL_light_workxCERT_low * ( (DrugA_CERT == 2) * (DrugA_QOL == 1)) +
b_QOL_light_workxCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_QOL == 1)) +
b_QOL_light_workxCERT_high * ( (DrugA_CERT == 4) * (DrugA_QOL == 1)) +
# b_QOL_no_workxCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_QOL == 2)) +
# b_QOL_no_workxCERT_low * ( (DrugA_CERT == 2) * (DrugA_QOL == 2)) +
# b_QOL_no_workxCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_QOL == 2)) +
# b_QOL_no_workxCERT_high * ( (DrugA_CERT == 4) * (DrugA_QOL == 2)) +
b_QOL_bed_50xCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_QOL == 3)) +
b_QOL_bed_50xCERT_low * ( (DrugA_CERT == 2) * (DrugA_QOL == 3)) +
b_QOL_bed_50xCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_QOL == 3)) +
b_QOL_bed_50xCERT_high * ( (DrugA_CERT == 4) * (DrugA_QOL == 3)) +
b_QOL_comp_disabledxCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_QOL == 4)) +
b_QOL_comp_disabledxCERT_low * ( (DrugA_CERT == 2) * (DrugA_QOL == 4)) +
b_QOL_comp_disabledxCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_QOL == 4)) +
b_QOL_comp_disabledxCERT_high * ( (DrugA_CERT == 4) * (DrugA_QOL == 4)) +
b_BENEFIT_smallxCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_low * ( (DrugA_CERT == 2) * (DrugA_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_high * ( (DrugA_CERT == 4) * (DrugA_BENEFIT == 1)) +
# b_BENEFIT_moderatexCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_low * ( (DrugA_CERT == 2) * (DrugA_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_high * ( (DrugA_CERT == 4) * (DrugA_BENEFIT == 2)) +
b_BENEFIT_substantialxCERT_very_low * ( (DrugA_CERT == 1) * (DrugA_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_low * ( (DrugA_CERT == 2) * (DrugA_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_moderate * ( (DrugA_CERT == 3) * (DrugA_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_high * ( (DrugA_CERT == 4) * (DrugA_BENEFIT == 3))
V[["DrugB"]] =
b_CERT_very_low * (DrugB_CERT == 1) +
b_CERT_low * (DrugB_CERT == 2) +
b_CERT_moderate * (DrugB_CERT == 3) +
b_CERT_high * (DrugB_CERT == 4) +
b_WAIT * DrugB_WAIT +
#Interactions with WAIT:
b_LIFExWAIT * DrugB_WAIT * DrugB_LIFE +
b_QOL_light_workxWAIT * ( DrugB_WAIT * (DrugB_QOL == 1)) +
b_QOL_no_workxWAIT * ( DrugB_WAIT * (DrugB_QOL == 2)) +
b_QOL_bed_50xWAIT * ( DrugB_WAIT * (DrugB_QOL == 3)) +
b_QOL_comp_disabledxWAIT * ( DrugB_WAIT * (DrugB_QOL == 4)) +
b_BENEFIT_smallxWAIT * ( DrugB_WAIT * (DrugB_BENEFIT == 1)) +
b_BENEFIT_moderatexWAIT * ( DrugB_WAIT * (DrugB_BENEFIT == 2)) +
b_BENEFIT_substantialxWAIT * ( DrugB_WAIT * (DrugB_BENEFIT == 3)) +
#Interactions with CERT (all categorical levels):
b_LIFExCERT_very_low * (DrugB_CERT == 1) * DrugB_LIFE +
b_LIFExCERT_low * (DrugB_CERT == 2) * DrugB_LIFE +
b_LIFExCERT_moderate * (DrugB_CERT == 3) * DrugB_LIFE +
b_LIFExCERT_high * (DrugB_CERT == 4) * DrugB_LIFE +
b_QOL_light_workxCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_QOL == 1)) +
b_QOL_light_workxCERT_low * ( (DrugB_CERT == 2) * (DrugB_QOL == 1)) +
b_QOL_light_workxCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_QOL == 1)) +
b_QOL_light_workxCERT_high * ( (DrugB_CERT == 4) * (DrugB_QOL == 1)) +
# b_QOL_no_workxCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_QOL == 2)) +
# b_QOL_no_workxCERT_low * ( (DrugB_CERT == 2) * (DrugB_QOL == 2)) +
# b_QOL_no_workxCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_QOL == 2)) +
# b_QOL_no_workxCERT_high * ( (DrugB_CERT == 4) * (DrugB_QOL == 2)) +
b_QOL_bed_50xCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_QOL == 3)) +
b_QOL_bed_50xCERT_low * ( (DrugB_CERT == 2) * (DrugB_QOL == 3)) +
b_QOL_bed_50xCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_QOL == 3)) +
b_QOL_bed_50xCERT_high * ( (DrugB_CERT == 4) * (DrugB_QOL == 3)) +
b_QOL_comp_disabledxCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_QOL == 4)) +
b_QOL_comp_disabledxCERT_low * ( (DrugB_CERT == 2) * (DrugB_QOL == 4)) +
b_QOL_comp_disabledxCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_QOL == 4)) +
b_QOL_comp_disabledxCERT_high * ( (DrugB_CERT == 4) * (DrugB_QOL == 4)) +
b_BENEFIT_smallxCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_low * ( (DrugB_CERT == 2) * (DrugB_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_BENEFIT == 1)) +
b_BENEFIT_smallxCERT_high * ( (DrugB_CERT == 4) * (DrugB_BENEFIT == 1)) +
# b_BENEFIT_moderatexCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_low * ( (DrugB_CERT == 2) * (DrugB_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_BENEFIT == 2)) +
# b_BENEFIT_moderatexCERT_high * ( (DrugB_CERT == 4) * (DrugB_BENEFIT == 2)) +
b_BENEFIT_substantialxCERT_very_low * ( (DrugB_CERT == 1) * (DrugB_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_low * ( (DrugB_CERT == 2) * (DrugB_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_moderate * ( (DrugB_CERT == 3) * (DrugB_BENEFIT == 3)) +
b_BENEFIT_substantialxCERT_high * ( (DrugB_CERT == 4) * (DrugB_BENEFIT == 3))
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(DrugA=1, DrugB=2),
choiceVar = CHOICE,
utilities = V
)
### 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)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### What does apollo_probabilities return? ####
# ################################################################# #
#This shows you the prob of that choice being made, given a number of alternatives and choice tasks
#For example, prob of A chosen 11 time = 0.5^11 = 0.00048828125
#Not needed - just for illustration purposes
apollo_probabilities(apollo_beta, apollo_inputs, functionality="estimate")
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Output:
Is my use of a reference category for interactions suitable?
Code: Select all
Model name : Pilot_v18_SimulationData_Randomised_05_July_2023
Model description : Simple MNL model with categorical data using example (simulation) (n = ~90)
Model run at : 2023-07-05 13:07:40.812633
Estimation method : bfgs
Model diagnosis : successful convergence
Optimisation diagnosis : Maximum found
hessian properties : Negative definitive
maximum eigenvalue : -0.672076
Number of individuals : 113
Number of rows in database : 1243
Number of modelled outcomes : 1243
Number of cores used : 1
Model without mixing
LL(start) : -861.58
LL at equal shares, LL(0) : -861.58
LL at observed shares, LL(C) : -860.77
LL(final) : -847.46
Rho-squared vs equal shares : 0.0164
Adj.Rho-squared vs equal shares : -0.0161
Rho-squared vs observed shares : 0.0155
Adj.Rho-squared vs observed shares : -0.0159
AIC : 1750.91
BIC : 1894.42
Estimated parameters : 28
Time taken (hh:mm:ss) : 00:00:1.88
pre-estimation : 00:00:0.97
estimation : 00:00:0.23
initial estimation : 00:00:0.21
estimation after rescaling : 00:00:0.01
post-estimation : 00:00:0.68
Iterations : 36
initial estimation : 35
estimation after rescaling : 1
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_DrugA 0.00000 NA NA NA NA
b_CERT_very_low 0.00000 NA NA NA NA
b_CERT_low 1.11608 0.69570 1.6043 0.66063 1.6894
b_CERT_moderate -0.38250 0.74456 -0.5137 0.75688 -0.5054
b_CERT_high -0.48795 0.72078 -0.6770 0.68442 -0.7129
b_WAIT 0.13674 0.36267 0.3770 0.36367 0.3760
b_LIFExWAIT -0.04542 0.08890 -0.5109 0.08602 -0.5281
b_QOL_light_workxWAIT 0.00000 NA NA NA NA
b_QOL_no_workxWAIT -0.04045 0.10128 -0.3994 0.09384 -0.4310
b_QOL_bed_50xWAIT 0.03902 0.10375 0.3761 0.09892 0.3945
b_QOL_comp_disabledxWAIT -0.02047 0.10155 -0.2016 0.09242 -0.2215
b_BENEFIT_smallxWAIT 0.00000 NA NA NA NA
b_BENEFIT_moderatexWAIT -0.04383 0.11060 -0.3963 0.11628 -0.3769
b_BENEFIT_substantialxWAIT -0.05969 0.09432 -0.6329 0.09449 -0.6317
b_LIFExCERT_very_low 0.00000 NA NA NA NA
b_LIFExCERT_low -0.10255 0.17848 -0.5746 0.16100 -0.6370
b_LIFExCERT_moderate 0.18427 0.19081 0.9657 0.19251 0.9572
b_LIFExCERT_high 0.33144 0.18911 1.7526 0.17597 1.8835
b_QOL_light_workxCERT_very_low 0.00000 NA NA NA NA
b_QOL_light_workxCERT_low -0.53769 0.32843 -1.6372 0.34599 -1.5541
b_QOL_light_workxCERT_moderate -0.22553 0.41085 -0.5489 0.41772 -0.5399
b_QOL_light_workxCERT_high -0.41790 0.44510 -0.9389 0.44234 -0.9447
b_QOL_bed_50xCERT_very_low 0.00000 NA NA NA NA
b_QOL_bed_50xCERT_low -0.52223 0.35169 -1.4849 0.37476 -1.3935
b_QOL_bed_50xCERT_moderate -0.15698 0.44449 -0.3532 0.49876 -0.3147
b_QOL_bed_50xCERT_high -0.50979 0.55683 -0.9155 0.56118 -0.9084
b_QOL_comp_disabledxCERT_very_low 0.00000 NA NA NA NA
b_QOL_comp_disabledxCERT_low -0.34493 0.34869 -0.9892 0.36309 -0.9500
b_QOL_comp_disabledxCERT_moderate -0.13042 0.46721 -0.2791 0.50525 -0.2581
b_QOL_comp_disabledxCERT_high -0.39992 0.55313 -0.7230 0.57989 -0.6896
b_BENEFIT_smallxCERT_very_low 0.00000 NA NA NA NA
b_BENEFIT_smallxCERT_low -0.54424 0.33151 -1.6417 0.33711 -1.6144
b_BENEFIT_smallxCERT_moderate -0.39646 0.50499 -0.7851 0.51274 -0.7732
b_BENEFIT_smallxCERT_high -0.49908 0.65097 -0.7667 0.68440 -0.7292
b_BENEFIT_substantialxCERT_very_low 0.00000 NA NA NA NA
b_BENEFIT_substantialxCERT_low -0.26905 0.29143 -0.9232 0.29456 -0.9134
b_BENEFIT_substantialxCERT_moderate 0.19418 0.38620 0.5028 0.38873 0.4995
b_BENEFIT_substantialxCERT_high 0.19482 0.43793 0.4449 0.43326 0.4497