ANA_EM has no covariance matrix

[kiki]

Hi,

I use the EM algorithm to model the LC_ANA model.

I try the following ways to adjust the model. All of them present the covariance matrix is not available. So maybe it is not the problem about the model.

1) Consider the parameters in class membership model according to attributes
2) Consider the parameters in class membership model according to class
3) Consider the parameters in class membership model according to attributes and class
4) Remove one ASC (i.e., J-1).
5) Not remove one ASC.
6) Remove one delta_d (constant in class membership)
7) Not remove delta_d

The codes and the results are as follows. Thank you for your check in advance!

Best regards,
kiki
-------------------------------------------------------------------------------------------------------------------------------------------------
# ################################################################# #
#### 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 = "EM_LC_with_covariates_ANA",
modelDescr = "LC model with class allocation model on Swiss route choice data, EM algorithm, ANA",
indivID = "ID",
noValidation = TRUE,
noDiagnostics = TRUE,
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_swissRouteChoiceData
### for data dictionary, use ?apollo_swissRouteChoiceData

### create weights column for use in EM algorithm
database$weights = 1

# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #


# Create the matrix
Q <- matrix(c(0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L), nrow = 2, byrow = TRUE)

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1_a = 0,
asc_1_b = 0,
asc_1_c = 0,
asc_1_d = 0,
asc_2_a = 0,
asc_2_b = 0,
asc_2_c = 0,
asc_2_d = 0,
beta_tt_a = 0,
beta_tt_b = 0,
beta_tt_c = 0,
beta_tt_d = 0,
beta_tc_a = 0,
beta_tc_b = 0,
beta_tc_c = 0,
beta_tc_d = 0,
# beta_hw_a =-0.0396,
# beta_hw_b =-0.0479,
# beta_ch_a =-0.7624,
# beta_ch_b =-2.1725,
delta_a_k1 = 0, # 0.0329,
delta_b_k1 = 0, # 0.0329,
delta_c_k1 = 0, # 0.0329,
delta_d_k1 =0, # 0.0329,
gamma_commute_a_k1 = 0,
gamma_commute_b_k1 = 0,
gamma_commute_c_k1 = 0,
gamma_commute_d_k1 = 0,
gamma_car_av_a_k1 = 0,
gamma_car_av_b_k1 = 0,
gamma_car_av_c_k1 = 0,
gamma_car_av_d_k1 = 0,
delta_a_k2 = 0,
delta_b_k2 = 0,
delta_c_k2 = 0,
delta_d_k2 = 0,
gamma_commute_a_k2 = 0,
gamma_commute_b_k2 = 0,
gamma_commute_c_k2 = 0,
gamma_commute_d_k2 = 0,
gamma_car_av_a_k2 = 0,
gamma_car_av_b_k2 = 0,
gamma_car_av_c_k2 = 0,
gamma_car_av_d_k2 = 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_2_a", "asc_2_b","delta_b","gamma_commute_b","gamma_car_av_b")

# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["asc_1" ]] = list( asc_1_a, asc_1_b,asc_1_c,asc_1_d)
lcpars[["asc_2" ]] = list( asc_2_a, asc_2_b,asc_2_c,asc_2_d)
lcpars[["beta_tt"]] = list(beta_tt_a, beta_tt_b,beta_tt_c,beta_tt_d)
lcpars[["beta_tc"]] = list(beta_tc_a, beta_tc_b,beta_tc_c,beta_tc_d)
# lcpars[["beta_hw"]] = list(beta_hw_a, beta_hw_b)
# lcpars[["beta_ch"]] = list(beta_ch_a, beta_ch_b)

V=list()
# V[["class_a"]] = delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability
# V[["class_b"]] = delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability

# V[["class_a"]] = (exp((delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)*Q[1,1])*exp((delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)*Q[2,1]))/((exp(delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)+1)*(exp(delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)+1))
# V[["class_b"]] = (exp((delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)*Q[1,2])*exp((delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)*Q[2,2]))/((exp(delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)+1)*(exp(delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)+1))
# V[["class_c"]] =(exp((delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)*Q[1,3])*exp((delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)*Q[2,3]))/((exp(delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)+1)*(exp(delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)+1))
# V[["class_d"]] =(exp((delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)*Q[1,4])*exp((delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)*Q[2,4]))/((exp(delta_a + gamma_commute_a*commute + gamma_car_av_a*car_availability)+1)*(exp(delta_b + gamma_commute_b*commute + gamma_car_av_b*car_availability)+1))

V[["class_a"]] = (exp((delta_a_k1 + gamma_commute_a_k1*commute + gamma_car_av_a_k1*car_availability)*Q[1,1])*exp((delta_a_k2 + gamma_commute_a_k2*commute + gamma_car_av_a_k2*car_availability)*Q[2,1]))/((exp(delta_a_k1 + gamma_commute_a_k1*commute + gamma_car_av_a_k1*car_availability)+1)*(exp(delta_a_k1 + gamma_commute_a_k1*commute + gamma_car_av_a_k1*car_availability)+1))
V[["class_b"]] = (exp((delta_b_k1 + gamma_commute_b_k1*commute + gamma_car_av_b_k1*car_availability)*Q[1,2])*exp((delta_b_k2 + gamma_commute_b_k2*commute + gamma_car_av_b_k2*car_availability)*Q[2,2]))/((exp(delta_b_k1 + gamma_commute_b_k1*commute + gamma_car_av_b_k1*car_availability)+1)*(exp(delta_b_k1 + gamma_commute_b_k1*commute + gamma_car_av_b_k1*car_availability)+1))
V[["class_c"]] =(exp((delta_c_k1 + gamma_commute_c_k1*commute + gamma_car_av_c_k1*car_availability)*Q[1,3])*exp((delta_c_k2 + gamma_commute_c_k2*commute + gamma_car_av_c_k2*car_availability)*Q[2,3]))/((exp(delta_c_k1 + gamma_commute_c_k1*commute + gamma_car_av_c_k1*car_availability)+1)*(exp(delta_c_k1 + gamma_commute_c_k1*commute + gamma_car_av_c_k1*car_availability)+1))
V[["class_d"]] =(exp((delta_d_k1+gamma_commute_d_k1*commute + gamma_car_av_d_k1*car_availability)*Q[1,4])*exp((delta_d_k2+gamma_commute_d_k2*commute + gamma_car_av_d_k2*car_availability)*Q[2,4]))/((exp(delta_d_k1+gamma_commute_d_k1*commute + gamma_car_av_d_k1*car_availability)+1)*(exp(delta_d_k2 + gamma_commute_d_k2*commute + gamma_car_av_d_k2*car_availability)+1))


### Settings for class allocation models
classAlloc_settings = list(
classes = c(class_a=1, class_b=2,class_c=3,class_d=4),
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()

### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2),
avail = list(alt1=1, alt2=1),
choiceVar = choice
)

### Loop over classes
for(s in 1:length(pi_values)){

### Compute class-specific utilities
V=list()
V[["alt1"]] = asc_1[[s]] + beta_tc[[s]]*tc1*Q[1,s] + beta_tt[[s]]*tt1*Q[2,s] #+ beta_hw[[s]]*hw1 + beta_ch[[s]]*ch1
V[["alt2"]] = asc_2[[s]] + beta_tc[[s]]*tc2*Q[1,s] + beta_tt[[s]]*tt2*Q[2,s] #+ beta_hw[[s]]*hw2 + beta_ch[[s]]*ch2

mnl_settings$utilities = V
mnl_settings$componentName = paste0("Class_",s)

### 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, 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)
}

# ################################################################# #
#### EM ESTIMATION FOR COVARIANCE MATRIX ####
# ################################################################# #

model=apollo_lcEM(apollo_beta, apollo_fixed, apollo_probabilities,
apollo_inputs, lcEM_settings = list(EMmaxIterations=100))

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model)

-------------------------------------------------------------------------------------------------------------------------------------------------
Model name : EM_LC_with_covariates_ANA
Model description : LC model with class allocation model on Swiss route choice data, EM algorithm, ANA
Model run at : 2024-04-05 19:00:00
Estimation method : EM algorithm (bfgs) -> Maximum likelihood (bfgs)
Model diagnosis : successful convergence
Optimisation diagnosis : Saddle point found
hessian properties : Some eigenvalues are positive and others negative
maximum eigenvalue : 0.778227
reciprocal of condition number : not calculated (Hessian is not negative definite)
Number of individuals : 388
Number of rows in database : 3492
Number of modelled outcomes : 3492

Number of cores used : 1
Model without mixing

LL(start) : -2420.47
LL (whole model) at equal shares, LL(0) : -2420.47
LL (whole model) at observed shares, LL(C) : -2420.39
LL(final, whole model) : -2295.6
Rho-squared vs equal shares : 0.0516
Adj.Rho-squared vs equal shares : 0.0351
Rho-squared vs observed shares : 0.0516
Adj.Rho-squared vs observed shares : 0.0392
AIC : 4671.19
BIC : 4917.52

LL(0,Class_1) : -2420.47
LL(final,Class_1) : -2428.88
LL(0,Class_2) : -2420.47
LL(final,Class_2) : -2574.11
LL(0,Class_3) : -2420.47
LL(final,Class_3) : -2621.55
LL(0,Class_4) : -2420.47
LL(final,Class_4) : -3487.33

Estimated parameters : 40
Time taken (hh:mm:ss) : 00:01:58.77
pre-estimation : 00:00:9.59
estimation : 00:00:58.6
post-estimation : 00:00:50.57
Iterations : 3 (EM) & 82 (bfgs)

Unconstrained optimisation.

Estimates:
Estimate s.e.
asc_1_a 0.06289 NA
asc_1_b -0.29744 NA
asc_1_c 0.03818 NA
asc_1_d 0.09051 NA
asc_2_a -0.06289 NA
asc_2_b 0.29744 NA
asc_2_c -0.03818 NA
asc_2_d -0.09050 NA
beta_tt_a 0.00000 NA
beta_tt_b 0.00000 NA
beta_tt_c -0.06053 NA
beta_tt_d -0.16335 NA
beta_tc_a 0.00000 NA
beta_tc_b -0.02635 NA
beta_tc_c 0.00000 NA
beta_tc_d -0.83472 NA
delta_a_k1 -8.77122 NA
delta_b_k1 0.00000 NA
delta_c_k1 -1.48338 NA
delta_d_k1 7.75759 NA
gamma_commute_a_k1 12.69659 NA
gamma_commute_b_k1 0.00000 NA
gamma_commute_c_k1 2.86563 NA
gamma_commute_d_k1 15.99644 NA
gamma_car_av_a_k1 24.42493 NA
gamma_car_av_b_k1 0.00000 NA
gamma_car_av_c_k1 -1.65684 NA
gamma_car_av_d_k1 -12.35994 NA
delta_a_k2 0.00000 NA
delta_b_k2 0.00000 NA
delta_c_k2 -0.71514 NA
delta_d_k2 7.75758 NA
gamma_commute_a_k2 0.00000 NA
gamma_commute_b_k2 0.00000 NA
gamma_commute_c_k2 2.25553 NA
gamma_commute_d_k2 15.99646 NA
gamma_car_av_a_k2 0.00000 NA
gamma_car_av_b_k2 0.00000 NA
gamma_car_av_c_k2 -0.20253 NA
gamma_car_av_d_k2 -12.35996 NA
t.rat.(0) Rob.s.e.
asc_1_a NA NA
asc_1_b NA NA
asc_1_c NA NA
asc_1_d NA NA
asc_2_a NA NA
asc_2_b NA NA
asc_2_c NA NA
asc_2_d NA NA
beta_tt_a NA NA
beta_tt_b NA NA
beta_tt_c NA NA
beta_tt_d NA NA
beta_tc_a NA NA
beta_tc_b NA NA
beta_tc_c NA NA
beta_tc_d NA NA
delta_a_k1 NA NA
delta_b_k1 NA NA
delta_c_k1 NA NA
delta_d_k1 NA NA
gamma_commute_a_k1 NA NA
gamma_commute_b_k1 NA NA
gamma_commute_c_k1 NA NA
gamma_commute_d_k1 NA NA
gamma_car_av_a_k1 NA NA
gamma_car_av_b_k1 NA NA
gamma_car_av_c_k1 NA NA
gamma_car_av_d_k1 NA NA
delta_a_k2 NA NA
delta_b_k2 NA NA
delta_c_k2 NA NA
delta_d_k2 NA NA
gamma_commute_a_k2 NA NA
gamma_commute_b_k2 NA NA
gamma_commute_c_k2 NA NA
gamma_commute_d_k2 NA NA
gamma_car_av_a_k2 NA NA
gamma_car_av_b_k2 NA NA
gamma_car_av_c_k2 NA NA
gamma_car_av_d_k2 NA NA
Rob.t.rat.(0)
asc_1_a NA
asc_1_b NA
asc_1_c NA
asc_1_d NA
asc_2_a NA
asc_2_b NA
asc_2_c NA
asc_2_d NA
beta_tt_a NA
beta_tt_b NA
beta_tt_c NA
beta_tt_d NA
beta_tc_a NA
beta_tc_b NA
beta_tc_c NA
beta_tc_d NA
delta_a_k1 NA
delta_b_k1 NA
delta_c_k1 NA
delta_d_k1 NA
gamma_commute_a_k1 NA
gamma_commute_b_k1 NA
gamma_commute_c_k1 NA
gamma_commute_d_k1 NA
gamma_car_av_a_k1 NA
gamma_car_av_b_k1 NA
gamma_car_av_c_k1 NA
gamma_car_av_d_k1 NA
delta_a_k2 NA
delta_b_k2 NA
delta_c_k2 NA
delta_d_k2 NA
gamma_commute_a_k2 NA
gamma_commute_b_k2 NA
gamma_commute_c_k2 NA
gamma_commute_d_k2 NA
gamma_car_av_a_k2 NA
gamma_car_av_b_k2 NA
gamma_car_av_c_k2 NA
gamma_car_av_d_k2 NA


Summary of class allocation for model
component :
Mean prob.
Class_1 0.2474
Class_2 0.1970
Class_3 0.2369
Class_4 0.3186

[stephanehess]

Hi

it's difficult to diagnose the issue here, but the large number of post-EM iterations is suspicious. Did you try estimation via apollo_estimate instead of EM first? That might better highlight any identification issues

Stephane