I am estimating a latent class model using RP/SP data with both classical and Bayesian estimation methods. However, when estimating the model with Bayesian estimation, I am encountering NaN values for model fit indices such as AIC and BIC. I have even tried using starting values obtained from the classical estimation of the latent class model, but the issue persists.
I have also tried using (i) non-random (NR) priors for all the parameters, as well as (ii) negative lognormal (LN-) priors for alternative-specific variables and non-random priors for all other socio-economic variables. However, this did not help much.
Could you please suggest a way to resolve this?
The code is as follows:
Code: Select all
# ################################################################# #
#### 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 = "Apollo_IPA",
modelDescr = "RP_SP_HCM_3",
indivID = "ID",
HB = TRUE,
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
database = read.csv("RP_SP_LCCM_MXL.csv")
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc_local = 0,
asc_metro = 1,
asc_car = 0,
asc_tw = 0,
b_tt_a = -0.023225,
b_tc_a = -0.168560,
b_at_a = -0.028592,
b_cc_a = -0.106122,
b_tt_b = -0.061838,
b_tc_b = -0.014420,
b_at_b = 0.004892,
b_cc_b = -0.060425,
delta_a = 1,
delta_b = 0,
female_a = 1,
Income_0_50k_a = 1,
Income_above_100k_a = 1,
Age_18_30_a = 1,
Age_30_45_a = 1,
Education_PG_a = 1,
Children_yes_a = 1,
Sat_air_a = 1,
female_b = 0,
Income_0_50k_b = 1,
Income_above_100k_b = 0,
Age_18_30_b = 0,
Age_30_45_b = 0,
Education_PG_b = 0,
Children_yes_b = 0,
Sat_air_b = 0,
mu_RP = 1,
mu_SP = 1)
### 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("mu_RP",
"asc_metro",
"delta_a",
"female_a",
"Income_0_50k_a",
"Income_above_100k_a",
"Children_yes_a",
"Age_18_30_a",
"Age_30_45_a",
"Education_PG_a",
"Sat_air_a")
### Vector of parameters, including any that are kept fixed in estimation
apollo_HB = list(
hbDist = c(asc_local = "NR",
asc_metro = "NR",
asc_car = "NR",
asc_tw = "NR",
b_tt_a = "LN-",
b_tc_a = "LN-",
b_at_a = "LN-",
b_cc_a = "LN-",
b_tt_b = "LN-",
b_tc_b = "LN-",
b_at_b = "LN-",
b_cc_b = "LN-",
delta_a = "NR",
delta_b = "NR",
female_a = "NR",
Income_0_50k_a = "NR",
Income_above_100k_a = "NR",
Age_18_30_a = "NR",
Age_30_45_a = "NR",
Education_PG_a = "NR",
Children_yes_a = "NR",
Sat_air_a = "NR",
female_b = "NR",
Income_0_50k_b = "NR",
Income_above_100k_b = "NR",
Age_18_30_b = "NR",
Age_30_45_b = "NR",
Education_PG_b = "NR",
Children_yes_b = "NR",
Sat_air_b = "NR",
mu_RP = "NR",
mu_SP = "NR"),
gNCREP = 50000, # burn-in iterations
gNEREP = 20000, # post burn-in iterations
gINFOSKIP = 500,
gFULLCV = FALSE
)
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_tt"]] = list(b_tt_a, b_tt_b)
lcpars[["b_tc"]] = list(b_tc_a, b_tc_b)
lcpars[["b_at"]] = list(b_at_a, b_at_b)
lcpars[["b_cc"]] = list(b_cc_a, b_cc_b)
### Utilities of class allocation model
V=list()
V[["class_a"]] = delta_a + female_a*female + Income_0_50k_a*Income_0_50k + Income_above_100k_a*Income_above_100k + Age_18_30_a*Age_18_30 + Age_30_45_a*Age_30_45 +
Education_PG_a*Education_PG + Children_yes_a*Children_yes + Sat_air_a*Sat_air
V[["class_b"]] = delta_b + female_b*female + Income_0_50k_b*Income_0_50k + Income_above_100k_b*Income_above_100k + Age_18_30_b*Age_18_30 + Age_30_45_b*Age_30_45 +
Education_PG_b*Education_PG + Children_yes_b*Children_yes + Sat_air_b*Sat_air
### 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"){
### 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()
P_temp = list()
### Loop over classes
for(s in 1:2){
### List of utilities (before applying scales): these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[["local"]] = asc_local + b_tt[[s]]*tt_local + b_tc[[s]]*tc_local + b_at[[s]]*at_local
V[["metro"]] = asc_metro + b_tt[[s]]*tt_metro + b_tc[[s]]*tc_metro + b_at[[s]]*at_metro
V[["car"]] = asc_car + b_tt[[s]]*tt_car + b_tc[[s]]*tc_car + b_cc[[s]]*cc_car
V[["tw"]] = asc_tw + b_tt[[s]]*tt_tw + b_tc[[s]]*tc_tw + b_cc[[s]]*cc_tw
### Compute probabilities for the RP part of the data using MNL model
mnl_settings_RP = list(
alternatives = c(local=1, metro=2, car=3, tw=4),
avail = list(local=av_local, metro=av_metro, car=av_car, tw=av_tw),
choiceVar = Choice,
utilities = list(local = mu_RP*V[["local"]],
metro = mu_RP*V[["metro"]],
car = mu_RP*V[["car"]],
tw = mu_RP*V[["tw"]]),
rows = (RP==1)
)
P_temp[[paste0("RP_class_",s)]] = apollo_mnl(mnl_settings_RP, functionality)
### Compute probabilities for the RP part of the data using MNL model
mnl_settings_SP = list(
alternatives = c(local=1, metro=2, car=3, tw=4),
avail = list(local=av_local, metro=av_metro, car=av_car, tw=av_tw),
choiceVar = Choice,
utilities = list(local = mu_SP*V[["local"]],
metro = mu_SP*V[["metro"]],
car = mu_SP*V[["car"]],
tw = mu_SP*V[["tw"]]),
rows = (SP==1)
)
P_temp[[paste0("SP_class_",s)]] = apollo_mnl(mnl_settings_SP, functionality)
### Combined model
P[[paste0("Class_",s)]] = apollo_combineModels(P_temp, apollo_inputs, functionality,components=c(paste0("RP_class_",s),paste0("SP_class_",s)),asList=FALSE)
}
### Compute latent class model probabilities
lc_settings = list(inClassProb = P, classProb=pi_values)
P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
### Estimate model
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
### Show output in screen
apollo_modelOutput(model)
### Save output to file(s)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- LATENT CLASSES ----
# ----------------------------------------------------------------- #
### Compute unconditionals
unconditionals = apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
### Compute conditionals
conditionals = apollo_conditionals(model,apollo_probabilities, apollo_inputs)
summary(conditionals)
summary(as.data.frame(unconditionals[["pi_values"]]))The estimation results are as follows:
Code: Select all
Model name : Apollo_IPA
Model description : RP_SP_HCM_3
Model run at : 2024-10-05 13:41:34.66498
Estimation method : Hierarchical Bayes
Number of individuals : 740
Number of rows in database : 7400
Number of modelled outcomes : 29600
model : 14800
RP_class_1 : 740
SP_class_1 : 6660
RP_class_1 : 740
SP_class_1 : 6660
RP_class_2 : 740
SP_class_2 : 6660
Number of cores used : 1
Estimation carried out using RSGHB
Burn-in iterations : 50000
Post burn-in iterations : 20000
Classical model fit statistics were calculated at parameter values obtained using averaging across
the post burn-in iterations.
LL(start) : -85188.35
LL (whole model) at equal shares, LL(0) : -8485.5
LL (whole model) at observed shares, LL(C) : -7690.5
LL(final, whole model) : NaN
Rho-squared vs equal shares : NaN
Adj.Rho-squared vs equal shares : NaN
Rho-squared vs observed shares : NaN
Adj.Rho-squared vs observed shares : NaN
AIC : NaN
BIC : NaN
Equiv. estimated parameters : 57
(non-random parameters : 13)
(means of random parameters : 8)
(covariance matrix terms : 36)
LL(0,Class_1) : -8485.5
LL(final,Class_1) : -5141.92
LL(0,Class_2) : -8485.5
LL(final,Class_2) : NaN Likelihood equal to zero for at least
one individual in this component.
Time taken (hh:mm:ss) : 00:28:34.94
pre-estimation : 00:00:1.46
estimation : 00:28:26.79
post-estimation : 00:00:6.69
Summary of parameter chains
Non-random coefficients
Mean SD
asc_local 0.6879 0.0481
asc_metro 1.0000 NA
asc_car 5.1753 0.3206
asc_tw 2.2240 0.1086
delta_a 1.0000 NA
delta_b 4.2102 0.5932
female_a 1.0000 NA
Income_0_50k_a 1.0000 NA
Income_above_100k_a 1.0000 NA
Age_18_30_a 1.0000 NA
Age_30_45_a 1.0000 NA
Education_PG_a 1.0000 NA
Children_yes_a 1.0000 NA
Sat_air_a 1.0000 NA
female_b 2.7746 0.6828
Income_0_50k_b -1.7494 1.0562
Income_above_100k_b 3.0297 0.7986
Age_18_30_b -6.2924 1.3248
Age_30_45_b -2.1368 0.8120
Education_PG_b -12.1400 1.7952
Children_yes_b -5.6329 0.9288
Sat_air_b -4.3782 1.2769
mu_RP 1.0000 NA
mu_SP 1.2713 0.0996
Results for posterior means for random coefficients
Mean SD
b_tt_a -0.0754 0.0384
b_tc_a -0.1085 0.1274
b_at_a -0.0034 0.0014
b_cc_a -0.1360 0.0706
b_tt_b -0.2428 0.0193
b_tc_b -1.7398 0.1866
b_at_b -1.2562 0.1841
b_cc_b -0.2188 0.0596
Summary of distributions of random coeffients (after distributional transforms)
Mean SD
b_tt_a -0.0754 0.0659
b_tc_a -0.1044 0.2583
b_at_a -0.0033 0.0034
b_cc_a -0.1355 0.1175
b_tt_b -0.2291 0.1472
b_tc_b -1.6998 0.8322
b_at_b -1.1643 0.7054
b_cc_b -0.2054 0.1653
Upper level model results for mean parameters for underlying Normals
Mean SD
b_tt_a -2.8771 0.0872
b_tc_a -3.3289 0.1105
b_at_a -6.0654 0.3263
b_cc_a -2.2840 0.1218
b_tt_b -1.6424 0.3778
b_tc_b 0.4277 0.2013
b_at_b -0.0024 0.4088
b_cc_b -1.8354 0.3716