Multinomial logit estimation gets stuck
Posted: 17 Sep 2021, 23:45
I am estimating a multinomial logit model. I use the simplest specification similar to Example 3, just to see if I specify everything correctly. My data has many individuals choosing many locations, and so I define the utilities and the availability objects (V and A) using a lot of get() statements, because there are hundreds of choices. Also, my data is "sparse", meaning that there are some alternatives that are never chosen.
The estimation takes about 30 minutes on a High Performance Cluster. The estimation converges without raising any problems After that, it is stuck at the following step (the last line):
After 15 hours of waiting, I killed the process because that did not seem reasonable (given that estimation itself took only 30 minutes).
My questions are:
- what does the LL(c) calculation do, exactly? is there a way to disable it in the seetings
- do you have a suggestion for why it might take up so much time/cause issues?
My code is below
### Load Apollo library
library(apollo)
library(data.table)
### Initialise code
apollo_initialise()
parallel::detectCores() #check how many cores the system has
### Set core controls
apollo_control = list(
modelName = "Logit_sch_location",
modelDescr = "Plain logit model",
indivID = "ind",
nCores = 12, #set the # of cores, should be number of cpus 1
outputDirectory = paste0(output_dir, "res_apollo/") #set the output directory
)
# # ################################################################# #
# #### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# # ################################################################# #
database = fread(paste0(data_dir, "mydata.csv"))
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(bprice = 1,
# family chars
bage = 1,
bagesq = -1,
bsex = -1,
bl1prealproppc = 1,
bschool = 1,
# location choice
bldistance = -1,
bcorn = -1,
bwheat = -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()
# ################################################################# #
#### 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))
# lists of alternatives (created directly from database)
list_alt <- gsub( "PRICE_", "", unique(grep("^PRICE_[^_]+$", colnames(apollo_inputs$database), perl = F, value=T)))
alternatives_set <- list_alt
names(alternatives_set) <- as.character(list_alt)
list_alt50 <- grep("^50.+$", alternatives_set, perl = F, value=T)
list_alt60 <- grep("^60.+$", alternatives_set, perl = F, value=T)
list_alt70 <- grep("^70.+$", alternatives_set, perl = F, value=T)
list_loc <- unique(substr(list_alt, 1, nchar(list_alt)-1)) # each alternative exists twice: school/non-sch, hence unique()
list_loc50 <- unique(substr(list_alt50, 1, nchar(list_alt)-1))
list_loc60 <- unique(substr(list_alt60, 1, nchar(list_alt)-1))
list_loc70 <- unique(substr(list_alt70, 1, nchar(list_alt)-1))
bpl_dummies <- grep("dBPL_MOM", names(apollo_inputs$database), perl = F, fixed = T, value=T)
betas_bpl_dummies <- gsub("dBPL_MOM", "bbpl", bpl_dummies)
expression_dummies <- paste0(paste(bpl_dummies, betas_bpl_dummies, sep= "*"), collapse = " + ")
### Create list of probabilities P
P = list() #probabilities of choice, will be filled automatically
V = list() #utility (except for the type-1 component), need to populate
A = list() # adjust availability based on year
### Create alternative specific constants and coefficients using interactions with socio-demographics
#inspired by this thread on Apollo forum: http://www.apollochoicemodelling.com/fo ... ?f=11&t=27
for (ylocsch_iter in list_alt){
sch_iter <- as.integer(substr(ylocsch_iter, nchar(ylocsch_iter), nchar(ylocsch_iter))) #last digit
loc_iter <- substr(ylocsch_iter, 1, nchar(ylocsch_iter)-1) #everything but the last digit
A[[ylocsch_iter]] = get( paste0("av_", ylocsch_iter) )
if (sch_iter == 0) {
V[[ylocsch_iter ]] = bldistance*get(paste0("LOG_DISTANCE_", ylocsch_iter)) + bcorn*get(paste0("distCorntoPot_", ylocsch_iter)) +
bwheat*get(paste0("distWheatoPot_", ylocsch_iter))
}
else if (sch_iter == 1) {
# no bsch because the origin FE saturate the model
V[[ylocsch_iter]] = bprice*get(paste0("PRICE_", ylocsch_iter)) +
bschool +
bage*AGE + bagesq*AGESQ + bsex*SEX +
bl1prealproppc*l1prealprop_pc +
bldistance*get(paste0("LOG_DISTANCE_", ylocsch_iter)) + bcorn*get(paste0("distCorntoPot_", ylocsch_iter)) +
bwheat*get(paste0("distWheatoPot_", ylocsch_iter))
}
}
### Define settings for MNL model component
mnl_settings = list(
alternatives = alternatives_set,
avail = A,
choiceVar = ylocsch_choice,
V = V
)
### Compute probabilities using NL model
P[["model"]] = apollo_mnl(mnl_settings, functionality)
### Take product across observations 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)
}
## In case you need to check the utility values, for some options, use this bit
# cols_50s009s1 <- grep("50s009s1", names(database), perl = F, fixed = T, value=T)
# cols_small <- c("ind", "BPL_MOM", "AGE", "l1prealprop_pc", "SEX", cols_50s009s1)
# cens_small <- database[,..cols_small]
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs,
estimate_settings = list(maxIterations=500))
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
The estimation takes about 30 minutes on a High Performance Cluster. The estimation converges without raising any problems After that, it is stuck at the following step (the last line):
Code: Select all
Computing covariance matrix using numerical methods (numDeriv).
0%....25%....50%....75%....100%
Negative definite Hessian with maximum eigenvalue: -27.102437
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
My questions are:
- what does the LL(c) calculation do, exactly? is there a way to disable it in the seetings
- do you have a suggestion for why it might take up so much time/cause issues?
My code is below
### Load Apollo library
library(apollo)
library(data.table)
### Initialise code
apollo_initialise()
parallel::detectCores() #check how many cores the system has
### Set core controls
apollo_control = list(
modelName = "Logit_sch_location",
modelDescr = "Plain logit model",
indivID = "ind",
nCores = 12, #set the # of cores, should be number of cpus 1
outputDirectory = paste0(output_dir, "res_apollo/") #set the output directory
)
# # ################################################################# #
# #### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# # ################################################################# #
database = fread(paste0(data_dir, "mydata.csv"))
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(bprice = 1,
# family chars
bage = 1,
bagesq = -1,
bsex = -1,
bl1prealproppc = 1,
bschool = 1,
# location choice
bldistance = -1,
bcorn = -1,
bwheat = -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()
# ################################################################# #
#### 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))
# lists of alternatives (created directly from database)
list_alt <- gsub( "PRICE_", "", unique(grep("^PRICE_[^_]+$", colnames(apollo_inputs$database), perl = F, value=T)))
alternatives_set <- list_alt
names(alternatives_set) <- as.character(list_alt)
list_alt50 <- grep("^50.+$", alternatives_set, perl = F, value=T)
list_alt60 <- grep("^60.+$", alternatives_set, perl = F, value=T)
list_alt70 <- grep("^70.+$", alternatives_set, perl = F, value=T)
list_loc <- unique(substr(list_alt, 1, nchar(list_alt)-1)) # each alternative exists twice: school/non-sch, hence unique()
list_loc50 <- unique(substr(list_alt50, 1, nchar(list_alt)-1))
list_loc60 <- unique(substr(list_alt60, 1, nchar(list_alt)-1))
list_loc70 <- unique(substr(list_alt70, 1, nchar(list_alt)-1))
bpl_dummies <- grep("dBPL_MOM", names(apollo_inputs$database), perl = F, fixed = T, value=T)
betas_bpl_dummies <- gsub("dBPL_MOM", "bbpl", bpl_dummies)
expression_dummies <- paste0(paste(bpl_dummies, betas_bpl_dummies, sep= "*"), collapse = " + ")
### Create list of probabilities P
P = list() #probabilities of choice, will be filled automatically
V = list() #utility (except for the type-1 component), need to populate
A = list() # adjust availability based on year
### Create alternative specific constants and coefficients using interactions with socio-demographics
#inspired by this thread on Apollo forum: http://www.apollochoicemodelling.com/fo ... ?f=11&t=27
for (ylocsch_iter in list_alt){
sch_iter <- as.integer(substr(ylocsch_iter, nchar(ylocsch_iter), nchar(ylocsch_iter))) #last digit
loc_iter <- substr(ylocsch_iter, 1, nchar(ylocsch_iter)-1) #everything but the last digit
A[[ylocsch_iter]] = get( paste0("av_", ylocsch_iter) )
if (sch_iter == 0) {
V[[ylocsch_iter ]] = bldistance*get(paste0("LOG_DISTANCE_", ylocsch_iter)) + bcorn*get(paste0("distCorntoPot_", ylocsch_iter)) +
bwheat*get(paste0("distWheatoPot_", ylocsch_iter))
}
else if (sch_iter == 1) {
# no bsch because the origin FE saturate the model
V[[ylocsch_iter]] = bprice*get(paste0("PRICE_", ylocsch_iter)) +
bschool +
bage*AGE + bagesq*AGESQ + bsex*SEX +
bl1prealproppc*l1prealprop_pc +
bldistance*get(paste0("LOG_DISTANCE_", ylocsch_iter)) + bcorn*get(paste0("distCorntoPot_", ylocsch_iter)) +
bwheat*get(paste0("distWheatoPot_", ylocsch_iter))
}
}
### Define settings for MNL model component
mnl_settings = list(
alternatives = alternatives_set,
avail = A,
choiceVar = ylocsch_choice,
V = V
)
### Compute probabilities using NL model
P[["model"]] = apollo_mnl(mnl_settings, functionality)
### Take product across observations 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)
}
## In case you need to check the utility values, for some options, use this bit
# cols_50s009s1 <- grep("50s009s1", names(database), perl = F, fixed = T, value=T)
# cols_small <- c("ind", "BPL_MOM", "AGE", "l1prealprop_pc", "SEX", cols_50s009s1)
# cens_small <- database[,..cols_small]
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs,
estimate_settings = list(maxIterations=500))
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)