ANA_EM has no covariance matrix

[kiki]

Hi Stephane,

Thank you! The mathematical explanation about my desired model can be referred to "A discrete choice model with endogenous attribute attendance"
and "It’s not that I don’t care, I just don’t care very much: confounding between attribute non-attendance and taste heterogeneity".

As I understand, it propose an entry relating to attribute k in class s since MNL formulas for class membership becomes unpractical as the number of attributes grows. The class membership is constructed by a product of binary logit formulas defined by socio-demographics or other covariates labeled k. When k is attended, the k'th logit is exp(..)/(exp(...)+1), otherwise is 1/(exp(...)+1).

More details are elaborated in "It’s not that I don’t care, I just don’t care very much: confounding between attribute non-attendance and taste heterogeneity", in Methodology Section.

Best regards,
kiki

[stephanehess]

Hi

in the paper, the product of logit probabilities in equation 2 can be replaced in an equivalent way by a logit model where the utilities are sums of utilities from the individual attend/non-attend models.

The code below seems to work, you'd just need to expand it to add the continuous heterogeneity:

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 = "LC_with_covariates_ANA",
  modelDescr = "LC model with class allocation model on Swiss route choice data, ANA",
  indivID = "ID",
  nCores = 3,
  debug = TRUE
)

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

### Loading data from package
database = apollo_swissRouteChoiceData

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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_1           = 0,
                beta_attend_tt  = -0.1,
                beta_attend_tc  = -0.1,
                beta_attend_hw  = -0.1,
                beta_attend_ch  = -0.1,
                delta_attend_tt = 0,
                delta_attend_tc = 0,
                delta_attend_hw = 0,
                delta_attend_ch = 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()

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

apollo_lcPars=function(apollo_beta, apollo_inputs){
  
  lcpars = list()
  lcpars[["beta_tt"]] = list()
  for(s in 1:nrow(apollo_inputs$attend_matrix)) lcpars[["beta_tt"]][[s]] = apollo_inputs$attend_matrix[s,1] * beta_attend_tt
  lcpars[["beta_tc"]] = list()
  for(s in 1:nrow(apollo_inputs$attend_matrix)) lcpars[["beta_tc"]][[s]] = apollo_inputs$attend_matrix[s,2] * beta_attend_tc
  lcpars[["beta_hw"]] = list()
  for(s in 1:nrow(apollo_inputs$attend_matrix)) lcpars[["beta_hw"]][[s]] = apollo_inputs$attend_matrix[s,3] * beta_attend_hw
  lcpars[["beta_ch"]] = list()
  for(s in 1:nrow(apollo_inputs$attend_matrix)) lcpars[["beta_ch"]][[s]] = apollo_inputs$attend_matrix[s,4] * beta_attend_ch
  
  V=list()
  for(s in 1:nrow(apollo_inputs$attend_matrix)) V[[paste("Class_",s)]] = apollo_inputs$attend_matrix[s,1] * delta_attend_tt + 
    apollo_inputs$attend_matrix[s,2] * delta_attend_tc + 
    apollo_inputs$attend_matrix[s,3] * delta_attend_hw + 
    apollo_inputs$attend_matrix[s,4] * delta_attend_ch 
  
  ### Settings for class allocation models
  classAlloc_settings = list(
    utilities = V
  )
  
  lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
  
  return(lcpars)
}



# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #

# Create the ANA matrix
attribute_count = 4
attend_vector = c(0,1)
attend_matrix = expand.grid(replicate(attribute_count,attend_vector,simplify=FALSE))
apollo_inputs <- list(attend_matrix=attend_matrix)

apollo_inputs = apollo_validateInputs(recycle=TRUE)



# ################################################################# #
#### 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:16){

    ### Compute class-specific utilities
    V=list()
    V[["alt1"]] = asc_1 + beta_tt[[s]]*tt1 + beta_tc[[s]]*tc1 + beta_hw[[s]]*hw1 + beta_ch[[s]]*ch1
    V[["alt2"]] =         beta_tt[[s]]*tt2 + beta_tc[[s]]*tc2 + 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)
}

# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #

## Estimate model
model=apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities,
                      apollo_inputs)

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

apollo_modelOutput(model)

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

apollo_saveOutput(model)

[johnathan]

Hi

you need to put Q into apollo_inputs after the line

Code: Select all

apollo_inputs = apollo_validateInputs()

using

Code: Select all

apollo_inputs$Q = Q

and then everywhere else in apollo_lcPars and apollo_probabilities, you need to replace Q by apollo_inputs$Q

Stephane

Hi Stephane,

This recommendation addressed an issue I was encountering when trying to develop a script to estimate an LCL with a variable number of classes & pass the parameter n_classes to write variables dynamically for each LCL class. However, I now encounter an issue with lcpars during initialization, which needs to access apollo_inputs$n_classes, which is not defined until after apollo_validateInputs().

Alternatively, it was working when I assigned database$n_classes, but I don't think that makes sense to do for a parameter & gave the warning "Warning: numerical expression has 11010 elements: only the first usedWarning: numerical expression has 11010 elements: only the first used".

Overall, I am trying to set up code to be able to run over a range of LCL classes. How can I define n_classes in a way that it can be passed to apollo_lcPars & apollo_probabilities?

Thank you for your assistance & for developing Apollo.

You can find my code, output, & the error below:

# Initialize Apollo
apollo_initialise()

# Set up core control parameters
apollo_control <- list(
modelName = "LCL_CarChoice_Full",
modelDescr = "Latent Class Logit Model for Car Choice",
indivID = "ID",
panelData = TRUE,
nCores = 8,
mixing = FALSE, # No mixing in LCL
weights = 'Weight2018',
outputDirectory = "output/"
)

# Assign the dataframe to 'database'
database <- cbc_short_car_df

# Number of classes (to be changed dynamically)
n_classes <- 2

# Initialize the parameters list for each class
apollo_beta <- c()

# Loop to define parameters for each class
for (i in 1:n_classes) {
# Defining parameters for each class (b_price, b_acceleration, etc.)
apollo_beta <- c(apollo_beta,
setNames(i*-.1, paste0("b_price_", i)),
setNames(0, paste0("b_acceleration_", i)),
setNames(0, paste0("b_opcost_", i)),
setNames(0, paste0("b_bev_", i)),
setNames(0, paste0("b_bevRangeRel_", i)),
setNames(0, paste0("b_phev20_", i)),
setNames(0, paste0("b_phev40_", i)),
setNames(0, paste0("b_hev_", i)),
setNames(0, paste0("class_intercept_", i))
)
}

# Fixed parameters
apollo_fixed <- c("class_intercept_1")

apollo_lcPars <- function(apollo_beta, apollo_inputs){
lcpars = list()
# Loop to define the class-specific parameters for vehicle choice
print("Print value for 'apollo_inputs$n_classes':")
print(apollo_inputs$n_classes)
n <- apollo_inputs$n_classes
for (i in 1:n) {
lcpars[['b_price']][] <- get(paste0('b_price_', i))
lcpars[['b_acceleration']][] <- get(paste0('b_acceleration_', i))
lcpars[['b_opcost']][] <- get(paste0('b_opcost_', i))
lcpars[['b_bev']][] <- get(paste0('b_bev_', i))
lcpars[['b_bevRangeRel']][] <- get(paste0('b_bevRangeRel_', i))
lcpars[['b_phev20']][] <- get(paste0('b_phev20_', i))
lcpars[['b_phev40']][] <- get(paste0('b_phev40_', i))
lcpars[['b_hev']][] <- get(paste0('b_hev_', i))
# Parameters for the probability of class membership
lcpars[['class_intercept']][] <- get(paste0('class_intercept_', i))
}

# Define utilities for class allocation model (for class membership probability)
V = list()

# Class allocation
class_arg <- c()
for (j in 1:n) {
V[[paste0('class_', j)]] <- 0
V[[paste0('class_', j)]] <- V[[paste0('class_', j)]] + get(paste0('class_intercept_', j))
class_arg[paste0('class_', j)] <- j
}

# Settings for the class allocation model
classAlloc_settings = list(
classes = class_arg, # Dynamic class names
utilities = V
)

# Compute the class probabilities (pi_values)
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)

return(lcpars)
}

# Validate inputs
apollo_inputs <- apollo_validateInputs()
apollo_inputs$n_classes <- n_classes

# Define probabilities function for LCL with 2 classes
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()

# Call apollo_lcPars to get the class parameters and pi_values
lcpars = apollo_lcPars(apollo_beta, apollo_inputs)

# Loop over classes
print("Printing value of 'apollo_inputs$n_classes'")
print(apollo_inputs$n_classes)
for (s in 1:apollo_inputs$n_classes) {

# Define class-specific utilities for each alternative
V = list()

# For each alternative, define its utility for class `s`
V[["C1"]] =
get(paste0("b_price_", s)) * PriceC1 +
get(paste0("b_acceleration_", s)) * AccelerationC1 +
get(paste0("b_opcost_", s)) * OpCostC1 +
get(paste0("b_bev_", s)) * bevC1 +
get(paste0("b_bevRangeRel_", s)) * bevRangeRelC1 +
get(paste0("b_phev20_", s)) * phev20C1 +
get(paste0("b_phev40_", s)) * phev40C1 +
get(paste0("b_hev_", s)) * hevC1

V[["C2"]] =
get(paste0("b_price_", s)) * PriceC2 +
get(paste0("b_acceleration_", s)) * AccelerationC2 +
get(paste0("b_opcost_", s)) * OpCostC2 +
get(paste0("b_bev_", s)) * bevC2 +
get(paste0("b_bevRangeRel_", s)) * bevRangeRelC2 +
get(paste0("b_phev20_", s)) * phev20C2 +
get(paste0("b_phev40_", s)) * phev40C2 +
get(paste0("b_hev_", s)) * hevC2

V[["C3"]] =
get(paste0("b_price_", s)) * PriceC3 +
get(paste0("b_acceleration_", s)) * AccelerationC3 +
get(paste0("b_opcost_", s)) * OpCostC3 +
get(paste0("b_bev_", s)) * bevC3 +
get(paste0("b_bevRangeRel_", s)) * bevRangeRelC3 +
get(paste0("b_phev20_", s)) * phev20C3 +
get(paste0("b_phev40_", s)) * phev40C3 +
get(paste0("b_hev_", s)) * hevC3

# Define settings for MNL model component
mnl_settings = list(
alternatives = c(C1=1, C2=2, C3=3),
choiceVar = Choice, # Assuming `Choice` is the choice variable in your dataset
utilities = V
)

# 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 (combining all classes)
lc_settings = list(
inClassProb = P, # Probabilities from each class
classProb = pi_values # Class probabilities (e.g., `delta_1`, `delta_2`)
)

# Calculate the overall model probabilities using latent class model
P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)

# Weight individuals
P <- apollo_weighting(P, apollo_inputs, functionality)

### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}

# Estimate the model
# apollo_beta <- apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
lcl_model <- apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)

# Print results
apollo_modelOutput(lcl_model)

Apollo ignition sequence completed
All checks on apollo_control completed.
All checks on database completed.
[1] "Print value for 'apollo_inputs$n_classes':"
NULL
Error in 1:n : argument of length 0

[dpalma]

Hi,

I modified your code by adding line 48 (inside apollo_lcPars). It is kind of a nasty trick, but it should allow you to run your model. I cannot test the code as I don't have the database, but let us know if it doesn't work.

Code: Select all

# Initialize Apollo
apollo_initialise()

# Set up core control parameters
apollo_control <- list(
  modelName = "LCL_CarChoice_Full",
  modelDescr = "Latent Class Logit Model for Car Choice",
  indivID = "ID",
  panelData = TRUE,
  nCores = 8,
  mixing = FALSE, # No mixing in LCL
  weights = 'Weight2018',
  outputDirectory = "output/"
)

# Assign the dataframe to 'database'
database <- cbc_short_car_df

# Number of classes (to be changed dynamically)
n_classes <- 2

# Initialize the parameters list for each class
apollo_beta <- c()

# Loop to define parameters for each class
for (i in 1:n_classes) {
  # Defining parameters for each class (b_price, b_acceleration, etc.)
  apollo_beta <- c(apollo_beta,
                   setNames(i*-.1, paste0("b_price_", i)),
                   setNames(0, paste0("b_acceleration_", i)),
                   setNames(0, paste0("b_opcost_", i)),
                   setNames(0, paste0("b_bev_", i)),
                   setNames(0, paste0("b_bevRangeRel_", i)),
                   setNames(0, paste0("b_phev20_", i)),
                   setNames(0, paste0("b_phev40_", i)),
                   setNames(0, paste0("b_hev_", i)),
                   setNames(0, paste0("class_intercept_", i))
  )
}

# Fixed parameters
apollo_fixed <- c("class_intercept_1")

apollo_lcPars <- function(apollo_beta, apollo_inputs){
  lcpars = list()
  # Loop to define the class-specific parameters for vehicle choice
  print("Print value for 'apollo_inputs$n_classes':")
  if(is.null(apollo_inputs$n_classes)) apollo_inputs$n_classes <- get("n_classes", envir=.GlobalEnv)
  print(apollo_inputs$n_classes)
  n <- apollo_inputs$n_classes
  for (i in 1:n) {
    lcpars[['b_price']][] <- get(paste0('b_price_', i))
    lcpars[['b_acceleration']][] <- get(paste0('b_acceleration_', i))
    lcpars[['b_opcost']][] <- get(paste0('b_opcost_', i))
    lcpars[['b_bev']][] <- get(paste0('b_bev_', i))
    lcpars[['b_bevRangeRel']][] <- get(paste0('b_bevRangeRel_', i))
    lcpars[['b_phev20']][] <- get(paste0('b_phev20_', i))
    lcpars[['b_phev40']][] <- get(paste0('b_phev40_', i))
    lcpars[['b_hev']][] <- get(paste0('b_hev_', i))
    # Parameters for the probability of class membership
    lcpars[['class_intercept']][] <- get(paste0('class_intercept_', i))
  }
  
  # Define utilities for class allocation model (for class membership probability)
  V = list()
  
  # Class allocation
  class_arg <- c()
  for (j in 1:n) {
    V[[paste0('class_', j)]] <- 0
    V[[paste0('class_', j)]] <- V[[paste0('class_', j)]] + get(paste0('class_intercept_', j))
    class_arg[paste0('class_', j)] <- j
  }
  
  # Settings for the class allocation model
  classAlloc_settings = list(
    classes = class_arg, # Dynamic class names
    utilities = V
  )
  
  # Compute the class probabilities (pi_values)
  lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
  
  return(lcpars)
}

# Validate inputs
apollo_inputs <- apollo_validateInputs()
apollo_inputs$n_classes <- n_classes

# Define probabilities function for LCL with 2 classes
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()
  
  # Call apollo_lcPars to get the class parameters and pi_values
  lcpars = apollo_lcPars(apollo_beta, apollo_inputs)
  
  # Loop over classes
  print("Printing value of 'apollo_inputs$n_classes'")
  print(apollo_inputs$n_classes)
  for (s in 1:apollo_inputs$n_classes) {
    
    # Define class-specific utilities for each alternative
    V = list()
    
    # For each alternative, define its utility for class `s`
    V[["C1"]] =
      get(paste0("b_price_", s)) * PriceC1 +
      get(paste0("b_acceleration_", s)) * AccelerationC1 +
      get(paste0("b_opcost_", s)) * OpCostC1 +
      get(paste0("b_bev_", s)) * bevC1 +
      get(paste0("b_bevRangeRel_", s)) * bevRangeRelC1 +
      get(paste0("b_phev20_", s)) * phev20C1 +
      get(paste0("b_phev40_", s)) * phev40C1 +
      get(paste0("b_hev_", s)) * hevC1
    
    V[["C2"]] =
      get(paste0("b_price_", s)) * PriceC2 +
      get(paste0("b_acceleration_", s)) * AccelerationC2 +
      get(paste0("b_opcost_", s)) * OpCostC2 +
      get(paste0("b_bev_", s)) * bevC2 +
      get(paste0("b_bevRangeRel_", s)) * bevRangeRelC2 +
      get(paste0("b_phev20_", s)) * phev20C2 +
      get(paste0("b_phev40_", s)) * phev40C2 +
      get(paste0("b_hev_", s)) * hevC2
    
    V[["C3"]] =
      get(paste0("b_price_", s)) * PriceC3 +
      get(paste0("b_acceleration_", s)) * AccelerationC3 +
      get(paste0("b_opcost_", s)) * OpCostC3 +
      get(paste0("b_bev_", s)) * bevC3 +
      get(paste0("b_bevRangeRel_", s)) * bevRangeRelC3 +
      get(paste0("b_phev20_", s)) * phev20C3 +
      get(paste0("b_phev40_", s)) * phev40C3 +
      get(paste0("b_hev_", s)) * hevC3
    
    # Define settings for MNL model component
    mnl_settings = list(
      alternatives = c(C1=1, C2=2, C3=3),
      choiceVar = Choice, # Assuming `Choice` is the choice variable in your dataset
      utilities = V
    )
    
    # 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 (combining all classes)
  lc_settings = list(
    inClassProb = P, # Probabilities from each class
    classProb = pi_values # Class probabilities (e.g., `delta_1`, `delta_2`)
  )
  
  # Calculate the overall model probabilities using latent class model
  P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
  
  # Weight individuals
  P <- apollo_weighting(P, apollo_inputs, functionality)
  
  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}

# Estimate the model
# apollo_beta <- apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs)
lcl_model <- apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)

# Print results
apollo_modelOutput(lcl_model)

Best wishes,
David

[dpalma]

Hi,

I had a quick look at the code, and I believe you should not comment out the line including n_classes inside apollo_inputs after calling apollo_validateInputs. Please modify the code as shown as follows.

Code: Select all

# Validate inputs
apollo_inputs <- apollo_validateInputs()
apollo_inputs$n_classes <- n_classes

If that doesn't work, and you can share the data, please send it to my email D.Palma@leeds.ac.uk and I will have a closer look at the issue. In the longer term we will look into how to make things like this easier in Apollo.

Best wishes,
David

[johnathan]

Hi David,

Thank you for your response. I commented that line out during troubleshooting, but uncommenting it does not change the error.

I plan to reach out via email after organizing my work. Thanks again for your help.