Could the LL of a full model worse than that for the reduced model?
Posted: 01 Jul 2020, 12:31
Dear Stephane and David,
Thank you very much for inventing the Apollo package and creating the forum allowing our users to communicate.
I recently come across a question in choice modelling using the Apollo package. The setting of my DCE scenario is: there is a health attribute and it is described as either an improvement or a deterioration from the current health situation. The model setting is: I have an asymmetric model in which one parameter (i.e., health) is split into two according to whether the attribute is described as an improvement (i.e., health_increase) or a deterioration (i.e., health_decrease) from a reference point.The idea is similar to Hess et al. (2008). The reduced model is the symmetric model where the health parameter is specified as linear regardless of its relative position to the reference point.
Theoretically, the two models should be nested, as the asymmetric model can be restricted to the symmetric model by imposing a constraint of beta[health_increase]=-beta[health_decrease]. It is expected that the LL (likelihood) for the asymmetric model is better than that for the symmetric model, when both models are specified as basic MNL. I did empirically find this using the MNL.
HOWEVER, I expected that this is also true (i.e., the LL for the asymmetric model is better than that for the symmetric model) when two models are specified as RPL where health parameters follow normal distributions. but the result suggests that the LL for the asymmetric model (in RPL form) is WORSE . This is a bit of a surprise, because given the two models are nested, I thought the model fit of the asymmetric model should be better anyway (although may not be significantly better).
I don't know if it is because the two models are not nested anymore when they are specified as RPL (e.g., because draws for the linear health parameter in the symmetric model and for the health_increase and health_decrease parameters in the asymmetric model are from different independent normal distributions) ? or because some estimation issues (e.g., one of the models reaches a local maximum or I should try different distributions?). Codes are attached. Thank you very much.
Reference
Hess, S., Rose, J. M., & Hensher, D. A. (2008). Asymmetric preference formation in willingness to pay estimates in discrete choice models. Transportation Research Part E: Logistics and Transportation Review, 44(5), 847-863.
Codes
Asymmetric model
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="RPL_asymmetric",
modelDescr ="RPL_asymmetric",
mixing = TRUE,
indivID ="ID",
nCores = 3
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("R_T3.csv",header=TRUE)
##renames##
names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$same_answersG==0&same_answersQ==0)
#
names(database)[names(database)=="E_h_1_inc"]<-"h_1_inc"
names(database)[names(database)=="E_h_2_inc"]<-"h_2_inc"
names(database)[names(database)=="E_h_1_dec"]<-"h_1_dec"
names(database)[names(database)=="E_h_2_dec"]<-"h_2_dec"
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1 = 0,
asc_3_mu = 0,
asc_3_sig = 0,
h_inc_mu = 0,
h_inc_sig = 0,
h_dec_mu = 0,
h_dec_sig = 0,
v_mu = 0,
v_sig = 0,
c = 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_1") #(for cost_ANA: don't set delta the same as ANA class)
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="mlhs",
interNDraws=500,
interUnifDraws=c(),
interNormDraws=c("draws_v",
"draws_h_inc","draws_h_dec","draws_asc_3"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["hinc"]] = h_inc_mu + h_inc_sig*draws_h_inc
randcoeff[["hdec"]] = h_dec_mu + h_dec_sig*draws_h_dec
randcoeff[["v"]] = v_mu + v_sig*draws_v
randcoeff[["asc_3"]] = asc_3_mu + asc_3_sig*draws_asc_3
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()
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_1 + hinc*h_1_inc + hdec*h_1_dec + v*v_1+ c*c_1
V[['alt2']] = asc_1 + hinc*h_2_inc + hdec*h_2_dec + v*v_2 + c*c_2
V[['alt3']] = asc_3 + v*v_3 + c*c_3
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice,
V = 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)
### 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 ####
# ################################################################# #
### Optional starting values search
# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Symmetric model
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="RPL_symmetric",
modelDescr ="RPL_symmetric",
mixing = TRUE,
indivID ="ID",
nCores = 3
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("R_T3.csv",header=TRUE)
##renames##
names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$same_answersG==0&same_answersQ==0)
#
names(database)[names(database)=="E_h_1_inc"]<-"h_1_inc"
names(database)[names(database)=="E_h_2_inc"]<-"h_2_inc"
names(database)[names(database)=="E_h_1_dec"]<-"h_1_dec"
names(database)[names(database)=="E_h_2_dec"]<-"h_2_dec"
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1 = 0,
asc_3_mu = 0,
asc_3_sig = 0,
h = 0,
v_mu = 0,
v_sig = 0,
c = 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_1") #(for cost_ANA: don't set delta the same as ANA class)
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="mlhs",
interNDraws=500,
interUnifDraws=c(),
interNormDraws=c("draws_v",
"draws_h","draws_asc_3"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["h"]] = h_mu + h_sig*draws_h
randcoeff[["v"]] = v_mu + v_sig*draws_v
randcoeff[["asc_3"]] = asc_3_mu + asc_3_sig*draws_asc_3
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()
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_1 + h*h_1 + v*v_1 + c*c_1
V[['alt2']] = asc_1 + h*h_2 + v*v_2 + c*c_2
V[['alt3']] = asc_3 + h*h_3 + v*v_3 + c*c_3
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice,
V = 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)
### 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 ####
# ################################################################# #
### Optional starting values search
# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
----------
Tim
Thank you very much for inventing the Apollo package and creating the forum allowing our users to communicate.
I recently come across a question in choice modelling using the Apollo package. The setting of my DCE scenario is: there is a health attribute and it is described as either an improvement or a deterioration from the current health situation. The model setting is: I have an asymmetric model in which one parameter (i.e., health) is split into two according to whether the attribute is described as an improvement (i.e., health_increase) or a deterioration (i.e., health_decrease) from a reference point.The idea is similar to Hess et al. (2008). The reduced model is the symmetric model where the health parameter is specified as linear regardless of its relative position to the reference point.
Theoretically, the two models should be nested, as the asymmetric model can be restricted to the symmetric model by imposing a constraint of beta[health_increase]=-beta[health_decrease]. It is expected that the LL (likelihood) for the asymmetric model is better than that for the symmetric model, when both models are specified as basic MNL. I did empirically find this using the MNL.
HOWEVER, I expected that this is also true (i.e., the LL for the asymmetric model is better than that for the symmetric model) when two models are specified as RPL where health parameters follow normal distributions. but the result suggests that the LL for the asymmetric model (in RPL form) is WORSE . This is a bit of a surprise, because given the two models are nested, I thought the model fit of the asymmetric model should be better anyway (although may not be significantly better).
I don't know if it is because the two models are not nested anymore when they are specified as RPL (e.g., because draws for the linear health parameter in the symmetric model and for the health_increase and health_decrease parameters in the asymmetric model are from different independent normal distributions) ? or because some estimation issues (e.g., one of the models reaches a local maximum or I should try different distributions?). Codes are attached. Thank you very much.
Reference
Hess, S., Rose, J. M., & Hensher, D. A. (2008). Asymmetric preference formation in willingness to pay estimates in discrete choice models. Transportation Research Part E: Logistics and Transportation Review, 44(5), 847-863.
Codes
Asymmetric model
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="RPL_asymmetric",
modelDescr ="RPL_asymmetric",
mixing = TRUE,
indivID ="ID",
nCores = 3
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("R_T3.csv",header=TRUE)
##renames##
names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$same_answersG==0&same_answersQ==0)
#
names(database)[names(database)=="E_h_1_inc"]<-"h_1_inc"
names(database)[names(database)=="E_h_2_inc"]<-"h_2_inc"
names(database)[names(database)=="E_h_1_dec"]<-"h_1_dec"
names(database)[names(database)=="E_h_2_dec"]<-"h_2_dec"
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1 = 0,
asc_3_mu = 0,
asc_3_sig = 0,
h_inc_mu = 0,
h_inc_sig = 0,
h_dec_mu = 0,
h_dec_sig = 0,
v_mu = 0,
v_sig = 0,
c = 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_1") #(for cost_ANA: don't set delta the same as ANA class)
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="mlhs",
interNDraws=500,
interUnifDraws=c(),
interNormDraws=c("draws_v",
"draws_h_inc","draws_h_dec","draws_asc_3"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["hinc"]] = h_inc_mu + h_inc_sig*draws_h_inc
randcoeff[["hdec"]] = h_dec_mu + h_dec_sig*draws_h_dec
randcoeff[["v"]] = v_mu + v_sig*draws_v
randcoeff[["asc_3"]] = asc_3_mu + asc_3_sig*draws_asc_3
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()
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_1 + hinc*h_1_inc + hdec*h_1_dec + v*v_1+ c*c_1
V[['alt2']] = asc_1 + hinc*h_2_inc + hdec*h_2_dec + v*v_2 + c*c_2
V[['alt3']] = asc_3 + v*v_3 + c*c_3
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice,
V = 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)
### 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 ####
# ################################################################# #
### Optional starting values search
# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Symmetric model
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="RPL_symmetric",
modelDescr ="RPL_symmetric",
mixing = TRUE,
indivID ="ID",
nCores = 3
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("R_T3.csv",header=TRUE)
##renames##
names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$same_answersG==0&same_answersQ==0)
#
names(database)[names(database)=="E_h_1_inc"]<-"h_1_inc"
names(database)[names(database)=="E_h_2_inc"]<-"h_2_inc"
names(database)[names(database)=="E_h_1_dec"]<-"h_1_dec"
names(database)[names(database)=="E_h_2_dec"]<-"h_2_dec"
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1 = 0,
asc_3_mu = 0,
asc_3_sig = 0,
h = 0,
v_mu = 0,
v_sig = 0,
c = 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_1") #(for cost_ANA: don't set delta the same as ANA class)
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="mlhs",
interNDraws=500,
interUnifDraws=c(),
interNormDraws=c("draws_v",
"draws_h","draws_asc_3"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["h"]] = h_mu + h_sig*draws_h
randcoeff[["v"]] = v_mu + v_sig*draws_v
randcoeff[["asc_3"]] = asc_3_mu + asc_3_sig*draws_asc_3
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()
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_1 + h*h_1 + v*v_1 + c*c_1
V[['alt2']] = asc_1 + h*h_2 + v*v_2 + c*c_2
V[['alt3']] = asc_3 + h*h_3 + v*v_3 + c*c_3
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice,
V = 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)
### 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 ####
# ################################################################# #
### Optional starting values search
# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
----------
Tim