Error in the model command when using EM algorithms

[nimslkg]

Dear Prof. Hess & Palma,

I am trying to estimate a Mixed Multinomial Logit model in WTP space using EM algorithms because I want all the parameters to be random. The model has 11 parameters. 10 parameters are normally distributed and the cost attribute has a negative log-normal distribution. I used the estimates I obtained from the preference space MNL model as starting values for means.

Here is my code:

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)
library(dplyr)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  ="RPL_WTPspace_EM",
  modelDescr ="RPL model1_EM algorithm ",
  indivID    ="ID",
  mixing     = TRUE,
  nCores     = 10
)


# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database <- read_csv("Project data.csv") 


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

### Vector of parameters, including any that are kept fixed in estimation

apollo_beta <- c(
  asc_1                       = -0.4,
  asc_2                       = 0,
  mu_c                        = 0.01,
  mu_m2                       = 0.5, 
  mu_m3                       = 0.6,
  mu_t2                       = 0.2,
  mu_t3                       = 0.4,
  mu_s2                       = 0.2,
  mu_s3                       = 0.2,
  mu_s4                       = 0.1,
  mu_s5                       = -0.2,
  mu_s6                       = -0.5,
  mu_cost                     = -0.007,
  sigma_c                     = 0.1,
  sigma_m2                    = 0.3, 
  sigma_m3                    = 0.2,
  sigma_t2                    = 0.3,
  sigma_t3                    = 0.2,
  sigma_s2                    = 0.3,
  sigma_s3                    = 0.2,
  sigma_s4                    = 0.3,
  sigma_s5                    = 0.2,
  sigma_s6                    = 0.3,
  sigma_cost                  = 2.5
)


### 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")

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 1000,
  interNormDraws = c("xi_c", "xi_m2", "xi_m3", "xi_t2", "xi_t3", 
                     "xi_s2", "xi_s3", "xi_s4", "xi_s5", "xi_s6", "xi_cost")
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["b_c"]]   = mu_c + sigma_c * xi_c
  randcoeff[["b_m2"]]  = mu_m2 + sigma_m2 * xi_m2
  randcoeff[["b_m3"]]  = mu_m3 + sigma_m3 * xi_m3
  randcoeff[["b_t2"]]  = mu_t2 + sigma_t2 * xi_t2
  randcoeff[["b_t3"]]  = mu_t3 + sigma_t3 * xi_t3
  randcoeff[["b_s2"]]  = mu_s2 + sigma_s2 * xi_s2
  randcoeff[["b_s3"]]  = mu_s3 + sigma_s3 * xi_s3
  randcoeff[["b_s4"]]  = mu_s4 + sigma_s4 * xi_s4
  randcoeff[["b_s5"]]  = mu_s5 + sigma_s5 * xi_s5
  randcoeff[["b_s6"]]  = mu_s6 + sigma_s6 * xi_s6
  randcoeff[["b_cost"]]= -exp( mu_cost  + sigma_cost * xi_cost )
  
  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()
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  
  V[['alt1']] = asc_1 + b_cost * (b_c*a1_c + b_m2*a1_m2 + b_m3*a1_m3 + b_t2*a1_t2 + b_t3*a1_t3 
                                 + b_s2*a1_s2 + b_s3*a1_s3 + b_s4*a1_s4 + b_s5*a1_s5 + b_s6*a1_s6 + cost1)
  
  V[['alt2']] = asc_2 + b_cost * (b_c*a2_c + b_m2*a2_m2 + b_m3*a2_m3 + b_t2*a2_t2 + b_t3*a2_t3 
                                 + b_s2*a2_s2 + b_s3*a2_s3 + b_s4*a2_s4 + b_s5*a2_s5 + b_s6*a2_s6 + cost2)

  ### Define settings for MNL model component
  mnl_settings <- list(
    alternatives = c(alt1 = 1, alt2 = 2),
    avail        = list(alt1 = 1, alt2 = 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)
}

# ################################################################# #
#### EM ESTIMATION                                               ####
# ################################################################# #
mixEM_settings = list(transforms=list(function(x) (x), 
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) log(-x)))
model = apollo_mixEM(apollo_beta, apollo_fixed, apollo_probabilities,
                     apollo_inputs, mixEM_settings)

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

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

apollo_modelOutput(model)

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

apollo_saveOutput(model)

With the model command, I am getting an error like this:

Code: Select all

model = apollo_mixEM(apollo_beta, apollo_fixed, apollo_probabilities,
+                      apollo_inputs, mixEM_settings)
INFORMATION: The use of apollo_mixEM has a number of requirements. No checks are run for these, so the user needs
  to ensure these conditions are met by their model: 
1: This function is only suitable for single component models, i.e. no use of apollo_combineModels or
  manual multiplication of model components.
2: All parameters need to be random, and a full covariance matrix needs to be estimated/specified in
  apollo_randCoeff.
3: All random parameters need to be based on Normal distributions or transformations thereof.
4: With K random parameters, the order of the elements in 'apollo_beta' needs to be as follows: K
  means for the underlying Normals, followed by the elements of the lower triangle of the Cholesky
  matrix, by row.

Initialising EM algorithm
Validating inputs of likelihood function (apollo_probabilities)

Error in bgw::bgw_mle(calcR = f, betaStart = beta_var_val, calcJ = grad,  : 
  
User-provided value for bgw_settings[["maxFunctionEvals"]] is not valid.

As suggested in another post, I tried installing "bgw" package and then running the model. But, it did not work. Any help would be greatly appreciated.

Many thanks,
Manori

[stephanehess]

Hi

thanks for spotting this issue. This has been fixed in the development version (0.3.3) which you can download from http://apollochoicemodelling.com/code.html

Best wishes

Stephane

[nimslkg]

Dear Prof. Hess,

Thank you for your reply.

I downloaded the development version (0.3.3) and ran the apollo_mixEM example first. It ran without any errors or warnings. But when I ran my code, the model did converge; however, I am getting a warning as follows:

Code: Select all

WARNING: Singular Hessian, cannot calculate s.e. 
Could not write hessian to a file.
Warning messages:
1: In apollo_beta[1:4] <- mu :
  number of items to replace is not a multiple of replacement length
2: In apollo_beta[5:14] <- cholesky :
  number of items to replace is not a multiple of replacement length
3: In apollo_beta[1:4] <- mu :
  number of items to replace is not a multiple of replacement length
4: In apollo_beta[5:14] <- cholesky :
  number of items to replace is not a multiple of replacement length
  

Here is my code:

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

setwd("C:/Users/nimsl/Desktop")

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)
library(tidyverse)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  ="RPL_utilityspace_EM",
  modelDescr ="RPL model1_EM algorithm ",
  indivID    ="ID",
  nCores     = 10,
  outputDirectory = "output"
)


# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database <- read_csv("Project data.csv") 


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

### Vector of parameters, including any that are kept fixed in estimation

apollo_beta <- c(
  mu_c                        = 0.01,
  mu_m2                       = 0.5, 
  mu_m3                       = 0.6,
  mu_t2                       = 0.2,
  mu_t3                       = 0.4,
  mu_s2                       = 0.2,
  mu_s3                       = 0.2,
  mu_s4                       = 0.1,
  mu_s5                       = -0.2,
  mu_s6                       = -0.5,
  mu_logcost                  = -3,
  sigma_c                     = 0.1,
  sigma_c_m2                  = 0,
  sigma_m2                    = 0.3,
  sigma_c_m3                  = 0,
  sigma_m2_m3                 = 0,
  sigma_m3                    = 0.2,
  sigma_c_t2                  = 0,
  sigma_m2_t2                 = 0,
  sigma_m3_t2                 = 0,
  sigma_t2                    = 0.3,
  sigma_c_t3                  = 0,
  sigma_m2_t3                 = 0,
  sigma_m3_t3                 = 0,
  sigma_t2_t3                 = 0,
  sigma_t3                    = 0.2,
  sigma_c_s2                  = 0,
  sigma_m2_s2                 = 0,
  sigma_m3_s2                 = 0,
  sigma_t2_s2                 = 0,
  sigma_t3_s2                 = 0,
  sigma_s2                    = 0.3,
  sigma_c_s3                  = 0,
  sigma_m2_s3                 = 0,
  sigma_m3_s3                 = 0,
  sigma_t2_s3                 = 0,
  sigma_t3_s3                 = 0,
  sigma_s2_s3                 = 0,
  sigma_s3                    = 0.2,
  sigma_c_s4                  = 0,
  sigma_m2_s4                 = 0,
  sigma_m3_s4                 = 0,
  sigma_t2_s4                 = 0,
  sigma_t3_s4                 = 0,
  sigma_s2_s4                 = 0,
  sigma_s3_s4                 = 0,
  sigma_s4                    = 0.3,
  sigma_c_s5                  = 0,
  sigma_m2_s5                 = 0,
  sigma_m3_s5                 = 0,
  sigma_t2_s5                 = 0,
  sigma_t3_s5                 = 0,
  sigma_s2_s5                 = 0,
  sigma_s3_s5                 = 0,
  sigma_s4_s5                 = 0,
  sigma_s5                    = 0.2,
  sigma_c_s6                  = 0,
  sigma_m2_s6                 = 0,
  sigma_m3_s6                 = 0,
  sigma_t2_s6                 = 0,
  sigma_t3_s6                 = 0,
  sigma_s2_s6                 = 0,
  sigma_s3_s6                 = 0,
  sigma_s4_s6                 = 0,
  sigma_s5_s6                 = 0,
  sigma_s6                    = 0.3,
  sigma_c_logcost             = 0,
  sigma_m2_logcost            = 0,
  sigma_m3_logcost            = 0,
  sigma_t2_logcost            = 0,
  sigma_t3_logcost            = 0,
  sigma_s2_logcost            = 0,
  sigma_s3_logcost            = 0,
  sigma_s4_logcost            = 0,
  sigma_s5_logcost            = 0,
  sigma_s6_logcost            = 0.3,
  sigma_logcost               = 2.5
)


### 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()

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 1000,
  interUnifDraws = c(),
  interNormDraws = c("xi_c", "xi_m2", "xi_m3", "xi_t2", "xi_t3", 
                     "xi_s2", "xi_s3", "xi_s4", "xi_s5", "xi_s6", "xi_logcost")
)


### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["b_c"]]   = mu_c  + sigma_c    * xi_c
  randcoeff[["b_m2"]]  = mu_m2 + sigma_c_m2 * xi_c + sigma_m2    * xi_m2
  randcoeff[["b_m3"]]  = mu_m3 + sigma_c_m3 * xi_c + sigma_m2_m3 * xi_m2 + sigma_m3    * xi_m3
  randcoeff[["b_t2"]]  = mu_t2 + sigma_c_t2 * xi_c + sigma_m2_t2 * xi_m2 + sigma_m3_t2 * xi_m3 + sigma_t2    * xi_t2
  randcoeff[["b_t3"]]  = mu_t3 + sigma_c_t3 * xi_c + sigma_m2_t3 * xi_m2 + sigma_m3_t3 * xi_m3 + sigma_t2_t3 * xi_t2 + sigma_t3    * xi_t3
  randcoeff[["b_s2"]]  = mu_s2 + sigma_c_s2 * xi_c + sigma_m2_s2 * xi_m2 + sigma_m3_s2 * xi_m3 + sigma_t2_s2 * xi_t2 + sigma_t3_s2 * xi_t3 + sigma_s2    * xi_s2
  randcoeff[["b_s3"]]  = mu_s3 + sigma_c_s3 * xi_c + sigma_m2_s3 * xi_m2 + sigma_m3_s3 * xi_m3 + sigma_t2_s3 * xi_t2 + sigma_t3_s3 * xi_t3 + sigma_s2_s3 * xi_s2 + sigma_s3    * xi_s3
  randcoeff[["b_s4"]]  = mu_s4 + sigma_c_s4 * xi_c + sigma_m2_s4 * xi_m2 + sigma_m3_s4 * xi_m3 + sigma_t2_s4 * xi_t2 + sigma_t3_s4 * xi_t3 + sigma_s2_s4 * xi_s2 + sigma_s3_s4 * xi_s3 + sigma_s4    * xi_s4
  randcoeff[["b_s5"]]  = mu_s5 + sigma_c_s5 * xi_c + sigma_m2_s5 * xi_m2 + sigma_m3_s5 * xi_m3 + sigma_t2_s5 * xi_t2 + sigma_t3_s5 * xi_t3 + sigma_s2_s5 * xi_s2 + sigma_s3_s5 * xi_s3 + sigma_s4_s5 * xi_s4 + sigma_s5    * xi_s5
  randcoeff[["b_s6"]]  = mu_s6 + sigma_c_s6 * xi_c + sigma_m2_s6 * xi_m2 + sigma_m3_s6 * xi_m3 + sigma_t2_s6 * xi_t2 + sigma_t3_s6 * xi_t3 + sigma_s2_s6 * xi_s2 + sigma_s3_s6 * xi_s3 + sigma_s4_s6 * xi_s4 + sigma_s5_s6 * xi_s5 + sigma_s6 * xi_s6
  randcoeff[["b_cost"]]= -exp(mu_logcost  + sigma_c_logcost * xi_c + sigma_m2_logcost * xi_m2 + sigma_m3_logcost * xi_m3 + sigma_t2_logcost * xi_t2 + sigma_t3_logcost * xi_t3 + sigma_s2_logcost * xi_s2 + sigma_s3_logcost * xi_s3 + sigma_s4_logcost * xi_s4 + sigma_s5_logcost * xi_s5 + sigma_s6_logcost * xi_s6 + sigma_logcost * xi_logcost )
  
  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()
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  
  V[['alt1']] = b_c*a1_c + b_m2*a1_m2 + b_m3*a1_m3 + b_t2*a1_t2 + b_t3*a1_t3 
                                 + b_s2*a1_s2 + b_s3*a1_s3 + b_s4*a1_s4 + b_s5*a1_s5 + b_s6*a1_s6 + b_cost*cost1
  
  V[['alt2']] = b_c*a2_c + b_m2*a2_m2 + b_m3*a2_m3 + b_t2*a2_t2 + b_t3*a2_t3 
                                 + b_s2*a2_s2 + b_s3*a2_s3 + b_s4*a2_s4 + b_s5*a2_s5 + b_s6*a2_s6 + b_cost*cost2

  ### Define settings for MNL model component
  mnl_settings <- list(
    alternatives = c(alt1 = 1, alt2 = 2),
    avail        = list(alt1 = 1, alt2 = 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)
}


# ################################################################# #
#### EM ESTIMATION                                               ####
# ################################################################# #
mixEM_settings = list(transforms=list(function(x) (x), 
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) (x),
                                      function(x) log(-x)))
model = apollo_mixEM(apollo_beta, apollo_fixed, apollo_probabilities,
                     apollo_inputs, mixEM_settings)


# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

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

apollo_modelOutput(model)

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

apollo_saveOutput(model)

I would be grateful if you could help me identifying the issue in my code. Also, if cost is specified as negative lognormal, how do we obtain the mean and standard deviation of the original parameter. Should we use the delta method as follows:

Code: Select all

deltaMethod_settings <- list(operation="lognormal", parName1="mu_logcost", parName2="sigma_logcost")
apollo_deltaMethod(model, deltaMethod_settings)

Can we still use "lognormal" operation even if we specify the parameter as negative lognormal?

Many thanks,
Manori

[stephanehess]

Hi

this was another bug in apollo_mixEM which has now been fixed too.

and yes, that's correct re the delta method

Stephane