Errors when estimating a latent variable model

[1123088145@qq.com]

Dear Stephane and David,

Thank you very much for inventing the Apollo package and creating the forum allowing our users to communicate.
Recently, I was trying to run an mnl model containing latent variables, but the code encountered the following error.
Overview of choices for model component indic_Benefit3:
1 2 3 4 5
Times chosen 24.00 176.00 608.00 1456.00 1448.00
Percentage chosen overall 0.65 4.74 16.38 39.22 39.01

Overview of choices for MNL model component choice:
V_bike1 V_bike2 V_bike3 V_bus V_taxi V_sharecar
Times available 3712.00 3712.00 3712.00 3712.00 3712.00 3712.00
Times chosen 587.00 814.00 117.00 545.00 187.00 235.00
Percentage chosen overall 15.81 21.93 3.15 14.68 5.04 6.33
Percentage chosen when available 15.81 21.93 3.15 14.68 5.04 6.33
V_car V_walk
Times available 3712.00 3712.00
Times chosen 655.00 572.00
Percentage chosen overall 17.65 15.41
Percentage chosen when available 17.65 15.41
Error in apollo_avgInterDraws(P, apollo_inputs, functionality) :
No Inter-individuals draws to average over!
My question is
Why is this situation happening?

I have attached my R code below.
Thank you very much.
Ren.

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       = "Hybrid_with_OL",
  modelDescr      = "Hybrid choice model on drug choice data, using ordered measurement model for indicators",
  indivID         = "id",
  mixing          = TRUE,
  nCores          = 4, 
  outputDirectory = "output"
)

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

### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv), 
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("huizong.csv",header=TRUE)
#database = apollo_drugChoiceData
### for data dictionary, use ?apollo_drugChoiceData

# ################################################################# #
#### ANALYSIS OF CHOICES                                         ####
# ################################################################# #

### Illustration of how to use apollo_choiceAnalysis with user-defined alternatives.
### This is useful in cases where the alternatives in the data differ 
### across tasks. The same approach can then also be used with unlabelled data



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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
  
  
  lambda   = 1,
  lambda_at1=0,
  lambda_at2=0,
  lambda_at3=0,
  gamma_occupation =0,
  gamma_age_35_less =0,
  gamma_male = 0,
  gamma_highEducation =0,
  gamma_high_income = 0,
  gamma_more_car = 0,
  gamma_no_car =  0,
  gamma_more_bike =0, 
  gamma_no_bike = 0,
  
  gamma_occupation1 =0,
  gamma_age_35_less1 =0,
  gamma_male1 = 0,
  gamma_highEducation1 =0,
  gamma_high_income1 = 0,
  gamma_more_car1 = 0,
  gamma_no_car1 =  0,
  gamma_more_bike1 =0, 
  gamma_no_bike1 = 0,
  
  gamma_occupation2 =0,
  gamma_age_35_less2 =0,
  gamma_male2 = 0,
  gamma_highEducation2 =0,
  gamma_high_income2 = 0,
  gamma_more_car2 = 0,
  gamma_no_car2 =  0,
  gamma_more_bike2 =0, 
  gamma_no_bike2 = 0,
  zeta_Availability1  = 1, 
  zeta_Tangibles1    = 1, 
  zeta_Benefit1       = 1, 
  zeta_Availability2  = 1, 
  zeta_Tangibles2     = 1, 
  zeta_Benefit2       = 1, 
  zeta_Availability3  = 1, 
  zeta_Tangibles3     = 1, 
  zeta_Benefit3       = 1, 
  zeta_Availability4  = 1, 
  zeta_Tangibles4     = 1, 
  tau_Availability_1  =-2, 
  tau_Availability_2   =-1, 
  tau_Availability_3   = 1, 
  tau_Availability_4    = 2, 
  
  tau_Tangibles_1  =-2, 
  tau_Tangibles_2  =-1, 
  tau_Tangibles_3  = 1, 
  tau_Tangibles_4  = 2, 
  tau_Benefit_1       =-2, 
  tau_Benefit_2       =-1, 
  tau_Benefit_3       = 1, 
  tau_Benefit_4       = 2, 
  
  ASC_bike1 = 0,
  ASC_bike2 = 0,
  ASC_bike3 = 0,
  ASC_bus = 0,
  ASC_taxi = 1,
  ASC_sharecar = 0,
  ASC_car = 0,
  ASC_walk =0 ,
  
  Beta_travel_time= 0,
  Beta_travel_time1= 0,
  Beta_travel_time2= 0,
  Beta_travel_time3= 0,
  Beta_travel_time4= 0,
  Beta_travel_time5= 0,
  Beta_travel_time6= 0,
  Beta_travel_time7= 0,
  Beta_travel_purpose= 0,
  Beta_travel_purpose1= 0,
  Beta_travel_purpose2= 0,
  Beta_travel_purpose3=0,
  Beta_travel_purpose4=0,
  Beta_travel_purpose5= 0,
  Beta_travel_purpose6= 0,
  Beta_travel_purpose7= 0,
  Beta_travel_purpose_income= 0,
  Beta_travel_purpose_income1= 0,
  Beta_travel_purpose_income2= 0,
  Beta_travel_purpose_income3=0,
  Beta_travel_purpose_income4=0,
  Beta_travel_purpose_income5= 0,
  Beta_travel_purpose_income6= 0,
  Beta_travel_purpose_income7= 0,
  
  Beta_travel_purpose_education= 0,
  Beta_travel_purpose_education1= 0,
  Beta_travel_purpose_education2= 0,
  Beta_travel_purpose_education3=0,
  Beta_travel_purpose_education4=0,
  Beta_travel_purpose_education5= 0,
  Beta_travel_purpose_education6= 0,
  Beta_travel_purpose_education7= 0,
  Beta_weather= 0,
  Beta_weather1=0,
  Beta_weather2= 0,
  Beta_weather3= 0,
  Beta_weather4=0,
  Beta_weather5= 0,
  Beta_weather6= 0,
  Beta_weather7= 0,
  Beta_weather_male= 0,
  Beta_weather_male1=0,
  Beta_weather_male2= 0,
  Beta_weather_male3= 0,
  Beta_weather_male4=0,
  Beta_weather_male5= 0,
  Beta_weather_male6= 0,
  Beta_weather_male7= 0,
  Beta_bike_waittime=0,
  Beta_bike_intime= 0,
  Beta_bike1_travel_time=0,
  Beta_bike1_cost= 0,
  Beta_bike2_walkt_time=0,
  Beta_bike2_intime= 0,
  Beta_bike2_travel_time= 0,
  Beta_bike2_cost= 0,
  Beta_bike3_walktime= 0,
  Beta_bike3_intime= 0,
  Beta_bike3_travel_time=0,
  Beta_bike3_cost=0,
  Beta_bus_waittime= 0,
  Beta_bus_intime= 0,
  Beta_bus_traveltime= 0,
  Beta_bus_cost= 0,
  Beta_Taxi_waittime= 0,
  Beta_Taxi_intime= 0,
  Beta_Taxi_traveltime=0,
  Beta_Taxi_cost= 0,
  Beta_share_car_waittime= 0,
  Beta_share_car_intime= 0,
  Beta_share_car_travel_time=0,
  Beta_share_car_cost=0,
  Beta_carpostcar=0,
  Beta_car_traveltime= 0,
  Beta_walk_time=0
)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c("ASC_walk")

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

### Set parameters for generating draws
apollo_draws = list(
  ##
  interDrawsType="mlhs", 
  interNDraws=100,          
  interUnifDraws=c(),      
  interNormDraws=c("eta","eta1","eta2") ,
  ## 
  intraDrawsType="",
  intraNDraws=0,          
  intraUnifDraws=c(),     
  intraNormDraws=c() 
)

### Create random parameters  结构方程
apollo_randCoeff=function(apollo_beta, apollo_inputs){
  randcoeff = list()
  randcoeff[["LV"]] =   gamma_high_income * ( income == 0 )
  + gamma_highEducation * (edu == 0)
  + gamma_occupation * (occupation == 3)
  + gamma_age_35_less * (age == 1)
  + gamma_male *  (gender == 1)
  + gamma_more_car* (own_car>=2)
  + gamma_no_car*  (own_car==0)
  + gamma_more_bike* (own_bike >= 2)
  + gamma_no_bike* (own_bike == 1)
  + eta
  randcoeff[["LV1"]] =   gamma_high_income1 * ( income == 0 )
  + gamma_highEducation1 * (edu == 0)
  + gamma_occupation1 * (occupation == 3)
  + gamma_age_35_less1 * (age == 1)
  + gamma_male1 *  (gender == 1)
  + gamma_more_car1* (own_car>=2)
  + gamma_no_car1*  (own_car==0)
  + gamma_more_bike1* (own_bike >= 2)
  + gamma_no_bike1* (own_bike == 1)
  + eta1
  randcoeff[["LV2"]] =   gamma_high_income2 * ( income == 0 )
  + gamma_highEducation2 * (edu == 0)
  + gamma_occupation2 * (occupation == 3)
  + gamma_age_35_less2 * (age == 1)
  + gamma_male2 *  (gender == 1)
  + gamma_more_car2* (own_car>=2)
  + gamma_no_car2*  (own_car==0)
  + gamma_more_bike2* (own_bike >= 2)
  + gamma_no_bike2* (own_bike == 1)
  + eta2
  return(randcoeff)
}

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

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION      测量方程                  ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P = list()
  
  ### Likelihood of indicators 
  ## 
  ol_settings1 = list(outcomeOrdered = A3, 
                      V              = zeta_Availability1*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4),
                      componentName  = "indic_Availability1")
  ol_settings2 = list(outcomeOrdered = A2, 
                      V              = zeta_Availability2*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4),
                      componentName  = "indic_Availability2")
  ol_settings3 = list(outcomeOrdered = A21, 
                      V              = zeta_Availability3*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4),
                      
                      componentName  = "indic_Availability3")
  ol_settings4 = list(outcomeOrdered = A16, 
                      V              = zeta_Availability4*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4),
                      componentName  = "indic_Availability4")
  ol_settings5 = list(outcomeOrdered = A10, 
                      V              = zeta_Tangibles1*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4), 
                      
                      componentName  = "indic_Tangibles1")
  ol_settings6 = list(outcomeOrdered = A11, 
                      V              = zeta_Tangibles2*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4), 
                      componentName  = "indic_Tangibles2")
  ol_settings7 = list(outcomeOrdered = A4, 
                      V              = zeta_Tangibles3*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4), 
                      componentName  = "indic_Tangibles3")
  ol_settings8 = list(outcomeOrdered = A17, 
                      V              = zeta_Tangibles4*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4), 
                      componentName  = "indic_Tangibles4")
  ol_settings9 = list(outcomeOrdered = A22, 
                      V              = zeta_Benefit1*LV2, 
                      tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4), 
                      componentName  = "indic_Benefit1")
  ol_settings10 = list(outcomeOrdered = A24, 
                       V              = zeta_Benefit2*LV2, 
                       tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4), 
                       componentName  = "indic_Benefit2")
  ol_settings11 = list(outcomeOrdered = A18, 
                       V              = zeta_Benefit3*LV2, 
                       tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4),  
                       componentName  = "indic_Benefit3")
  P[["indic_Availability1"]]= apollo_ol(ol_settings1, functionality)
  P[["indic_Availability2"]]= apollo_ol(ol_settings2, functionality)
  P[["indic_Availability3"]]= apollo_ol(ol_settings3, functionality)
  P[["indic_Availability4"]] = apollo_ol(ol_settings4, functionality)
  P[["indic_Tangibles1"]] = apollo_ol(ol_settings5, functionality)
  P[["indic_Tangibles2"]] = apollo_ol(ol_settings6, functionality)
  P[["indic_Tangibles3"]] = apollo_ol(ol_settings7, functionality)
  P[["indic_Tangibles4"]] = apollo_ol(ol_settings8, functionality)
  P[["indic_Benefit1"]]   = apollo_ol(ol_settings9, functionality)
  P[["indic_Benefit2"]]   = apollo_ol(ol_settings10, functionality)
  P[["indic_Benefit3"]]   = apollo_ol(ol_settings11, functionality)
  
  ### Likelihood of choices
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["V_bike1"]]= (ASC_bike1 + Beta_travel_purpose * (travel_purpose == 1)
                   + Beta_weather * (weather == 1)
                   + (travel_time == 1) * Beta_travel_time
                   + Beta_bike_waittime * bike_waittime
                   + Beta_bike_intime * bike_intime
                   + Beta_bike1_travel_time * bike1_travel_time
                   + Beta_bike1_cost * bike1_cost
                   + Beta_travel_purpose_income * (travel_purpose == 1) * (income == 1)
                   + Beta_weather_male * (weather == 1) * (gender == 1)
                   + lambda_at1*LV+lambda_at2*LV1 + lambda_at3 * LV2
  )
  
  
  V[["V_bike2"]] = (ASC_bike2 + Beta_travel_purpose1 * (travel_purpose == 1)
                    + Beta_weather1 * (weather == 1)
                    + (travel_time == 1) * Beta_travel_time1
                    + Beta_bike2_walkt_time * bike2_walkt_time
                    + Beta_bike2_intime * bike2_intime
                    + Beta_bike2_travel_time * bike2_travel_time
                    + Beta_bike2_cost * bike2_cost
                    + Beta_travel_purpose_income1 * (travel_purpose == 1) * (income == 1)
                    + Beta_weather_male1 * (weather == 1) * (gender == 1)
                    + lambda*LV)
  V[["V_bike3"]] = (ASC_bike3 + Beta_travel_purpose2 * (travel_purpose == 1)
                    + (travel_time == 1) * Beta_travel_time2
                    + Beta_weather2 * (weather == 1)
                    + Beta_bike3_walktime * bike3_walktime
                    + Beta_bike3_intime * bike3_intime
                    + Beta_bike3_travel_time * bike3_travel_time
                    + Beta_bike3_cost * bike3_cost
                    + Beta_travel_purpose_income2 * (travel_purpose == 1) * (income == 1)
                    + Beta_weather_male2 * (weather == 1) * (gender == 1)
                    + lambda*LV)
  
  V[["V_bus"]]= (ASC_bus + Beta_travel_purpose3 * (travel_purpose==1) + Beta_weather3 * (weather==1)
                 + (travel_time==1)  * Beta_travel_time3
                 + Beta_bus_waittime * bus_waittime
                 + Beta_bus_intime * bus_intime
                 + Beta_bus_traveltime * bus_traveltime
                 + Beta_bus_cost * bus_cost
                 + Beta_travel_purpose_income3 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male3 * (weather == 1) * (gender == 1) )
  
  V[["V_taxi"]]=(ASC_taxi + Beta_travel_purpose4 * (travel_purpose==1) + Beta_weather4 * (weather==1)
                 + (travel_time==1)  * Beta_travel_time4
                 + Beta_Taxi_waittime * Taxi_waittime
                 + Beta_Taxi_intime * Taxi_intime
                 + Beta_Taxi_traveltime * Taxi_traveltime
                 + Beta_Taxi_cost * Taxi_cost
                 + Beta_travel_purpose_income4 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male4 * (weather == 1) *(gender == 1)
  )
  
  V[["V_sharecar"]] =(ASC_sharecar + Beta_travel_purpose5 * (travel_purpose==1) + Beta_weather5 * (weather==1)
                      + (travel_time==1)  * Beta_travel_time5
                      + Beta_share_car_waittime * share_car_waittime
                      + Beta_share_car_intime * share_car_intime
                      + Beta_share_car_travel_time * share_car_travel_time
                      + Beta_share_car_cost * share_car_cost
                      + Beta_travel_purpose_income5 * (travel_purpose == 1) * (income==1)
                      + Beta_weather_male5 * (weather == 1) * (gender == 1)
  )
  
  V[["V_car"]] =(ASC_walk + Beta_travel_purpose7 * (travel_purpose==1)
                 + (travel_time==1) * Beta_travel_time7
                 + Beta_weather7 * (weather==1)
                 + Beta_walk_time * walk_time
                 + Beta_travel_purpose_income6 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male7 * (weather == 1) * (gender == 1))
  
  V[["V_walk"]] =(ASC_walk+  Beta_travel_purpose7* travel_purpose + Beta_weather7* weather +Beta_walk_time*walk_time 
                  +Beta_travel_purpose_income6*(travel_purpose==1)*( income==1)
                  +Beta_weather_male7* (weather==1) *  gender == 1)
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(
      V_bike1=1,
      V_bike2=2,
      V_bike3=3,
      V_bus=4,
      V_taxi=5,
      V_sharecar=6,
      V_car=7,
      V_walk=8), 
    avail         = 1,
    choiceVar     = choice,
    utilities     = V,
    componentName = "choice"
  )
  
  ### Compute probabilities for MNL model component
  P[["choice"]] = apollo_mnl(mnl_settings, functionality)
  
  ### Likelihood of the whole 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                                            ####
# ################################################################# #

### Optional: calculate LL before model estimation
# apollo_llCalc(apollo_beta, apollo_probabilities, apollo_inputs)

### Estimate model
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)               ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model)

# ################################################################# #
##### POST-PROCESSING                                            ####
# ################################################################# #

### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
apollo_sink()

# ----------------------------------------------------------------- #
#---- MODEL PREDICTIONS                                          ----
# ----------------------------------------------------------------- #


# ----------------------------------------------------------------- #
#---- CONDITIONALS AND UNCONDITIONALS                            ----
# ----------------------------------------------------------------- #

conditionals <- apollo_conditionals(model,apollo_probabilities,apollo_inputs)

summary(conditionals)

unconditionals <- apollo_unconditionals(model,apollo_probabilities,apollo_inputs)

mean(unconditionals[[1]])
sd(unconditionals[[1]])

# ----------------------------------------------------------------- #
#---- switch off writing to file                                 ----
# ----------------------------------------------------------------- #

apollo_sink()

[stephanehess]

Hi

could you please share the code and data offline with David and me and one of us will debug this for you

Stephane

[1123088145@qq.com]

Dear Stephane and David,
Thank you for your reply!
I would be glad to send you the code.
May I ask if it can be sent to the following email address?
s.hess@its.leeds.ac.uk

Thank you very much.
Ren.

[dpalma]

Hi Ren,

Please send it to both

• S.Hess [at] leeds.ac.uk

and

• D.Palma [at] leeds.ac.uk

Best wishes
David

[dpalma]

Hi Ren,

In R, when you are writing an expression across multiple lines, you need to put the connecting operator at the end of an incomplete line. For example, the following expression:

Code: Select all

x = 2 + 3 +
    5

Will cause x to contain the value 10.

Instead putting the connecting operator in the second line, as follows:

Code: Select all

x = 2 + 3
    + 5

will make R think these are two independent lines, and x will store the value 5.

In your code, you make that mistake when defining the latent variables inside apollo_randCoeff. Your definition of latent variables should look as below.

Code: Select all

apollo_randCoeff=function(apollo_beta, apollo_inputs){
  randcoeff = list()
  randcoeff[["LV"]] =   gamma_high_income * ( income == 0 ) + 
    gamma_highEducation * (edu == 0) + 
    gamma_occupation * (occupation == 3) + 
    gamma_age_35_less * (age == 1) + 
    gamma_male *  (gender == 1) + 
    gamma_more_car* (own_car>=2) + 
    gamma_no_car*  (own_car==0) + 
    gamma_more_bike* (own_bike >= 2) + 
    gamma_no_bike* (own_bike == 1) + 
    eta
  randcoeff[["LV1"]] =   gamma_high_income1 * ( income == 0 ) + 
    gamma_highEducation1 * (edu == 0) + 
    gamma_occupation1 * (occupation == 3) + 
    gamma_age_35_less1 * (age == 1) + 
    gamma_male1 *  (gender == 1) + 
    gamma_more_car1* (own_car>=2) + 
    gamma_no_car1*  (own_car==0) + 
    gamma_more_bike1* (own_bike >= 2) + 
    gamma_no_bike1* (own_bike == 1) + 
    eta1
  randcoeff[["LV2"]] =   gamma_high_income2 * ( income == 0 ) + 
    gamma_highEducation2 * (edu == 0) + 
    gamma_occupation2 * (occupation == 3) + 
    gamma_age_35_less2 * (age == 1) + 
    gamma_male2 *  (gender == 1) + 
    gamma_more_car2* (own_car>=2) + 
    gamma_no_car2*  (own_car==0) + 
    gamma_more_bike2* (own_bike >= 2) + 
    gamma_no_bike2* (own_bike == 1) + 
    eta2
  return(randcoeff)
}

However, even after fixing this, there are still issues with your model. Notably, there are several parameters that do not influence the likelihood. For example "gamma_no_car", "gamma_no_car1" and "gamma_no_car2" don't influence the utility because they are multiplied by (own_car==0), but own_car never takes a value equal to zero. Apollo will tell you which parameters have this issue, so you should remove them from your model.

Best wishes
David

[1123088145@qq.com]

Hi David
Thank you for your prompt attention.
Your suggested changes have been very helpful to me.
I have already sent the code and data to S.Hess@leeds.ac.uk and D.Palma@leeds.ac.uk 。
I hope you can help me debug this program.
We look forward to your reply.


Thank you very much.
Ren.

[dpalma]

Hi Ren,

I applied the changes I mentioned in the previous post to your code, and it is running now. Besides moving the "+" signs to the lines above in apollo_randCoeff, I had to remove (comment out) several parameters because they were not influencing the likelihood. See the code below.

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       = "Hybrid_with_OL",
  modelDescr      = "Hybrid choice model on drug choice data, using ordered measurement model for indicators",
  indivID         = "id",
  mixing          = TRUE,
  nCores          = 4, 
  outputDirectory = "output"
)

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

### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv), 
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("huizong.csv",header=TRUE)
#database = apollo_drugChoiceData
### for data dictionary, use ?apollo_drugChoiceData

# ################################################################# #
#### ANALYSIS OF CHOICES                                         ####
# ################################################################# #

### Illustration of how to use apollo_choiceAnalysis with user-defined alternatives.
### This is useful in cases where the alternatives in the data differ 
### across tasks. The same approach can then also be used with unlabelled data



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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
  lambda_at1=0,
  lambda_at2=0,
  lambda_at3=0,
  gamma_occupation =0,
  gamma_age_35_less =0,
  gamma_male = 0,
  gamma_highEducation =0,
  gamma_high_income = 0,
  gamma_more_car = 0,
  #gamma_no_car =  0,
  gamma_more_bike =0, 
  gamma_no_bike = 0,
  
  gamma_occupation1 =0,
  gamma_age_35_less1 =0,
  gamma_male1 = 0,
  gamma_highEducation1 =0,
  gamma_high_income1 = 0,
  gamma_more_car1 = 0,
  #gamma_no_car1 =  0,
  gamma_more_bike1 =0, 
  gamma_no_bike1 = 0,
  
  gamma_occupation2 =0,
  gamma_age_35_less2 =0,
  gamma_male2 = 0,
  gamma_highEducation2 =0,
  gamma_high_income2 = 0,
  gamma_more_car2 = 0,
  #gamma_no_car2 =  0,
  gamma_more_bike2 =0, 
  gamma_no_bike2 = 0,
  zeta_Availability1  = 1, 
  zeta_Tangibles1    = 1, 
  zeta_Benefit1       = 1, 
  zeta_Availability2  = 1, 
  zeta_Tangibles2     = 1, 
  zeta_Benefit2       = 1, 
  zeta_Availability3  = 1, 
  zeta_Tangibles3     = 1, 
  zeta_Benefit3       = 1, 
  zeta_Availability4  = 1, 
  zeta_Tangibles4     = 1, 
  tau_Availability_1  =-2, 
  tau_Availability_2   =-1, 
  tau_Availability_3   = 1, 
  tau_Availability_4    = 2, 
  
  tau_Tangibles_1  =-2, 
  tau_Tangibles_2  =-1, 
  tau_Tangibles_3  = 1, 
  tau_Tangibles_4  = 2, 
  tau_Benefit_1       =-2, 
  tau_Benefit_2       =-1, 
  tau_Benefit_3       = 1, 
  tau_Benefit_4       = 2, 
  
  ASC_bike1 = 0,
  ASC_bike2 = 0,
  ASC_bike3 = 0,
  ASC_bus = 0,
  ASC_taxi = 0,
  ASC_sharecar = 0,
  ASC_car = 0,
  ASC_walk =0 ,
  
  Beta_travel_time= 0,
  Beta_travel_time1= 0,
  Beta_travel_time2= 0,
  Beta_travel_time3= 0,
  Beta_travel_time4= 0,
  Beta_travel_time5= 0,
  Beta_travel_time6= 0,
  #Beta_travel_time7= 0,
  Beta_travel_purpose= 0,
  Beta_travel_purpose1= 0,
  Beta_travel_purpose2= 0,
  Beta_travel_purpose3=0,
  Beta_travel_purpose4=0,
  Beta_travel_purpose5= 0,
  Beta_travel_purpose6= 0,
  #Beta_travel_purpose7= 0,
  Beta_travel_purpose_income= 0,
  Beta_travel_purpose_income1= 0,
  Beta_travel_purpose_income2= 0,
  Beta_travel_purpose_income3=0,
  Beta_travel_purpose_income4=0,
  Beta_travel_purpose_income5= 0,
  Beta_travel_purpose_income6= 0,
  #Beta_travel_purpose_income7= 0,
  
  #Beta_travel_purpose_education= 0,
  #Beta_travel_purpose_education1= 0,
  #Beta_travel_purpose_education2= 0,
  #Beta_travel_purpose_education3=0,
  #Beta_travel_purpose_education4=0,
  #Beta_travel_purpose_education5= 0,
  #Beta_travel_purpose_education6= 0,
  #Beta_travel_purpose_education7= 0,
  Beta_weather= 0,
  Beta_weather1=0,
  Beta_weather2= 0,
  Beta_weather3= 0,
  Beta_weather4=0,
  Beta_weather5= 0,
  Beta_weather6= 0,
  #Beta_weather7= 0,
  Beta_weather_male= 0,
  Beta_weather_male1=0,
  Beta_weather_male2= 0,
  Beta_weather_male3= 0,
  Beta_weather_male4=0,
  Beta_weather_male5= 0,
  Beta_weather_male6= 0,
  #Beta_weather_male7= 0,
  Beta_bike_waittime=0,
  Beta_bike_intime= 0,
  Beta_bike1_travel_time=0,
  Beta_bike1_cost= 0,
  Beta_bike2_walkt_time=0,
  Beta_bike2_intime= 0,
  Beta_bike2_travel_time= 0,
  Beta_bike2_cost= 0,
  Beta_bike3_walktime= 0,
  Beta_bike3_intime= 0,
  Beta_bike3_travel_time=0,
  Beta_bike3_cost=0,
  Beta_bus_waittime= 0,
  Beta_bus_intime= 0,
  Beta_bus_traveltime= 0,
  Beta_bus_cost= 0,
  Beta_Taxi_waittime= 0,
  Beta_Taxi_intime= 0,
  Beta_Taxi_traveltime=0,
  Beta_Taxi_cost= 0,
  Beta_share_car_waittime= 0,
  Beta_share_car_intime= 0,
  Beta_share_car_travel_time=0,
  Beta_share_car_cost=0,
  Beta_carpostcar=0,
  Beta_car_traveltime= 0#,
  #Beta_walk_time=0
)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c("ASC_walk")

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

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType="mlhs", 
  interNDraws=100,          
  interUnifDraws=c(),      
  interNormDraws=c("eta","eta1","eta2") ,

  intraDrawsType="",
  intraNDraws=0,          
  intraUnifDraws=c(),     
  intraNormDraws=c() 
)

### Create random parameters  
apollo_randCoeff=function(apollo_beta, apollo_inputs){
  randcoeff = list()
  randcoeff[["LV"]] =   gamma_high_income * ( income == 0 ) + 
    gamma_highEducation * (edu == 0) + 
    gamma_occupation * (occupation == 3) + 
    gamma_age_35_less * (age == 1) + 
    gamma_male *  (gender == 1) + 
    gamma_more_car* (own_car>=2) + 
    #gamma_no_car*  (own_car==0) + 
    gamma_more_bike* (own_bike >= 2) + 
    gamma_no_bike* (own_bike == 1) + 
    eta
  randcoeff[["LV1"]] =   gamma_high_income1 * ( income == 0 ) + 
    gamma_highEducation1 * (edu == 0) + 
    gamma_occupation1 * (occupation == 3) + 
    gamma_age_35_less1 * (age == 1) + 
    gamma_male1 *  (gender == 1) + 
    gamma_more_car1* (own_car>=2) + 
    #gamma_no_car1*  (own_car==0) + 
    gamma_more_bike1* (own_bike >= 2) + 
    gamma_no_bike1* (own_bike == 1) + 
    eta1
  randcoeff[["LV2"]] =   gamma_high_income2 * ( income == 0 ) + 
    gamma_highEducation2 * (edu == 0) + 
    gamma_occupation2 * (occupation == 3) + 
    gamma_age_35_less2 * (age == 1) + 
    gamma_male2 *  (gender == 1) + 
    gamma_more_car2* (own_car>=2) + 
    #gamma_no_car2*  (own_car==0) + 
    gamma_more_bike2* (own_bike >= 2) + 
    gamma_no_bike2* (own_bike == 1) + 
    eta2
  return(randcoeff)
}

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

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION      测量方程                  ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P = list()
  
  ### Likelihood of indicators 
  ## 
  ol_settings1 = list(outcomeOrdered = A3, 
                      V              = zeta_Availability1*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4))
                      #componentName  = "indic_Availability1")
  ol_settings2 = list(outcomeOrdered = A2, 
                      V              = zeta_Availability2*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4))
                      #componentName  = "indic_Availability2")
  ol_settings3 = list(outcomeOrdered = A21, 
                      V              = zeta_Availability3*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4))
                      
                      #componentName  = "indic_Availability3")
  ol_settings4 = list(outcomeOrdered = A16, 
                      V              = zeta_Availability4*LV, 
                      tau            = list(tau_Availability_1, tau_Availability_2, tau_Availability_3, tau_Availability_4))
                      #componentName  = "indic_Availability4")
  ol_settings5 = list(outcomeOrdered = A10, 
                      V              = zeta_Tangibles1*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4))
                      
                      #componentName  = "indic_Tangibles1")
  ol_settings6 = list(outcomeOrdered = A11, 
                      V              = zeta_Tangibles2*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4))
                      #componentName  = "indic_Tangibles2")
  ol_settings7 = list(outcomeOrdered = A4, 
                      V              = zeta_Tangibles3*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4))
                      #componentName  = "indic_Tangibles3")
  ol_settings8 = list(outcomeOrdered = A17, 
                      V              = zeta_Tangibles4*LV1, 
                      tau            = list(tau_Tangibles_1, tau_Tangibles_2, tau_Tangibles_3, tau_Tangibles_4))
                      #componentName  = "indic_Tangibles4")
  ol_settings9 = list(outcomeOrdered = A22, 
                      V              = zeta_Benefit1*LV2, 
                      tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4))
                      #componentName  = "indic_Benefit1")
  ol_settings10 = list(outcomeOrdered = A24, 
                       V              = zeta_Benefit2*LV2, 
                       tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4))
                       #componentName  = "indic_Benefit2")
  ol_settings11 = list(outcomeOrdered = A18, 
                       V              = zeta_Benefit3*LV2, 
                       tau            = list(tau_Benefit_1, tau_Benefit_2, tau_Benefit_3, tau_Benefit_4))  
                       #componentName  = "indic_Benefit3")
  P[["indic_Availability1"]]= apollo_ol(ol_settings1, functionality)
  P[["indic_Availability2"]]= apollo_ol(ol_settings2, functionality)
  P[["indic_Availability3"]]= apollo_ol(ol_settings3, functionality)
  P[["indic_Availability4"]] = apollo_ol(ol_settings4, functionality)
  P[["indic_Tangibles1"]] = apollo_ol(ol_settings5, functionality)
  P[["indic_Tangibles2"]] = apollo_ol(ol_settings6, functionality)
  P[["indic_Tangibles3"]] = apollo_ol(ol_settings7, functionality)
  P[["indic_Tangibles4"]] = apollo_ol(ol_settings8, functionality)
  P[["indic_Benefit1"]]   = apollo_ol(ol_settings9, functionality)
  P[["indic_Benefit2"]]   = apollo_ol(ol_settings10, functionality)
  P[["indic_Benefit3"]]   = apollo_ol(ol_settings11, functionality)
  
  ### Likelihood of choices
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["V_bike1"]]= (ASC_bike1 + Beta_travel_purpose * (travel_purpose == 1)
                   + Beta_weather * (weather == 1)
                   + (travel_time == 1) * Beta_travel_time
                   + Beta_bike_waittime * bike_waittime
                   + Beta_bike_intime * bike_intime
                   + Beta_bike1_travel_time * bike1_travel_time
                   + Beta_bike1_cost * bike1_cost
                   + Beta_travel_purpose_income * (travel_purpose == 1) * (income == 1)
                   + Beta_weather_male * (weather == 1) * (gender == 1)
                   + lambda_at1*LV
                   +lambda_at2*LV1 
                   + lambda_at3 * LV2
  )
  
  
  V[["V_bike2"]] = (ASC_bike2 + Beta_travel_purpose1 * (travel_purpose == 1)
                    + Beta_weather1 * (weather == 1)
                    + (travel_time == 1) * Beta_travel_time1
                    + Beta_bike2_walkt_time * bike2_walkt_time
                    + Beta_bike2_intime * bike2_intime
                    + Beta_bike2_travel_time * bike2_travel_time
                    + Beta_bike2_cost * bike2_cost
                    + Beta_travel_purpose_income1 * (travel_purpose == 1) * (income == 1)
                    + Beta_weather_male1 * (weather == 1) * (gender == 1)
                    + lambda_at1*LV
                    +lambda_at2*LV1 
                    + lambda_at3 * LV2)
  V[["V_bike3"]] = (ASC_bike3 + Beta_travel_purpose2 * (travel_purpose == 1)
                    + (travel_time == 1) * Beta_travel_time2
                    + Beta_weather2 * (weather == 1)
                    + Beta_bike3_walktime * bike3_walktime
                    + Beta_bike3_intime * bike3_intime
                    + Beta_bike3_travel_time * bike3_travel_time
                    + Beta_bike3_cost * bike3_cost
                    + Beta_travel_purpose_income2 * (travel_purpose == 1) * (income == 1)
                    + Beta_weather_male2 * (weather == 1) * (gender == 1)
                    + lambda_at1*LV
                    + lambda_at2*LV1 
                    + lambda_at3 * LV2)
  
  V[["V_bus"]]= (ASC_bus + Beta_travel_purpose3 * (travel_purpose==1) + Beta_weather3 * (weather==1)
                 + (travel_time==1)  * Beta_travel_time3
                 + Beta_bus_waittime * bus_waittime
                 + Beta_bus_intime * bus_intime
                 + Beta_bus_traveltime * bus_traveltime
                 + Beta_bus_cost * bus_cost
                 + Beta_travel_purpose_income3 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male3 * (weather == 1) * (gender == 1) )
  
  V[["V_taxi"]]=(ASC_taxi + Beta_travel_purpose4 * (travel_purpose==1) + Beta_weather4 * (weather==1)
                 + (travel_time==1)  * Beta_travel_time4
                 + Beta_Taxi_waittime * Taxi_waittime
                 + Beta_Taxi_intime * Taxi_intime
                 + Beta_Taxi_traveltime * Taxi_traveltime
                 + Beta_Taxi_cost * Taxi_cost
                 + Beta_travel_purpose_income4 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male4 * (weather == 1) *(gender == 1)
  )
  
  V[["V_sharecar"]] =(ASC_sharecar + Beta_travel_purpose5 * (travel_purpose==1) + Beta_weather5 * (weather==1)
                      + (travel_time==1)  * Beta_travel_time5
                      + Beta_share_car_waittime * share_car_waittime
                      + Beta_share_car_intime * share_car_intime
                      + Beta_share_car_travel_time * share_car_travel_time
                      + Beta_share_car_cost * share_car_cost
                      + Beta_travel_purpose_income5 * (travel_purpose == 1) * (income==1)
                      + Beta_weather_male5 * (weather == 1) * (gender == 1)
  )
  
  V[["V_car"]] =(ASC_car + Beta_travel_purpose6 * (travel_purpose==1)
                 + Beta_weather6 * (weather==1)
                 + (travel_time==1)  * Beta_travel_time6
                 + Beta_carpostcar * carpostcar + Beta_car_traveltime * car_traveltime
                 + Beta_travel_purpose_income6 * (travel_purpose == 1) * (income==1)
                 + Beta_weather_male6 * (weather == 1) * (gender == 1)
  )
  
  V[["V_walk"]] =ASC_walk
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(
      V_bike1=1,
      V_bike2=2,
      V_bike3=3,
      V_bus=4,
      V_taxi=5,
      V_sharecar=6,
      V_car=7,
      V_walk=8), 
    avail         = 1,
    choiceVar     = choice,
    utilities     = V,
    componentName = "choice"
  )
  
  ### Compute probabilities for MNL model component
  P[["choice"]] = apollo_mnl(mnl_settings, functionality)
  
  ### Likelihood of the whole 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                                            ####
# ################################################################# #

### Optional: calculate LL before model estimation
# apollo_llCalc(apollo_beta, apollo_probabilities, apollo_inputs)

### Estimate model
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)               ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model)

# ################################################################# #
##### POST-PROCESSING                                            ####
# ################################################################# #

### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
apollo_sink()

# ----------------------------------------------------------------- #
#---- MODEL PREDICTIONS                                          ----
# ----------------------------------------------------------------- #


# ----------------------------------------------------------------- #
#---- CONDITIONALS AND UNCONDITIONALS                            ----
# ----------------------------------------------------------------- #

conditionals <- apollo_conditionals(model,apollo_probabilities,apollo_inputs)

summary(conditionals)

unconditionals <- apollo_unconditionals(model,apollo_probabilities,apollo_inputs)

mean(unconditionals[[1]])
sd(unconditionals[[1]])

# ----------------------------------------------------------------- #
#---- switch off writing to file                                 ----
# ----------------------------------------------------------------- #

apollo_sink()

Best wishes
David