ApolloChoiceModelling forum • Latent Class Joint Model

Page 1 of 2

Latent Class Joint Model

Posted: 26 Feb 2021, 13:02

by cybey

Hello everyone,

I have a new data set with dual response, i.e. the data set contains forced choices on the one hand, but also free choices where the alternative selected first is compared to a none option. Now I have learned that the scale parameter may differ for these two data sets. Therefore, I estimated a MNL and MIXL model with a parameter "mu_CBC" and "mu_DR", respectively, to capture these differences. In fact, in these joint models mu_CBC is significantly greater than one, meaning that scale differs.

Now I wanted to estimate a Latent Class model (using HB or classical estimation). The problem I encountered is that unlike the models covered in the manual, the submodels can no longer simply be written into a list (namely "P"), but I have submodels in the classes (or do I actually have classes in the submodels? I hope my assumption is correct - otherwise please correct me).
I solved the problem by creating a second list “P.componentName” in the while loop, that is while iterating over the classes. This list in then used for apollo_combineModels().

My apollo_probabilities() function for HB looks like this:

Code: Select all

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()
  
  ### Loop over classes
  s=1
  while(s<=2){
    
    ### Create list of probabilities P.componentName
    P.componentName <- list()
    
    ### Compute class-specific utilities
    V = list()
    
    ### CBC (Forced Choices)
    V[['alt1']] = b_asc_1 +
      b_Quota_Count_15[[s]] * Quota_Count2.1 + b_Quota_Count_25[[s]] * Quota_Count3.1 + b_Quota_Count_35[[s]] * Quota_Count4.1 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.1 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.1 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    V[['alt2']] = b_asc_2 +
      b_Quota_Count_15[[s]] * Quota_Count2.2 + b_Quota_Count_25[[s]] * Quota_Count3.2 + b_Quota_Count_35[[s]] * Quota_Count4.2 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.2 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.2 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.2 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
      b_Quota_Compensation[[s]] * Quota_Compensation.2
    
    V[['alt3']] = b_asc_3 +
      b_Quota_Count_15[[s]] * Quota_Count2.3 + b_Quota_Count_25[[s]] * Quota_Count3.3 + b_Quota_Count_35[[s]] * Quota_Count4.3 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.3 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.3 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.3 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
      b_Quota_Compensation[[s]] * Quota_Compensation.3
    
    curr.componentName = paste("CBC_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt1=1, alt2=2, alt3=3),
      avail         = list(alt1=1, alt2=1, alt3=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", mu_CBC),
      rows          = (CBC_type=="R"),
      componentName = curr.componentName
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[curr.componentName]] = apollo_mnl(mnl_settings, functionality)
    
    # ### Take product across observation for same individual
    # P.componentName[[curr.componentName]] = apollo_panelProd(P[[s]], apollo_inputs ,functionality)
    # 
    # ### Average across inter-individual draws within classes
    # P.componentName[[curr.componentName]] = apollo_avgInterDraws(P[[s]], apollo_inputs, functionality)
    
    
    ### DR (Free Choices)
    V[['alt0']] = 0
    
    V[['alt1']] = b_DR_asc_1[[s]] +
      b_Quota_Count_15[[s]] * Quota_Count2.1 + b_Quota_Count_25[[s]] * Quota_Count3.1 + b_Quota_Count_35[[s]] * Quota_Count4.1 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.1 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.1 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    curr.componentName = paste("CBC_DR_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt0=0, alt1=1),
      avail         = list(alt0=1, alt1=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", mu_DR),
      rows          = (CBC_type=="RD"),
      componentName = curr.componentName
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[curr.componentName]] = apollo_mnl(mnl_settings, functionality)
    
    # ### Take product across observation for same individual
    # P.componentName[[curr.componentName]] = apollo_panelProd(P[[s]], apollo_inputs ,functionality)
    # 
    # ### Average across inter-individual draws within classes
    # P.componentName[[curr.componentName]] = apollo_avgInterDraws(P[[s]], apollo_inputs, functionality)
    
    ### Combined model
    P[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality)
    
    s=s+1
  }
  
  ### 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)
}

For classical model estimation using MSL, it looks like this:

Code: Select all

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()
  
  ### Loop over classes
  s=1
  while(s<=2){
    
    ### Create list of probabilities P.componentName
    P.componentName <- list()
    
    ### Compute class-specific utilities
    V = list()
    
    ### CBC (Forced Choices)
    V[['alt1']] = b_asc_1 +
      b_Quota_Count_15[[s]] * Quota_Count2.1 + b_Quota_Count_25[[s]] * Quota_Count3.1 + b_Quota_Count_35[[s]] * Quota_Count4.1 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.1 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.1 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    V[['alt2']] = b_asc_2 +
      b_Quota_Count_15[[s]] * Quota_Count2.2 + b_Quota_Count_25[[s]] * Quota_Count3.2 + b_Quota_Count_35[[s]] * Quota_Count4.2 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.2 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.2 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.2 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
      b_Quota_Compensation[[s]] * Quota_Compensation.2
    
    V[['alt3']] = b_asc_3 +
      b_Quota_Count_15[[s]] * Quota_Count2.3 + b_Quota_Count_25[[s]] * Quota_Count3.3 + b_Quota_Count_35[[s]] * Quota_Count4.3 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.3 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.3 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.3 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
      b_Quota_Compensation[[s]] * Quota_Compensation.3
    
    curr.componentName = paste("CBC_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt1=1, alt2=2, alt3=3),
      avail         = list(alt1=1, alt2=1, alt3=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", mu_CBC),
      rows          = (CBC_type=="R"),
      componentName = curr.componentName
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[curr.componentName]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[curr.componentName]] = apollo_panelProd(P.componentName[[curr.componentName]], apollo_inputs ,functionality)
    
    ### Average across inter-individual draws within classes
    P.componentName[[curr.componentName]] = apollo_avgInterDraws(P.componentName[[curr.componentName]], apollo_inputs, functionality)
    
    
    ### DR (Free Choices)
    V[['alt0']] = 0
    
    V[['alt1']] = b_DR_asc_1[[s]] +
      b_Quota_Count_15[[s]] * Quota_Count2.1 + b_Quota_Count_25[[s]] * Quota_Count3.1 + b_Quota_Count_35[[s]] * Quota_Count4.1 +
      b_Quota_Max_Duration_1h[[s]] * Quota_Max_Duration2.1 + b_Quota_Max_Duration_2h[[s]] * Quota_Max_Duration3.1 + b_Quota_Max_Duration_3h[[s]] * Quota_Max_Duration4.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    curr.componentName = paste("CBC_DR_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt0=0, alt1=1),
      avail         = list(alt0=1, alt1=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", mu_DR),
      rows          = (CBC_type=="RD"),
      componentName = curr.componentName
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[curr.componentName]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[curr.componentName]] = apollo_panelProd(P.componentName[[curr.componentName]], apollo_inputs ,functionality)
    
    ### Average across inter-individual draws within classes
    P.componentName[[curr.componentName]] = apollo_avgInterDraws(P.componentName[[curr.componentName]], apollo_inputs, functionality)
    
    ### Combined model
    P[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality)
    
    s=s+1
  }
  
  ### 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)
}

I just wanted to make sure I did not make a mistake in my model specification.
In particular, I wanted to ask if the multiplication by mu_CBC and mu_DR is still correct for taking into account differences in scale?

I am looking forward to your answers.

Best
Nico

Re: Latent Class Joint Model

Posted: 05 Mar 2021, 18:05

by stephanehess

Nico

yes, this should be fine. In newer versions of Apollo, you can also use an additional setting for apollo_combineModels that will not retain the subcomponents, which would avoid the need for separate lists, but what you are doing is fine too

Stephane

Re: Latent Class Joint Model

Posted: 30 Mar 2021, 14:54

by cybey

Hi Stephane,

thank you very much for your quick reply.
Unfortunately I have to consult you again, because I get an error message when estimating the model (classical estimation).

With Apollo v0.2.4:

Code: Select all

Pre-processing likelihood function...
Preparing workers for multithreading...
Error in checkForRemoteErrors(lapply(cl, recvResult)) : 
  27 nodes produced errors; first error: attempt to use zero-length variable name

With the "nightly build" v0.2.5:

Code: Select all

Computing covariance matrix using numerical methods (numDeriv).
 0%....25%....50%....75%....100%
Negative definite Hessian with maximum eigenvalue: -29.134597
Computing score matrix...

Summary of class allocation for LC model component :
         Mean prob.
Class_1      0.3463
Class_2      0.6537

Calculating LL(0)...
Error in log(x) : non-numeric argument to mathematical function

Estimation worked some time with v0.2.5. But I think after updating other R packages, I get the error message mentioned above.

Here is the whole code of the model:

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 = "flexQgrid_LC_2Classes",
  modelDescr = "Latent class model on flexQgrid panel data (classical estimation)",
  indivID ="sys_RespNum",  
  nCores = 27, # parallel::detectCores()
  HB = FALSE,
  mixing = FALSE,
  workInLogs = FALSE # TRUE for higher numeric stability at the expense of computational time; works with the sum of log-probabilities, rather than the product of probabilities
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #
input.file <- "./input/CBC_results_merged.wide.rds"
database <- readRDS(input.file)

### Dual response
database[which(database$CBC_type == "RD" & database$Choice <7), "Choice"] <- 0
database[which(database$CBC_type == "RD" & database$Choice >=7), "Choice"] <- 1

### Dummy for first row a repondent
database[, "sys_RespNum.first"] <- rep(0, times=nrow(database))
database[match( unique(database$sys_RespNum), database$sys_RespNum), "sys_RespNum.first"] <- 1

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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
  # ----------------------------------------------------------------- #
  #---- Generic parameters
  # ----------------------------------------------------------------- #
  ### Scaling parameters
  mu_CBC = 0,
  mu_DR = 0,
  
  ### Alternative specific constants (ASC)
  b_asc_1 = 0,
  b_asc_2 = 0,
  b_asc_3 = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Class Assignment (A)
  # ----------------------------------------------------------------- #
  delta_A = 0,
  
  ### Sociodemographics
  gamma_A_Gender = 0,
  gamma_A_Age = 0,
  gamma_A_Education = 0,
  gamma_A_IncomeClasses = 0,
  
  gamma_A_ExistingApps = 0,
  gamma_A_ApartmentOwner = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Class Assignment (B)
  # ----------------------------------------------------------------- #
  delta_B = 0.1,
  
  ### Sociodemographics
  gamma_B_Gender = 0,
  gamma_B_Age = 0,
  gamma_B_Education = 0,
  gamma_B_IncomeClasses = 0,
  
  gamma_B_ExistingApps = 0,
  gamma_B_ApartmentOwner = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Intra-Class-Parameters (A)
  # ----------------------------------------------------------------- #
  ### DR_asc
  A_b_DR_asc_0 = 0,
  A_b_DR_asc_1 = 0,
  
  
  ### Quota_Count
  A_b_Quota_Count = 0,
  
  
  ### Quota_Max_Duration
  A_b_Quota_Max_Duration = 0,
  
  
  ### Quota_Period_midmorning
  A_b_Quota_Period_midmorning = 0,
  
  ### Quota_Period_afternoon
  A_b_Quota_Period_afternoon = 0,
  
  ### Quota_Period_evening
  A_b_Quota_Period_evening = 0,
  
  ### Quota_Period_night
  A_b_Quota_Period_night = 0,
  
  
  ### Quota_Compensation
  A_b_Quota_Compensation = -0.1,
  
  
  # ----------------------------------------------------------------- #
  #---- Intra-Class-Parameters (B)
  # ----------------------------------------------------------------- #
  ### DR_asc
  B_b_DR_asc_0 = 0,
  B_b_DR_asc_1 = 0,
  
  
  ### Quota_Count
  B_b_Quota_Count = 0,
  
  
  ### Quota_Max_Duration
  B_b_Quota_Max_Duration = 0,
  
  
  ### Quota_Period_midmorning
  B_b_Quota_Period_midmorning = 0,
  
  ### Quota_Period_afternoon
  B_b_Quota_Period_afternoon = 0,
  
  ### Quota_Period_evening
  B_b_Quota_Period_evening = 0,
  
  ### Quota_Period_night
  B_b_Quota_Period_night = 0,
  
  
  ### Quota_Compensation
  B_b_Quota_Compensation = -0.1
)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_bwtp_fixed = c() if none
apollo_fixed = c("mu_DR",
                 "b_asc_3",
                 
                 "delta_A",
                 
                 "gamma_A_Gender",
                 "gamma_A_Age",
                 "gamma_A_Education",
                 "gamma_A_IncomeClasses",
                 
                 "gamma_A_ExistingApps",
                 "gamma_A_ApartmentOwner",
                 
                 "A_b_DR_asc_0",
                 "B_b_DR_asc_0"
                 )
  

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

apollo_lcPars = function(apollo_beta, apollo_inputs){
  lcpars = list()
  
  lcpars[["b_DR_asc_0"]] = list(A_b_DR_asc_0, B_b_DR_asc_0)
  lcpars[["b_DR_asc_1"]] = list(A_b_DR_asc_1, B_b_DR_asc_1)
  
  lcpars[["b_Quota_Count"]] = list(A_b_Quota_Count, B_b_Quota_Count)
  
  lcpars[["b_Quota_Max_Duration"]] = list(A_b_Quota_Max_Duration, B_b_Quota_Max_Duration)
  
  lcpars[["b_Quota_Period_midmorning"]] = list(A_b_Quota_Period_midmorning, B_b_Quota_Period_midmorning)
  lcpars[["b_Quota_Period_afternoon"]] = list(A_b_Quota_Period_afternoon, B_b_Quota_Period_afternoon)
  lcpars[["b_Quota_Period_evening"]] = list(A_b_Quota_Period_evening, B_b_Quota_Period_evening)
  lcpars[["b_Quota_Period_night"]] = list(A_b_Quota_Period_night, B_b_Quota_Period_night)
  
  lcpars[["b_Quota_Compensation"]] = list(A_b_Quota_Compensation, B_b_Quota_Compensation)
  
  V=list()
  V[["class_A"]] = delta_A +
    gamma_A_Gender * COV_Gender +
    gamma_A_Age * COV_Age_centered +
    gamma_A_Education * COV_Education_log_centered +
    gamma_A_IncomeClasses * COV_IncomeClasses_centered +
    
    gamma_A_ExistingApps * COV_ExistingApps +
    gamma_A_ApartmentOwner * COV_ApartmentOwner
  
  V[["class_B"]] = delta_B +
    gamma_B_Gender * COV_Gender +
    gamma_B_Age * COV_Age_centered +
    gamma_B_Education * COV_Education_log_centered +
    gamma_B_IncomeClasses * COV_IncomeClasses_centered +
    
    gamma_B_ExistingApps * COV_ExistingApps +
    gamma_B_ApartmentOwner * COV_ApartmentOwner
  
  mnl_settings = list(
    alternatives = c(class_A=1, class_B=2), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  
  lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
  lcpars[["pi_values"]] = apollo_firstRow(lcpars[["pi_values"]], apollo_inputs)
  
  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()
  
  ### Loop over classes
  for (s in 1:2) {
    
    ### Create list of probabilities P.componentName
    P.componentName <- list()
    
    ### Compute class-specific utilities
    
    ### CBC (Forced Choices)
    V = list()
    
    V[['alt1']] = b_asc_1 +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    V[['alt2']] = b_asc_2 +
      b_Quota_Count[[s]] * Quota_Count.2 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.2 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
      b_Quota_Compensation[[s]] * Quota_Compensation.2
    
    V[['alt3']] = b_asc_3 +
      b_Quota_Count[[s]] * Quota_Count.3 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.3 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
      b_Quota_Compensation[[s]] * Quota_Compensation.3
    
    componentName1 = paste("CBC_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt1=1, alt2=2, alt3=3),
      avail         = list(alt1=1, alt2=1, alt3=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", (1+mu_CBC)),
      rows          = (CBC_type=="R"),
      componentName = componentName1
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName1]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName1]] = apollo_panelProd(P.componentName[[componentName1]], apollo_inputs, functionality)
    
    
    ### CBC (Free Choices)
    V = list()
    
    V[['alt0']] = b_DR_asc_0[[s]]
    
    V[['alt1']] = b_DR_asc_1[[s]] +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    componentName2 = paste("CBC_DR_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt0=0, alt1=1),
      avail         = list(alt0=1, alt1=1),
      choiceVar     = Choice,
      V             = lapply(V, "*", (1+mu_DR)),
      rows          = (CBC_type=="RD"),
      componentName = componentName2
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName2]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName2]] = apollo_panelProd(P.componentName[[componentName2]], apollo_inputs, functionality)
    
    
    ### Combined model
    P[[s]] = apollo_combineModels(P.componentName, 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                                            ####
# ################################################################# #

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

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

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

apollo_modelOutput(model)

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

saveOutput_settings = list(
  printClassical = TRUE,
  printPVal = TRUE,
  printDiagnostics = TRUE,
  printCovar = TRUE,
  printCorr = TRUE,
  printOutliers = 100,
  printChange = TRUE,
  saveEst = TRUE,
  saveCov = TRUE,
  saveCorr = TRUE,
  saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)

If needed, I can also send you the data via e-mail.
Thanks in advance.

Nico

Re: Latent Class Joint Model

Posted: 08 Apr 2021, 11:04

by stephanehess

Nico

could you share your data and code with us offline, please. Then we can look into it

Stephane

Re: Latent Class Joint Model

Posted: 23 Apr 2021, 08:29

by cybey

Good morning Stephane,

I have sent you an email with the code + data. By the way, the latent class model works fine if I don't estimate a joint model, i.e. omit the dual response part. However, I am still a little lost on what is causing the problem with the joint model.

Nico

Re: Latent Class Joint Model

Posted: 18 May 2021, 17:09

by stephanehess

Hi Nico

the problem is caused by the fact that your classes inside the latent class model are themselves multi-component models, and Apollo then creates an error when it comes to the calculation of the log-likelihood at zero. Until we come up with a more automatic solution, there are two easy ways of avoiding the problem.

1. The first solution is that you use a different list from P for the within class model. This means three changes to your code:

Code: Select all

  ### Create list of probabilities P
  P_within <- list()
  P <- list()

instead of

Code: Select all

  ### Create list of probabilities P
  P <- list()

and

Code: Select all

P_within[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality)

instead of

Code: Select all

P[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality)

and

Code: Select all

lc_settings = list(inClassProb = P_within, classProb=pi_values)

instead of

Code: Select all

lc_settings = list(inClassProb = P, classProb=pi_values)

2. The second solution is easier to implement, and all you need to do is to replace

Code: Select all

P[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality)

by

Code: Select all

P[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality,asList=FALSE)

That second solution means that Apollo will only retain the combined model inside each class, rather than your two separate submodels within each class

Stephane

Re: Latent Class Joint Model

Posted: 19 May 2021, 17:06

by cybey

Hi Stephane,

thank you very much for your response. I have tried both ways and they work.

That I have not thought of it myself ...

Anyway: What I actually wanted to do is to link the class allocation and latent variables, with the latter being modelled as random coefficients. While I managed to implement a hybrid choice joint mixed logit model, I failed with the hybrid choice joint latent class model. Or let's say I mistrust my model specification.

My attempt looks as follows:

(1) First, I have defined the LV in apollo_randCoeff(). For example:

Code: Select all

randcoeff[["LV_NACount"]] = gamma_Age_NACount * COV_Age_centered +
    gamma_Education_NACount * COV_Education_log_centered +
    gamma_ExistingApps_NACount * COV_ExistingApps +
    sigma_NACount * eta_NACount

(2) Afterwards, I have defined the latent class components. For example:

Code: Select all

V=list()
  V[["class_A"]] = delta_A +
    gamma_A_NACount * LV_NACount

V[["class_B"]] = delta_B +
    gamma_B_NACount * LV_NACount

V[["class_C"]] = delta_C +
    gamma_C_NACount * LV_NACount

(3) Then it is time to define apollo_probabilities().

Measurement component:

Code: Select all

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

# ----------------------------------------------------------------- #
#---- Likelihood of indicators
# ----------------------------------------------------------------- #
### NACount
normalDensity_settings1 = list(outcomeNormal = NonAttendance_r1, 
                               xNormal       = zeta_NACount_r1_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r1, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r1_NACount")

normalDensity_settings2 = list(outcomeNormal = NonAttendance_r2, 
                               xNormal       = zeta_NACount_r2_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r2, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r2_NACount")

normalDensity_settings3 = list(outcomeNormal = NonAttendance_r3, 
                               xNormal       = zeta_NACount_r3_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r3, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r3_NACount")

normalDensity_settings4 = list(outcomeNormal = NonAttendance_r4, 
                               xNormal       = zeta_NACount_r4_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r4, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r4_NACount")

normalDensity_settings5 = list(outcomeNormal = NonAttendance_r5, 
                               xNormal       = zeta_NACountIrrelevance_r5_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r5, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r5_NACount")

normalDensity_settings6 = list(outcomeNormal = NonAttendance_r6, 
                               xNormal       = zeta_NACountIrrelevance_r6_NACount * LV_NACount, 
                               mu            = 0, 
                               sigma         = sigma_NACount_r6, 
                               rows=(sys_RespNum.first==1),
                               componentName = "NACount_r6_NACount")
                               
### NACount
  P[["indic_NACount_NACount_r1"]] = apollo_normalDensity(normalDensity_settings1, functionality)
  P[["indic_NACount_NACount_r2"]] = apollo_normalDensity(normalDensity_settings2, functionality)
  P[["indic_NACount_NACount_r3"]] = apollo_normalDensity(normalDensity_settings3, functionality)
  P[["indic_NACount_NACount_r4"]] = apollo_normalDensity(normalDensity_settings4, functionality)
  P[["indic_NACount_NACount_r5"]] = apollo_normalDensity(normalDensity_settings5, functionality)
  P[["indic_NACount_NACount_r6"]] = apollo_normalDensity(normalDensity_settings6, functionality)

### Combined model
P = apollo_combineModels(P, apollo_inputs, functionality)
  
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
  
### Rename model
names(P)[which(names(P)=="model")] <- "Measurement_model"

Latent Class component:

Code: Select all

# ----------------------------------------------------------------- #
#---- LC model
# ----------------------------------------------------------------- #
  ### Create list of probabilities P
  P_within <- list() # list for within class models
  
  ### Loop over classes
  for (s in 1:3) {
    
    ### Create list of probabilities P.componentName
    P.componentName <- list()
    
    ### Compute class-specific utilities
    
    ### CBC (Forced Choices)
    V = list()
    
    V[['alt1']] = b_asc_1 +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    V[['alt2']] = b_asc_2 +
      b_Quota_Count[[s]] * Quota_Count.2 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.2 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
      b_Quota_Compensation[[s]] * Quota_Compensation.2
    
    V[['alt3']] = b_asc_3 +
      b_Quota_Count[[s]] * Quota_Count.3 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.3 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
      b_Quota_Compensation[[s]] * Quota_Compensation.3
    
    componentName1 = paste("CBC_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt1=1, alt2=2, alt3=3),
      avail         = list(alt1=1, alt2=1, alt3=1),
      choiceVar     = Choice,
      V             = V,
      rows          = (CBC_type=="R"),
      componentName = componentName1
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName1]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName1]] = apollo_panelProd(P.componentName[[componentName1]], apollo_inputs, functionality)
    
    
    ### CBC (Free Choices)
    V = list()
    
    V[['alt0']] = b_DR_asc_0[[s]]
    
    V[['alt1']] = b_DR_asc_1[[s]] +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    componentName2 = paste("CBC_DR_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt0=0, alt1=1),
      avail         = list(alt0=1, alt1=1),
      choiceVar     = Choice,
      V             = V,
      rows          = (CBC_type=="RD"),
      componentName = componentName2
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName2]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName2]] = apollo_panelProd(P.componentName[[componentName2]], apollo_inputs, functionality)
    
    ### Combined model
    P_within[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality, asList=FALSE)
  }
  
  ### Compute latent class model probabilities
  lc_settings = list(inClassProb = P_within, classProb=pi_values)
  P[["LC_model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)

Combine the model components:

Code: Select all

# ----------------------------------------------------------------- #
#---- Combine model components
# ----------------------------------------------------------------- #
  ### Combined model
  P = apollo_combineModels(P, apollo_inputs, 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)

With this specification, Rstudio has started to estimate. However, as already mentioned, I suspect that something is wrong. I suspect that the error is in the order of the function calls.
I know that the estimation involves maximising the joint likelihood of the observed sequence of choices and the observed answers to the attitudinal questions. Yet, I do not know how the implement this in Apollo. I've tried several model specifications, looked in the manual, and at the Apollo examples.

I look forward to your tips!

Nico

Re: Latent Class Joint Model

Posted: 19 May 2021, 19:47

by stephanehess

Nico

so you want the class allocation model to be a function of the LV, correct?

can you copy your whole code in one block, please? Then I'll tell you how to adjust it.

Stephane

Re: Latent Class Joint Model

Posted: 10 Jun 2021, 07:47

by cybey

Hi Stephane,

sorry for my late response. I had not received an email that a reply to the forum entry had been created - or I deleted it unintentionally. Anyway: Yes, I would like the class allocation model to be a function of the LV.

But first a few words about the model: In addition to respondents' preferences (using the choice model), we would like to investigate why respondents ignored certain attributes or considered them less important in their choices (latent variable model). To do this, we queried three latent variables per attribute:
(1) The attribute is too complex to understand.
(2) It is not worth the effort (e.g., in terms of time) to think about the impact the attribute has on a respondent's daily life.
(3) The attribute is in fact irrelevant.

Factor analysis showed that the first two LVs cannot be separated, which is why we merged them. That leaves us with two latent variables:
(1) ComplexityCosts
(2) Irrelevance

We would like to specify the latent variable model with two stages, i.e. it should look like this:

Non-attendance is to be used for class allocation.

Below is the complete code:

Code: Select all

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

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

### Load Apollo library
library(apollo)
library(corrplot)
library(Hmisc)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName = "flexQgrid_ICLV_LC_3Classes_DR1_CutOff6_SO2000",
  modelDescr = "Latent class model on flexQgrid panel data (classical estimation)",
  indivID ="sys_RespNum",  
  nCores = 60, # parallel::detectCores()
  HB = FALSE,
  mixing = TRUE,
  workInLogs = FALSE # TRUE for higher numeric stability at the expense of computational time; works with the sum of log-probabilities, rather than the product of probabilities
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #
input.file <- "./input/CBC_results_merged.wide.rds"
database <- readRDS(input.file)

### Dual response
database[which(database$CBC_type == "RD" & database$Choice <6), "Choice"] <- 0
database[which(database$CBC_type == "RD" & database$Choice >=6), "Choice"] <- 1

### Dummy for first row of a respondent
database[, "sys_RespNum.first"] <- rep(0, times=nrow(database))
database[match( unique(database$sys_RespNum), database$sys_RespNum), "sys_RespNum.first"] <- 1

### Mean centering of indicators
NonAttendance_colIndices <- grep("NonAttendance", colnames(database))
database[,NonAttendance_colIndices] <- apply(database[,NonAttendance_colIndices], 2, scale, center=TRUE, scale=FALSE)

### Load priors from step 1
# priors.step1 = list()
# priors.step1[["FC"]] = readRDS(file = "./priors_step2/FC_step1.rds") # A vector of starting values for the fixed parameters
# priors.step1[["svN"]] = readRDS(file = "./priors_step2/svN_step1.rds") # A vector of starting values for the means of the underlying normals for the random parameters
# priors.step1[["pvMatrix"]] = readRDS(file = "./priors_step2/pvMatrix_step1.rds") # A custom prior variance-covariance matrix to be used in estimation

# ################################################################# #
#### PRE-PROCESSING DATA                                         ####
# ################################################################# #
covariates.colnames <- c(
  "COV_Gender",
  "COV_Age_centered",
  "COV_Education_log_centered",
  "COV_IncomeClasses_log_centered",
  
  "COV_HouseholdSize_centered",
  "COV_ApartmentOwner",
  "COV_ExistingApps",
  "COV_CurrentMix",
  "COV_EnergyDecisionMaker",
  "COV_CheapTalk"
)

covariates <- database[,c("sys_RespNum",covariates.colnames)]
covariates <- covariates[!duplicated(covariates$sys_RespNum),] # only keep one row per respondent
covariates <- subset(covariates, select=-c(sys_RespNum))

## Corrplot
covariates.correlation = cor(x = covariates, method="pearson")
corrplot(covariates.correlation, method = "number")

covariates.correlation.p = rcorr(as.matrix(covariates), type="pearson") # "pearson" or "spearman"

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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
  # ----------------------------------------------------------------- #
  #---- Generic parameters
  # ----------------------------------------------------------------- #
  ### Alternative specific constants (ASC)
  b_asc_1 = 0,
  b_asc_2 = 0,
  b_asc_3 = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Class Assignment (A)
  # ----------------------------------------------------------------- #
  delta_A = 0,
  
  ### Sociodemographics
  # gamma_A_Gender = 0,
  # gamma_A_Age = 0,
  # gamma_A_Education = 0,
  # gamma_A_IncomeClasses = 0,
  # 
  # gamma_A_ExistingApps = 0,
  # gamma_A_ApartmentOwner = 0,
  
  ### Influence of LV on class assignment
  gamma_A_NACount = 0,
  gamma_A_NADuration = 0,
  gamma_A_NAPeriod = 0,
  gamma_A_NACompensation = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Class Assignment (B)
  # ----------------------------------------------------------------- #
  delta_B = 0.1,
  
  ### Sociodemographics
  # gamma_B_Gender = 0,
  # gamma_B_Age = 0,
  # gamma_B_Education = 0,
  # gamma_B_IncomeClasses = 0,
  # 
  # gamma_B_ExistingApps = 0,
  # gamma_B_ApartmentOwner = 0,
  
  ### Influence of LV on class assignment
  gamma_B_NACount = 0,
  gamma_B_NADuration = 0,
  gamma_B_NAPeriod = 0,
  gamma_B_NACompensation = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Class Assignment (C)
  # ----------------------------------------------------------------- #
  delta_C = 0.1,
  
  ### Sociodemographics
  # gamma_C_Gender = 0,
  # gamma_C_Age = 0,
  # gamma_C_Education = 0,
  # gamma_C_IncomeClasses = 0,
  # 
  # gamma_C_ExistingApps = 0,
  # gamma_C_ApartmentOwner = 0,
  
  ### Influence of LV on class assignment
  gamma_C_NACount = 0,
  gamma_C_NADuration = 0,
  gamma_C_NAPeriod = 0,
  gamma_C_NACompensation = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Intra-Class-Parameters (A)
  # ----------------------------------------------------------------- #
  ### DR_asc
  A_b_DR_asc_0 = 0,
  A_b_DR_asc_1 = 0,
  
  
  ### Quota_Count
  A_b_Quota_Count = 0,
  
  
  ### Quota_Max_Duration
  A_b_Quota_Max_Duration = 0,
  
  
  ### Quota_Period_midmorning
  A_b_Quota_Period_midmorning = 0,
  
  ### Quota_Period_afternoon
  A_b_Quota_Period_afternoon = 0,
  
  ### Quota_Period_evening
  A_b_Quota_Period_evening = 0,
  
  ### Quota_Period_night
  A_b_Quota_Period_night = 0,
  
  
  ### Quota_Compensation
  A_b_Quota_Compensation = -0.1,
  
  
  # ----------------------------------------------------------------- #
  #---- Intra-Class-Parameters (B)
  # ----------------------------------------------------------------- #
  ### DR_asc
  B_b_DR_asc_0 = 0,
  B_b_DR_asc_1 = 0,
  
  
  ### Quota_Count
  B_b_Quota_Count = 0,
  
  
  ### Quota_Max_Duration
  B_b_Quota_Max_Duration = 0,
  
  
  ### Quota_Period_midmorning
  B_b_Quota_Period_midmorning = 0,
  
  ### Quota_Period_afternoon
  B_b_Quota_Period_afternoon = 0,
  
  ### Quota_Period_evening
  B_b_Quota_Period_evening = 0,
  
  ### Quota_Period_night
  B_b_Quota_Period_night = 0,
  
  
  ### Quota_Compensation
  B_b_Quota_Compensation = -0.1,
  
  
  # ----------------------------------------------------------------- #
  #---- Intra-Class-Parameters (C)
  # ----------------------------------------------------------------- #
  ### DR_asc
  C_b_DR_asc_0 = 0,
  C_b_DR_asc_1 = 0,
  
  
  ### Quota_Count
  C_b_Quota_Count = 0,
  
  
  ### Quota_Max_Duration
  C_b_Quota_Max_Duration = 0,
  
  
  ### Quota_Period_midmorning
  C_b_Quota_Period_midmorning = 0,
  
  ### Quota_Period_afternoon
  C_b_Quota_Period_afternoon = 0,
  
  ### Quota_Period_evening
  C_b_Quota_Period_evening = 0,
  
  ### Quota_Period_night
  C_b_Quota_Period_night = 0,
  
  
  ### Quota_Compensation
  C_b_Quota_Compensation = -0.2,
  
  
  # ----------------------------------------------------------------- #
  #---- Latent Variable Model
  # ----------------------------------------------------------------- #
  ### Influence of LV on other LV
  lambda_NACountComplexityCosts_NACount = 1,
  lambda_NACountIrrelevance_NACount = 1,
  
  lambda_NADurationComplexityCosts_NADuration = 1,
  lambda_NADurationIrrelevance_NADuration = 1,
  
  lambda_NAPeriodComplexityCosts_NAPeriod = 1,
  lambda_NAPeriodIrrelevance_NAPeriod = 1,
  
  lambda_NACompensationComplexityCosts_NACompensation = 1,
  lambda_NACompensationIrrelevance_NACompensation = 1,
  
  ### Continuous measurement model
  sigma_NACount_r1 = 1,
  sigma_NACount_r2 = 1,
  sigma_NACount_r3 = 1,
  sigma_NACount_r4 = 1,
  sigma_NACount_r5 = 1,
  sigma_NACount_r6 = 1,
  
  sigma_NADuration_r1 = 1,
  sigma_NADuration_r2 = 1,
  sigma_NADuration_r3 = 1,
  sigma_NADuration_r4 = 1,
  sigma_NADuration_r5 = 1,
  sigma_NADuration_r6 = 1,
  
  sigma_NAPeriod_r1 = 1,
  sigma_NAPeriod_r2 = 1,
  sigma_NAPeriod_r3 = 1,
  sigma_NAPeriod_r4 = 1,
  sigma_NAPeriod_r5 = 1,
  sigma_NAPeriod_r6 = 1,
  
  sigma_NACompensation_r1 = 1,
  sigma_NACompensation_r2 = 1,
  sigma_NACompensation_r3 = 1,
  sigma_NACompensation_r4 = 1,
  sigma_NACompensation_r5 = 1,
  sigma_NACompensation_r6 = 1,
  
  ### Marketing scales
  zeta_NACountComplexityCosts_r1_NACount= 1,
  zeta_NACountComplexityCosts_r2_NACount = 1,
  zeta_NACountComplexityCosts_r3_NACount = 1,
  zeta_NACountComplexityCosts_r4_NACount = 1,
  zeta_NACountIrrelevance_r5_NACount = 1,
  zeta_NACountIrrelevance_r6_NACount = 1,
  
  zeta_NADurationComplexityCosts_r1_NADuration= 1,
  zeta_NADurationComplexityCosts_r2_NADuration = 1,
  zeta_NADurationComplexityCosts_r3_NADuration = 1,
  zeta_NADurationComplexityCosts_r4_NADuration = 1,
  zeta_NADurationIrrelevance_r5_NADuration = 1,
  zeta_NADurationIrrelevance_r6_NADuration = 1,
  
  zeta_NAPeriodComplexityCosts_r1_NAPeriod= 1,
  zeta_NAPeriodComplexityCosts_r2_NAPeriod = 1,
  zeta_NAPeriodComplexityCosts_r3_NAPeriod = 1,
  zeta_NAPeriodComplexityCosts_r4_NAPeriod = 1,
  zeta_NAPeriodIrrelevance_r5_NAPeriod = 1,
  zeta_NAPeriodIrrelevance_r6_NAPeriod = 1,
  
  zeta_NACompensationComplexityCosts_r1_NACompensation= 1,
  zeta_NACompensationComplexityCosts_r2_NACompensation = 1,
  zeta_NACompensationComplexityCosts_r3_NACompensation = 1,
  zeta_NACompensationComplexityCosts_r4_NACompensation = 1,
  zeta_NACompensationIrrelevance_r5_NACompensation = 1,
  zeta_NACompensationIrrelevance_r6_NACompensation = 1,
  
  ### Sociodemographics
  gamma_Age_NACountComplexityCosts = 0,
  gamma_Education_NACountComplexityCosts = 0,
  gamma_ExistingApps_NACountComplexityCosts = 0,
  sigma_NACountComplexityCosts = 1,
  
  gamma_Age_NACountIrrelevance = 0,
  gamma_Education_NACountIrrelevance = 0,
  gamma_ExistingApps_NACountIrrelevance = 0,
  sigma_NACountIrrelevance = 1,
  
  
  gamma_Age_NADurationComplexityCosts = 0,
  gamma_Education_NADurationComplexityCosts = 0,
  gamma_ExistingApps_NADurationComplexityCosts = 0,
  sigma_NADurationComplexityCosts = 1,
  
  gamma_Age_NADurationIrrelevance = 0,
  gamma_Education_NADurationIrrelevance = 0,
  gamma_ExistingApps_NADurationIrrelevance = 0,
  sigma_NADurationIrrelevance = 1,
  
  
  gamma_Age_NAPeriodComplexityCosts = 0,
  gamma_Education_NAPeriodComplexityCosts = 0,
  gamma_ExistingApps_NAPeriodComplexityCosts = 0,
  sigma_NAPeriodComplexityCosts = 1,
  
  gamma_Age_NAPeriodIrrelevance = 0,
  gamma_Education_NAPeriodIrrelevance = 0,
  gamma_ExistingApps_NAPeriodIrrelevance = 0,
  sigma_NAPeriodIrrelevance = 1,
  
  
  gamma_Gender_NACompensationComplexityCosts = 0,
  gamma_Age_NACompensationComplexityCosts = 0,
  gamma_Education_NACompensationComplexityCosts = 0,
  gamma_Income_NACompensationComplexityCosts = 0,
  gamma_CheapTalk_NACompensationComplexityCosts = 0,
  sigma_NACompensationComplexityCosts = 1,
  
  gamma_Gender_NACompensationIrrelevance = 0,
  gamma_Age_NACompensationIrrelevance = 0,
  gamma_Education_NACompensationIrrelevance = 0,
  gamma_Income_NACompensationIrrelevance = 0,
  gamma_CheapTalk_NACompensationIrrelevance = 0,
  sigma_NACompensationIrrelevance = 1
)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_bwtp_fixed = c() if none
apollo_fixed = c("b_asc_3",
                 
                 "delta_A",
                 
                 "gamma_A_NACount",
                 "gamma_A_NADuration",
                 "gamma_A_NAPeriod",
                 "gamma_A_NACompensation",
                 
                 "lambda_NACountComplexityCosts_NACount",
                 "lambda_NADurationComplexityCosts_NADuration",
                 "lambda_NAPeriodComplexityCosts_NAPeriod",
                 "lambda_NACompensationComplexityCosts_NACompensation",
                 
                 "zeta_NACountComplexityCosts_r1_NACount",  # fix most reliable item if possible/known
                 "zeta_NACountIrrelevance_r5_NACount",
                 
                 "zeta_NADurationComplexityCosts_r1_NADuration",  # fix most reliable item if possible/known
                 "zeta_NADurationIrrelevance_r5_NADuration",
                 
                 "zeta_NAPeriodComplexityCosts_r1_NAPeriod",  # fix most reliable item if possible/known
                 "zeta_NAPeriodIrrelevance_r5_NAPeriod",
                 
                 "zeta_NACompensationComplexityCosts_r1_NACompensation",  # fix most reliable item if possible/known
                 "zeta_NACompensationIrrelevance_r5_NACompensation",
                 
                 "A_b_DR_asc_0",
                 "B_b_DR_asc_0",
                 "C_b_DR_asc_0"
                 )

### Read estimates from other models
apollo_beta = apollo_readBeta(apollo_beta, apollo_fixed, "flexQgrid_LC_3Classes_DR1_CutOff6", overwriteFixed=FALSE)
# apollo_beta_new = replace(apollo_beta, names(priors.step1$FC), priors.step1$FC)
# apollo_beta_new = replace(apollo_beta_new, names(priors.step1$svN), priors.step1$svN)
# apollo_beta = apollo_beta_new[names(apollo_beta)]

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

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "sobolOwen",
  interNDraws    = 2000,
  interUnifDraws = c(),
  interNormDraws = c("eta_NACountComplexityCosts", "eta_NACountIrrelevance",
                     "eta_NADurationComplexityCosts", "eta_NADurationIrrelevance",
                     "eta_NAPeriodComplexityCosts", "eta_NAPeriodIrrelevance",
                     "eta_NACompensationComplexityCosts", "eta_NACompensationIrrelevance"),
  intraDrawsType = "sobolOwen",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs) {
  randcoeff = list()
  
  # ----------------------------------------------------------------- #
  #---- Latent Variables
  # ----------------------------------------------------------------- #
  randcoeff[["LV_NACountComplexityCosts"]] = gamma_Age_NACountComplexityCosts * COV_Age_centered +
    gamma_Education_NACountComplexityCosts * COV_Education_log_centered +
    gamma_ExistingApps_NACountComplexityCosts * COV_ExistingApps +
    sigma_NACountComplexityCosts * eta_NACountComplexityCosts
  
  randcoeff[["LV_NACountIrrelevance"]] = gamma_Age_NACountIrrelevance * COV_Age_centered +
    gamma_Education_NACountIrrelevance * COV_Education_log_centered +
    gamma_ExistingApps_NACountIrrelevance * COV_ExistingApps +
    sigma_NACountIrrelevance * eta_NACountIrrelevance
  
  randcoeff[["LV_NACount"]] = lambda_NACountComplexityCosts_NACount *
    ( # LV_NACountComplexityCosts
      gamma_Age_NACountComplexityCosts * COV_Age_centered + 
      gamma_Education_NACountComplexityCosts * COV_Education_log_centered +
      gamma_ExistingApps_NACountComplexityCosts * COV_ExistingApps +
      sigma_NACountComplexityCosts * eta_NACountComplexityCosts
    ) +
    lambda_NACountIrrelevance_NACount *
    ( # LV_NACountIrrelevance
      gamma_Age_NACountIrrelevance * COV_Age_centered +
      gamma_Education_NACountIrrelevance * COV_Education_log_centered +
      gamma_ExistingApps_NACountIrrelevance * COV_ExistingApps +
      sigma_NACountIrrelevance * eta_NACountIrrelevance
    )
  
  
  randcoeff[["LV_NADurationComplexityCosts"]] = gamma_Age_NADurationComplexityCosts * COV_Age_centered +
    gamma_Education_NADurationComplexityCosts * COV_Education_log_centered +
    gamma_ExistingApps_NADurationComplexityCosts * COV_ExistingApps +
    sigma_NADurationComplexityCosts * eta_NADurationComplexityCosts
  
  randcoeff[["LV_NADurationIrrelevance"]] = gamma_Age_NADurationIrrelevance * COV_Age_centered +
    gamma_Education_NADurationIrrelevance * COV_Education_log_centered +
    gamma_ExistingApps_NADurationIrrelevance * COV_ExistingApps +
    sigma_NADurationIrrelevance * eta_NADurationIrrelevance
  
  randcoeff[["LV_NADuration"]] = lambda_NADurationComplexityCosts_NADuration *
    ( # LV_NADurationComplexityCosts
      gamma_Age_NADurationComplexityCosts * COV_Age_centered +
      gamma_Education_NADurationComplexityCosts * COV_Education_log_centered +
      gamma_ExistingApps_NADurationComplexityCosts * COV_ExistingApps +
      sigma_NADurationComplexityCosts * eta_NADurationComplexityCosts
    ) +
    lambda_NADurationIrrelevance_NADuration *
    ( # LV_NADurationIrrelevance
      gamma_Age_NADurationIrrelevance * COV_Age_centered +
      gamma_Education_NADurationIrrelevance * COV_Education_log_centered +
      gamma_ExistingApps_NADurationIrrelevance * COV_ExistingApps +
      sigma_NADurationIrrelevance * eta_NADurationIrrelevance
    )
  
  
  randcoeff[["LV_NAPeriodComplexityCosts"]] = gamma_Age_NAPeriodComplexityCosts * COV_Age_centered +
    gamma_Education_NAPeriodComplexityCosts * COV_Education_log_centered +
    gamma_ExistingApps_NAPeriodComplexityCosts * COV_ExistingApps +
    sigma_NAPeriodComplexityCosts * eta_NAPeriodComplexityCosts
  
  randcoeff[["LV_NAPeriodIrrelevance"]] = gamma_Age_NAPeriodIrrelevance * COV_Age_centered +
    gamma_Education_NAPeriodIrrelevance * COV_Education_log_centered +
    gamma_ExistingApps_NAPeriodIrrelevance * COV_ExistingApps +
    sigma_NAPeriodIrrelevance * eta_NAPeriodIrrelevance
  
  randcoeff[["LV_NAPeriod"]] = lambda_NAPeriodComplexityCosts_NAPeriod *
    ( # LV_NAPeriodComplexityCosts
      gamma_Age_NAPeriodComplexityCosts * COV_Age_centered +
      gamma_Education_NAPeriodComplexityCosts * COV_Education_log_centered +
      gamma_ExistingApps_NAPeriodComplexityCosts * COV_ExistingApps +
      sigma_NAPeriodComplexityCosts * eta_NAPeriodComplexityCosts
    ) +
    lambda_NAPeriodIrrelevance_NAPeriod *
    ( # LV_NAPeriodIrrelevance
      gamma_Age_NAPeriodIrrelevance * COV_Age_centered +
      gamma_Education_NAPeriodIrrelevance * COV_Education_log_centered +
      gamma_ExistingApps_NAPeriodIrrelevance * COV_ExistingApps +
      sigma_NAPeriodIrrelevance * eta_NAPeriodIrrelevance
    )
  
  
  randcoeff[["LV_NACompensationComplexityCosts"]] = gamma_Gender_NACompensationComplexityCosts * COV_Gender +
    gamma_Age_NACompensationComplexityCosts * COV_Age_centered +
    gamma_Education_NACompensationComplexityCosts * COV_Education_log_centered +
    (1-IncomeQuan_miss) * gamma_Income_NACompensationComplexityCosts * COV_IncomeClasses_log_centered +
    gamma_CheapTalk_NACompensationComplexityCosts * COV_CheapTalk +
    sigma_NACompensationComplexityCosts * eta_NACompensationComplexityCosts
  
  randcoeff[["LV_NACompensationIrrelevance"]] = gamma_Gender_NACompensationIrrelevance * COV_Gender +
    gamma_Age_NACompensationIrrelevance * COV_Age_centered +
    gamma_Education_NACompensationIrrelevance * COV_Education_log_centered +
    (1-IncomeQuan_miss) * gamma_Income_NACompensationIrrelevance * COV_IncomeClasses_log_centered +
    gamma_CheapTalk_NACompensationIrrelevance * COV_CheapTalk +
    sigma_NACompensationIrrelevance * eta_NACompensationIrrelevance
  
  randcoeff[["LV_NACompensation"]] = lambda_NACompensationComplexityCosts_NACompensation *
    ( # LV_NACompensationComplexityCosts
      gamma_Gender_NACompensationComplexityCosts * COV_Gender +
      gamma_Age_NACompensationComplexityCosts * COV_Age_centered +
      gamma_Education_NACompensationComplexityCosts * COV_Education_log_centered +
      (1-IncomeQuan_miss) * gamma_Income_NACompensationComplexityCosts * COV_IncomeClasses_log_centered +
      gamma_CheapTalk_NACompensationComplexityCosts * COV_CheapTalk +
      sigma_NACompensationComplexityCosts * eta_NACompensationComplexityCosts
    ) +
    lambda_NACompensationIrrelevance_NACompensation *
    ( # LV_NACompensationIrrelevance
      gamma_Gender_NACompensationIrrelevance * COV_Gender +
      gamma_Age_NACompensationIrrelevance * COV_Age_centered +
      gamma_Education_NACompensationIrrelevance * COV_Education_log_centered +
      (1-IncomeQuan_miss) * gamma_Income_NACompensationIrrelevance * COV_IncomeClasses_log_centered +
      gamma_CheapTalk_NACompensationIrrelevance * COV_CheapTalk +
      sigma_NACompensationIrrelevance * eta_NACompensationIrrelevance
    )
  
  return(randcoeff)
}
  

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

apollo_lcPars = function(apollo_beta, apollo_inputs){
  lcpars = list()
  
  lcpars[["b_DR_asc_0"]] = list(A_b_DR_asc_0, B_b_DR_asc_0, C_b_DR_asc_0)
  lcpars[["b_DR_asc_1"]] = list(A_b_DR_asc_1, B_b_DR_asc_1, C_b_DR_asc_1)
  
  lcpars[["b_Quota_Count"]] = list(A_b_Quota_Count, B_b_Quota_Count, C_b_Quota_Count)
  
  lcpars[["b_Quota_Max_Duration"]] = list(A_b_Quota_Max_Duration, B_b_Quota_Max_Duration, C_b_Quota_Max_Duration)
  
  lcpars[["b_Quota_Period_midmorning"]] = list(A_b_Quota_Period_midmorning, B_b_Quota_Period_midmorning, C_b_Quota_Period_midmorning)
  lcpars[["b_Quota_Period_afternoon"]] = list(A_b_Quota_Period_afternoon, B_b_Quota_Period_afternoon, C_b_Quota_Period_afternoon)
  lcpars[["b_Quota_Period_evening"]] = list(A_b_Quota_Period_evening, B_b_Quota_Period_evening, C_b_Quota_Period_evening)
  lcpars[["b_Quota_Period_night"]] = list(A_b_Quota_Period_night, B_b_Quota_Period_night, C_b_Quota_Period_night)
  
  lcpars[["b_Quota_Compensation"]] = list(A_b_Quota_Compensation, B_b_Quota_Compensation, C_b_Quota_Compensation)
  
  V=list()
  V[["class_A"]] = delta_A +
    gamma_A_NACount * LV_NACount +
    gamma_A_NADuration * LV_NADuration +
    gamma_A_NAPeriod * LV_NAPeriod +
    gamma_A_NACompensation * LV_NACompensation
  
  V[["class_B"]] = delta_B +
    gamma_B_NACount * LV_NACount +
    gamma_B_NADuration * LV_NADuration +
    gamma_B_NAPeriod * LV_NAPeriod +
    gamma_B_NACompensation * LV_NACompensation
  
  V[["class_C"]] = delta_C +
    gamma_C_NACount * LV_NACount +
    gamma_C_NADuration * LV_NADuration +
    gamma_C_NAPeriod * LV_NAPeriod +
    gamma_C_NACompensation * LV_NACompensation
  
  mnl_settings = list(
    alternatives = c(class_A=1, class_B=2, class_C=3), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  
  lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
  lcpars[["pi_values"]] = apollo_firstRow(lcpars[["pi_values"]], apollo_inputs)
  
  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()
  
  # ----------------------------------------------------------------- #
  #---- Likelihood of indicators
  # ----------------------------------------------------------------- #
  ### NACount
  normalDensity_settings1 = list(outcomeNormal = NonAttendance_r1, 
                                 xNormal       = zeta_NACountComplexityCosts_r1_NACount * LV_NACountComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r1, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r1_NACount")
  
  normalDensity_settings2 = list(outcomeNormal = NonAttendance_r2, 
                                 xNormal       = zeta_NACountComplexityCosts_r2_NACount * LV_NACountComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r2, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r2_NACount")
  
  normalDensity_settings3 = list(outcomeNormal = NonAttendance_r3, 
                                 xNormal       = zeta_NACountComplexityCosts_r3_NACount * LV_NACountComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r3, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r3_NACount")
  
  normalDensity_settings4 = list(outcomeNormal = NonAttendance_r4, 
                                 xNormal       = zeta_NACountComplexityCosts_r4_NACount * LV_NACountComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r4, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r4_NACount")
  
  normalDensity_settings5 = list(outcomeNormal = NonAttendance_r5, 
                                 xNormal       = zeta_NACountIrrelevance_r5_NACount * LV_NACountIrrelevance, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r5, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r5_NACount")
  
  normalDensity_settings6 = list(outcomeNormal = NonAttendance_r6, 
                                 xNormal       = zeta_NACountIrrelevance_r6_NACount * LV_NACountIrrelevance, 
                                 mu            = 0, 
                                 sigma         = sigma_NACount_r6, 
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NACount_r6_NACount")
  
  ###  NADuration
  normalDensity_settings7 = list(outcomeNormal = NonAttendance2_r1,
                                 xNormal       = zeta_NADurationComplexityCosts_r1_NADuration * LV_NADurationComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NADuration_r1,
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NADuration_r1_NADuration")
  
  normalDensity_settings8 = list(outcomeNormal = NonAttendance2_r2,
                                 xNormal       = zeta_NADurationComplexityCosts_r2_NADuration * LV_NADurationComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NADuration_r2,
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NADuration_r2_NADuration")
  
  normalDensity_settings9 = list(outcomeNormal = NonAttendance2_r3,
                                 xNormal       = zeta_NADurationComplexityCosts_r3_NADuration * LV_NADurationComplexityCosts, 
                                 mu            = 0, 
                                 sigma         = sigma_NADuration_r3,
                                 rows=(sys_RespNum.first==1),
                                 componentName = "NADuration_r3_NADuration")
  
  normalDensity_settings10 = list(outcomeNormal = NonAttendance2_r4,
                                  xNormal       = zeta_NADurationComplexityCosts_r4_NADuration * LV_NADurationComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NADuration_r4,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NADuration_r4_NADuration")
  
  normalDensity_settings11 = list(outcomeNormal = NonAttendance2_r5,
                                  xNormal       = zeta_NADurationIrrelevance_r5_NADuration * LV_NADurationIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NADuration_r5,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NADuration_r5_NADuration")
  
  normalDensity_settings12 = list(outcomeNormal = NonAttendance2_r6,
                                  xNormal       = zeta_NADurationIrrelevance_r6_NADuration * LV_NADurationIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NADuration_r6,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NADuration_r6_NADuration")
  
  ### NAPeriod
  normalDensity_settings13 = list(outcomeNormal = NonAttendance3_r1,
                                  xNormal       = zeta_NAPeriodComplexityCosts_r1_NAPeriod * LV_NAPeriodComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r1,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r1_NAPeriod")
  
  normalDensity_settings14 = list(outcomeNormal = NonAttendance3_r2,
                                  xNormal       = zeta_NAPeriodComplexityCosts_r2_NAPeriod * LV_NAPeriodComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r2,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r2_NAPeriod")
  
  normalDensity_settings15 = list(outcomeNormal = NonAttendance3_r3,
                                  xNormal       = zeta_NAPeriodComplexityCosts_r3_NAPeriod * LV_NAPeriodComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r3,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r3_NAPeriod")
  
  normalDensity_settings16 = list(outcomeNormal = NonAttendance3_r4,
                                  xNormal       = zeta_NAPeriodComplexityCosts_r4_NAPeriod * LV_NAPeriodComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r4,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r4_NAPeriod")
  
  normalDensity_settings17 = list(outcomeNormal = NonAttendance3_r5,
                                  xNormal       = zeta_NAPeriodIrrelevance_r5_NAPeriod * LV_NAPeriodIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r5,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r5_NAPeriod")
  
  normalDensity_settings18 = list(outcomeNormal = NonAttendance3_r6,
                                  xNormal       = zeta_NAPeriodIrrelevance_r6_NAPeriod * LV_NAPeriodIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NAPeriod_r6,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NAPeriod_r6_NAPeriod")
  
  ### NACompensation
  normalDensity_settings19 = list(outcomeNormal = NonAttendance4_r1,
                                  xNormal       = zeta_NACompensationComplexityCosts_r1_NACompensation * LV_NACompensationComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r1,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r1_NACompensation")
  
  normalDensity_settings20 = list(outcomeNormal = NonAttendance4_r2,
                                  xNormal       = zeta_NACompensationComplexityCosts_r2_NACompensation * LV_NACompensationComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r2,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r2_NACompensation")
  
  normalDensity_settings21 = list(outcomeNormal = NonAttendance4_r3,
                                  xNormal       = zeta_NACompensationComplexityCosts_r3_NACompensation * LV_NACompensationComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r3,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r3_NACompensation")
  
  normalDensity_settings22 = list(outcomeNormal = NonAttendance4_r4,
                                  xNormal       = zeta_NACompensationComplexityCosts_r4_NACompensation * LV_NACompensationComplexityCosts, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r4,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r4_NACompensation")
  
  normalDensity_settings23 = list(outcomeNormal = NonAttendance4_r5,
                                  xNormal       = zeta_NACompensationIrrelevance_r5_NACompensation * LV_NACompensationIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r5,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r5_NACompensation")
  
  normalDensity_settings24 = list(outcomeNormal = NonAttendance4_r6,
                                  xNormal       = zeta_NACompensationIrrelevance_r6_NACompensation * LV_NACompensationIrrelevance, 
                                  mu            = 0, 
                                  sigma         = sigma_NACompensation_r6,
                                  rows=(sys_RespNum.first==1),
                                  componentName = "NACompensation_r6_NACompensation")
  
  ### NACount
  P[["indic_NACount_NACount_r1"]] = apollo_normalDensity(normalDensity_settings1, functionality)
  P[["indic_NACount_NACount_r2"]] = apollo_normalDensity(normalDensity_settings2, functionality)
  P[["indic_NACount_NACount_r3"]] = apollo_normalDensity(normalDensity_settings3, functionality)
  P[["indic_NACount_NACount_r4"]] = apollo_normalDensity(normalDensity_settings4, functionality)
  P[["indic_NACount_NACount_r5"]] = apollo_normalDensity(normalDensity_settings5, functionality)
  P[["indic_NACount_NACount_r6"]] = apollo_normalDensity(normalDensity_settings6, functionality)
  
  ### NADuration
  P[["indic_NADuration_NADuration_r1"]] = apollo_normalDensity(normalDensity_settings7, functionality)
  P[["indic_NADuration_NADuration_r2"]] = apollo_normalDensity(normalDensity_settings8, functionality)
  P[["indic_NADuration_NADuration_r3"]] = apollo_normalDensity(normalDensity_settings9, functionality)
  P[["indic_NADuration_NADuration_r4"]] = apollo_normalDensity(normalDensity_settings10, functionality)
  P[["indic_NADuration_NADuration_r5"]] = apollo_normalDensity(normalDensity_settings11, functionality)
  P[["indic_NADuration_NADuration_r6"]] = apollo_normalDensity(normalDensity_settings12, functionality)
  
  ### NAPeriod
  P[["indic_NAPeriod_NAPeriod_r1"]] = apollo_normalDensity(normalDensity_settings13, functionality)
  P[["indic_NAPeriod_NAPeriod_r2"]] = apollo_normalDensity(normalDensity_settings14, functionality)
  P[["indic_NAPeriod_NAPeriod_r3"]] = apollo_normalDensity(normalDensity_settings15, functionality)
  P[["indic_NAPeriod_NAPeriod_r4"]] = apollo_normalDensity(normalDensity_settings16, functionality)
  P[["indic_NAPeriod_NAPeriod_r5"]] = apollo_normalDensity(normalDensity_settings17, functionality)
  P[["indic_NAPeriod_NAPeriod_r6"]] = apollo_normalDensity(normalDensity_settings18, functionality)
  
  
  ### NACompensation
  P[["indic_NACompensation_NACompensation_r1"]] = apollo_normalDensity(normalDensity_settings19, functionality)
  P[["indic_NACompensation_NACompensation_r2"]] = apollo_normalDensity(normalDensity_settings20, functionality)
  P[["indic_NACompensation_NACompensation_r3"]] = apollo_normalDensity(normalDensity_settings21, functionality)
  P[["indic_NACompensation_NACompensation_r4"]] = apollo_normalDensity(normalDensity_settings22, functionality)
  P[["indic_NACompensation_NACompensation_r5"]] = apollo_normalDensity(normalDensity_settings23, functionality)
  P[["indic_NACompensation_NACompensation_r6"]] = apollo_normalDensity(normalDensity_settings24, functionality)
  
  ### Combined model
  P = apollo_combineModels(P, apollo_inputs, 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)
  
  ### Rename model
  names(P)[which(names(P)=="model")] <- "Measurement_model"
  
  # ----------------------------------------------------------------- #
  #---- LC model
  # ----------------------------------------------------------------- #
  ### Create list of probabilities P
  P_within <- list() # list for within class models
  
  ### Loop over classes
  for (s in 1:3) {
    
    ### Create list of probabilities P.componentName
    P.componentName <- list()
    
    ### Compute class-specific utilities
    
    ### CBC (Forced Choices)
    V = list()
    
    V[['alt1']] = b_asc_1 +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    V[['alt2']] = b_asc_2 +
      b_Quota_Count[[s]] * Quota_Count.2 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.2 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
      b_Quota_Compensation[[s]] * Quota_Compensation.2
    
    V[['alt3']] = b_asc_3 +
      b_Quota_Count[[s]] * Quota_Count.3 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.3 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
      b_Quota_Compensation[[s]] * Quota_Compensation.3
    
    componentName1 = paste("CBC_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt1=1, alt2=2, alt3=3),
      avail         = list(alt1=1, alt2=1, alt3=1),
      choiceVar     = Choice,
      V             = V,
      rows          = (CBC_type=="R"),
      componentName = componentName1
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName1]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName1]] = apollo_panelProd(P.componentName[[componentName1]], apollo_inputs, functionality)
    
    
    ### CBC (Free Choices)
    V = list()
    
    V[['alt0']] = b_DR_asc_0[[s]]
    
    V[['alt1']] = b_DR_asc_1[[s]] +
      b_Quota_Count[[s]] * Quota_Count.1 +
      b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
      b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
      b_Quota_Compensation[[s]] * Quota_Compensation.1
    
    componentName2 = paste("CBC_DR_", toString(s), sep="")
    
    mnl_settings = list(
      alternatives  = c(alt0=0, alt1=1),
      avail         = list(alt0=1, alt1=1),
      choiceVar     = Choice,
      V             = V,
      rows          = (CBC_type=="RD"),
      componentName = componentName2
    )
    
    ### Compute within-class choice probabilities using MNL model
    P.componentName[[componentName2]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P.componentName[[componentName2]] = apollo_panelProd(P.componentName[[componentName2]], apollo_inputs, functionality)
    
    ### Combined model
    P_within[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality, asList=FALSE)
  }
  
  ### Compute latent class model probabilities
  lc_settings = list(inClassProb = P_within, classProb=pi_values)
  P[["LC_model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
  
  # ----------------------------------------------------------------- #
  #---- Combine model components
  # ----------------------------------------------------------------- #
  ### Combined model
  P = apollo_combineModels(P, apollo_inputs, 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)
}

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

estimate_settings = list(hessianRoutine = "maxLik",
                         estimationRoutine = "bfgs",
                         maxIterations = 300,
                         printLevel = 3)

model = apollo_estimate(apollo_beta,
                        apollo_fixed,
                        apollo_probabilities,
                        apollo_inputs,
                        estimate_settings = estimate_settings)

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

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

apollo_modelOutput(model)

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

saveOutput_settings = list(
  printClassical = TRUE,
  printPVal = TRUE,
  printDiagnostics = TRUE,
  printCovar = TRUE,
  printCorr = TRUE,
  printOutliers = 100,
  printChange = TRUE,
  saveEst = TRUE,
  saveCov = TRUE,
  saveCorr = TRUE,
  saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)

As always: Thanks in advance for the help!

Nico

Re: Latent Class Joint Model

Posted: 10 Jun 2021, 10:27

by stephanehess

Nico

thanks. So what happens if you run this?

Stephane