ApolloChoiceModelling forum

Code: Select all

WARNING: Object apollo_randCoeff from the global environment should not be used inside apollo_probabilities. If needed, save it inside apollo_inputs and call it as
  apollo_inputs$apollo_randCoeff. 
Preparing user-defined functions.
[1] 6.512178
[1] 6.512178
[1] 6.512078
[1] 6.512078
[1] 6.512078
[1] 6.512178
[1] 6.512178
[1] 6.512178
Error in value[[3L]](cond): argument is not used (cond) 错误于value[[3L]](cond): 参数没有用(cond)

Code: Select all

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

### Load Apollo library
library(apollo)

###Preparing environment
apollo_initialise()

###Options controlling the running of the code
apollo_control = list(
    modelName  ="30_0421Class_RUM_PRRM_Inertia_witnLV",
  modelDescr ="Latent class with 2 classes: RUM_PRRMwithInertiawithLV",
  indivID    ="ID",
  nCores     = 6 ,
  HB = FALSE,
  mixing = TRUE,
  outputDirectory = "C:/Users/zuoxi/Desktop/LC_ModelEstimation/LC_errorcopmonent/output"
)

###Reading in database
database = read.csv("C:/Users/zuoxi/Desktop/LC_modelestimation/LC_RUM_PRRM_2classes/data/250422_historical01_diff.csv", header = TRUE)

###Parameters to be estimated and their starting values
apollo_beta = c(
 # ----- Structural equations -----
 # socio-economic 
  coef_intercept_metrohabit = 5,
  coef_Gender_metrohabit = -1,
  coef_Age18_24_metrohabit = -0.7,
  coef_Age25_34_metrohabit = -1,
  coef_Age35_44_metrohabit = -1,
  coef_PrivateCarOwned0_metrohabit = -2,
  coef_PrivateCarOwned2_metrohabit  = 2,
  coef_NonmotorOwned_metrohabit = 0,
  coef_Income3k_metrohabit = 3.94,
  coef_Income3_6k_metrohabit = -2.02,
  coef_Income6_10k_metrohabit = -3.49,
  coef_CareerStudent_metrohabit = 5.21,
  coef_CareerEmployed_metrohabit = -0.81,
  coef_CareerFreelancer_metrohabit = -5.15,
  coef_TransitCard_metrohabit = 0.4,
  coef_Bike_sharingExperience_metrohabit = 0.004,
  sigma_metrohabit = 0,

 ### Influence of LV on class 1
  Alt1_CL_1 = 0.1,
  Alt2_CL_1 = 0.01,
  Alt3_CL_1 = 0.001,
   ### Influence of LV on class 1
  Alt1_CL_2 = 0.1,
  Alt2_CL_2 = 0.01,
  Alt3_CL_2 = 0.001,

 # ----- Measurement equations ----- 
  ## Habit indicators (binary)
  #MaxFrequenceModeMetro, MetroFrequence1_2, MetroFrequence7_10, TravelPurposeMetroWork, TravelPurposeMetroLeisure
  int_MaxFrequenceModeMetro = 0.002,
  int_MetroFrequence1_2 = 0.002,# fixed
  int_MetroFrequence7_10 = 0.003,
  int_TravelPurposeMetroWork = 0.001,
  int_TravelPurposeMetroLeisure = 0.002, 
  zeta_MaxFrequenceModeMetro = 0.002,   
  zeta_MetroFrequence1_2  = -1,# fixed
  zeta_MetroFrequence7_10 = 0.001,
  zeta_TravelPurposeMetroWork = 0.003,
  zeta_TravelPurposeMetroLeisure = -1,
  

 # -----intraClass a Alternative-spercific-------
  TIME_1 = -1.948067,
  COST_1 = -0.311116,
  TRANSFERDISTANCE_1  = -0.235411,
  TRANSFERWAITTIME_1  = -0.857142,
  
  ## Class a URM ASC
  ASC_A = 1.39974,
  ASC_B = 1.979026,
  ASC_C = -1.43927,
  ASC_D = 2.49566,
  ASC_E = 2.45503,
  ASC_F = 0,    #Reference
 
 # -----intraClass b Alternative probability--------
  TIME_2 = 0.760591,
  COST_2 = 0.98366,
  TRANSFERDISTANCE_2  = 0.735376,
  TRANSFERWAITTIME_2  = 0.966079,
  
  #Inertia effect
  BETA_Inertia_T_2  = -0.001,
  BETA_Inertia_C_2  = -0.001,
  #Inertiagamma_2  = 0.317839,

 # ----- Class membership parameters-----
  ## Class a
  delta_a = 0,
  GENDER_a = 0,
  AGE18_24_a = 0,
  AGE25_34_a = 0,
  AGE35_44_a = 0,
  PRICAROWN0_a = 0,
  PRICAROWN2_a = 0,
  NONMOTEROWN_a = 0,
  INCOME0_3k_a = 0,
  INCOME3_6k_a = 0,
  INCOME6_10k_a = 0,
  CAREER_STUDENT_a = 0,
  CAREER_EMPLOYED_a = 0,
  CAREER_FREELANCER_a = 0,
  TRANSITCARD_a = 0,
  BIKESHARINGEXP_a = 0,
  TAXIAPPEXP_a = 0,
  DISRUPTIONEXP_a = 0,
  ## Class b
  delta_b     = -1.09,
  GENDER_b = 0.2,
  AGE18_24_b = -0.430993,
  AGE25_34_b = -0.730195,
  AGE35_44_b = 0.3,
  PRICAROWN0_b = 0.070752,
  PRICAROWN2_b = 0.160588,
  NONMOTEROWN_b = -0.4,
  INCOME0_3k_b = 0.16,
  INCOME3_6k_b = -0.3,
  INCOME6_10k_b = -0.161534,
  CAREER_STUDENT_b = -0.300524,
  CAREER_EMPLOYED_b = 0.382501,
  CAREER_FREELANCER_b = -0.010985,
  TRANSITCARD_b = -0.044079,
  BIKESHARINGEXP_b = 0.208017,
  TAXIAPPEXP_b = -0.603442,
  DISRUPTIONEXP_b = -0.228155
)
###Fixed parameters: should be in quotes 
apollo_fixed = c("ASC_F",
                 "delta_a","GENDER_a","AGE18_24_a","AGE25_34_a",
                 "AGE35_44_a","PRICAROWN0_a","PRICAROWN2_a","NONMOTEROWN_a",
                 "INCOME0_3k_a","INCOME3_6k_a","INCOME6_10k_a",
                 "CAREER_STUDENT_a","CAREER_EMPLOYED_a","CAREER_FREELANCER_a",
                 "TRANSITCARD_a", "BIKESHARINGEXP_a","TAXIAPPEXP_a","DISRUPTIONEXP_a", 
                 "int_MetroFrequence1_2", "zeta_MetroFrequence1_2") 

### DEFINE RANDOM COMPONENTS ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "sobolOwen",
  interNDraws    = 600,
  interUnifDraws = c(),
  interNormDraws = c("draw_metrohabit"),
  intraDrawsType = "sobolOwen",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
 randcoeff = list()
 randcoeff[["LV_metrohabit"]] =     
      coef_intercept_metrohabit
    + coef_Gender_metrohabit *  Gender
    + coef_Age18_24_metrohabit * Age18_24
    + coef_Age25_34_metrohabit * Age25_34
    + coef_Age35_44_metrohabit * Age35_44
    + coef_PrivateCarOwned0_metrohabit * PrivateCarOwned0
    + coef_PrivateCarOwned2_metrohabit * PrivateCarOwned2
    + coef_NonmotorOwned_metrohabit * NonmotorOwned
    + coef_Income3k_metrohabit * Income3k
    + coef_Income3_6k_metrohabit * Income3_6k
    + coef_Income6_10k_metrohabit * Income6_10k
    + coef_CareerStudent_metrohabit * CareerStudent
    + coef_CareerEmployed_metrohabit * CareerEmployed
    + coef_CareerFreelancer_metrohabit * CareerFreelancer
    + coef_TransitCard_metrohabit * TransitCard
    + coef_Bike_sharingExperience_metrohabit * Bike_sharingExperience
    + sigma_metrohabit * draw_metrohabit
return(randcoeff)
}


###Grouping latent class parameters 不同组的相同参数在一个list
apollo_lcPars = function(apollo_beta, apollo_inputs){
  lcpars = list()
  lcpars[["TIME"]] = list(TIME_1, TIME_2)
  lcpars[["COST"]] = list(COST_1, COST_2)
  lcpars[["TRANSFERDISTANCE"]]  = list(TRANSFERDISTANCE_1, TRANSFERDISTANCE_2)
  lcpars[["TRANSFERWAITTIME"]]  = list(TRANSFERWAITTIME_1, TRANSFERWAITTIME_2)
  lcpars[["Alt1_CL"]]  = list(Alt1_CL_1, Alt1_CL_2)
  lcpars[["Alt2_CL"]]  = list(Alt2_CL_1, Alt2_CL_2)
  lcpars[["Alt3_CL"]]  = list(Alt3_CL_1, Alt3_CL_2)
  
  ###Class membership probabilities based on s_1, s_2: use of MNL fomula
  V=list()
  
  V[["class_a"]] = delta_a + GENDER_a * Gender + AGE18_24_a * Age18_24 +
    AGE25_34_a * Age25_34 + AGE35_44_a * Age35_44 + PRICAROWN0_a * PrivateCarOwned0 +
    PRICAROWN2_a * PrivateCarOwned2 + NONMOTEROWN_a * NonmotorOwned +
    INCOME0_3k_a * Income3k + INCOME3_6k_a * Income3_6k + INCOME6_10k_a * Income6_10k +
    CAREER_STUDENT_a * CareerStudent + CAREER_EMPLOYED_a * CareerEmployed +
    CAREER_FREELANCER_a * CareerFreelancer + TRANSITCARD_a * TransitCard+
    BIKESHARINGEXP_a * Bike_sharingExperience + TAXIAPPEXP_a * TaxiAppExperience + DISRUPTIONEXP_a * DisruptionExperience

  V[["class_b"]] = delta_b + GENDER_b * Gender + AGE18_24_b * Age18_24 +
    AGE25_34_b * Age25_34 + AGE35_44_b * Age35_44 + PRICAROWN0_b * PrivateCarOwned0 +
    PRICAROWN2_b * PrivateCarOwned2 + NONMOTEROWN_b * NonmotorOwned +
    INCOME0_3k_b * Income3k + INCOME3_6k_b * Income3_6k + INCOME6_10k_b * Income6_10k +
    CAREER_STUDENT_b * CareerStudent + CAREER_EMPLOYED_b * CareerEmployed +
     CAREER_FREELANCER_b * CareerFreelancer + TRANSITCARD_b * TransitCard +
    BIKESHARINGEXP_b * Bike_sharingExperience + TAXIAPPEXP_b * TaxiAppExperience + DISRUPTIONEXP_b * DisruptionExperience
  
  
  ###Settings for class membership probabilities 
  classAlloc_settings = list(
    classes      = c(class_a=1, class_b=2), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  ###Class membership probabilities
  lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
  
  return(lcpars)
}

###Search the user work space for all necessary input 
apollo_inputs = apollo_validateInputs()
apollo_inputs$apollo_randCoeff = apollo_randCoeff
###Probability function
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ###Attaches parameters and data so that variables can be referred by name
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  randcoeff = apollo_randCoeff(apollo_beta, apollo_inputs)
  print(head(randcoeff$LV_metrohabit))  # View the first few rows of the latent variable values
  P <- list()

# ----------------------------------------------------------------- #
#---- Likelihood of binary Probit indicators for LV_metrohabit ----#
# ----------------------------------------------------------------- #

mnl_settings1 = list(
  alternatives = c(chosen = 1, not_chosen = 0), # Must define two options
  avail = list(chosen = 1, not_chosen = 1),     # Both options are available
  choiceVar = MaxFrequenceModeMetro,            # Key correction: specify observed variable
  V = list(
    chosen = int_MaxFrequenceModeMetro + zeta_MaxFrequenceModeMetro * LV_metrohabit,
    not_chosen = 0  # Reference category is set to 0
  ),
  componentName = "indic_MaxFrequenceModeMetro"
)

mnl_settings2 = list(
  alternatives = c(chosen = 1, not_chosen = 0),
  avail = list(chosen = 1, not_chosen = 1),
  choiceVar = MetroFrequence1_2,
  V = list(
    chosen = int_MetroFrequence1_2 + zeta_MetroFrequence1_2 * LV_metrohabit,
    not_chosen = 0
  ),
  componentName = "indic_MetroFrequence1_2"
)

# Indicator 3: Metro usage frequency is 7–10 times/week
mnl_settings3 = list(
  alternatives = c(chosen = 1, not_chosen = 0),
  avail = list(chosen = 1, not_chosen = 1),
  choiceVar = MetroFrequence7_10,
  V = list(
    chosen = int_MetroFrequence7_10 + zeta_MetroFrequence7_10 * LV_metrohabit,
    not_chosen = 0
  ),
  componentName = "indic_MetroFrequence7_10"
)

# Indicator 4: Trip purpose is commuting
mnl_settings4 = list(
  alternatives = c(chosen = 1, not_chosen = 0),
  avail = list(chosen = 1, not_chosen = 1),
  choiceVar = TravelPurposeMetroWork,
  V = list(
    chosen = int_TravelPurposeMetroWork + zeta_TravelPurposeMetroWork * LV_metrohabit,
    not_chosen = 0
  ),
  componentName = "indic_TravelPurposeMetroWork"
)

# Indicator 5: Trip purpose is leisure
mnl_settings5 = list(
  alternatives = c(chosen = 1, not_chosen = 0),
  avail = list(chosen = 1, not_chosen = 1),
  choiceVar = TravelPurposeMetroLeisure,
  V = list(
    chosen = int_TravelPurposeMetroLeisure + zeta_TravelPurposeMetroLeisure * LV_metrohabit,
    not_chosen = 0
  ),
  componentName = "indic_TravelPurposeMetroLeisure"
)

# Add the five indicators to the probability function
# (functionality is a variable, usually passed as an argument)
P[["indic_MaxFrequenceModeMetro"]]     = apollo_mnl(mnl_settings1, functionality)
P[["indic_MetroFrequence1_2"]]         = apollo_mnl(mnl_settings2, functionality)
P[["indic_MetroFrequence7_10"]]        = apollo_mnl(mnl_settings3, functionality)
P[["indic_TravelPurposeMetroWork"]]    = apollo_mnl(mnl_settings4, functionality)
P[["indic_TravelPurposeMetroLeisure"]] = apollo_mnl(mnl_settings5, functionality)

# Combine the models
P = apollo_combineModels(P, apollo_inputs, functionality)

# Take the product of probabilities across observations for the same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
   ### Rename model
  names(P)[which(names(P)=="model")] <- "Measurement_model"
  # ----------------------------------------------------------------- #
  #---- LC model
  # ----------------------------------------------------------------- #
  ###Input for calculating MNL probabilities
  mnl_settings = list(
    alternatives  = c(Alt1=1, Alt2=2, Alt3=3, Alt4=4, Alt5=5, Alt6=6), 
    avail         = list(Alt1=1, Alt2=1, Alt3=1, Alt4=1, Alt5=1, Alt6=1),
    choiceVar     = CHOICE
  )

  ### Compute class-specific utilities
  V=list()
  V[['Alt1']]  = ASC_A + Alt1_CL_1 * LV_metrohabit + TIME_1 * time1_sc + COST_1 * cost1_sc + TRANSFERDISTANCE_1 * transferdistance1_sc + TRANSFERWAITTIME_1 * transferwaittime1_sc 
  V[['Alt2']]  = ASC_B + Alt2_CL_1 * LV_metrohabit + TIME_1 * time2_sc + COST_1 * cost2_sc + TRANSFERDISTANCE_1 * transferdistance2_sc + TRANSFERWAITTIME_1 * transferwaittime2_sc 
  V[['Alt3']]  = ASC_C + Alt3_CL_1 * LV_metrohabit + TIME_1 * time3_sc + COST_1 * cost3_sc + TRANSFERDISTANCE_1 * transferdistance3_sc + TRANSFERWAITTIME_1 * transferwaittime3_sc 
  V[['Alt4']]  = ASC_D + TIME_1 * time4_sc + COST_1 * cost4_sc + TRANSFERDISTANCE_1 * transferdistance4_sc + TRANSFERWAITTIME_1 * transferwaittime4_sc
  V[['Alt5']]  = ASC_E + TIME_1 * time5_sc + COST_1 * cost5_sc + TRANSFERDISTANCE_1 * transferdistance5_sc + TRANSFERWAITTIME_1 * transferwaittime5_sc
  V[['Alt6']]  = ASC_F
  
  ###Calculating probabilities based on MNL function for class 1
  mnl_settings$V = V
  P[[1]] = apollo_mnl(mnl_settings, functionality)
  P[[1]] = apollo_panelProd(P[[1]], apollo_inputs ,functionality)
  
  ### Compute class-specific regrets
  R=list()
  R[['Alt1']]  = TIME_2 * X_time1 + COST_2 * X_cost1 + TRANSFERDISTANCE_2 *  X_transferdistance1 + TRANSFERWAITTIME_2 * X_transferwaittime1 +
    X_Time_diff1 * BETA_Inertia_T_2 + X_Cost_diff1 * BETA_Inertia_C_2 + Alt1_CL_2 * LV_metrohabit
  R[['Alt2']]  = TIME_2 * X_time2 + COST_2 * X_cost2 + TRANSFERDISTANCE_2 *  X_transferdistance2 + TRANSFERWAITTIME_2 * X_transferwaittime2 +
    X_Time_diff2 * BETA_Inertia_T_2 + X_Cost_diff2 * BETA_Inertia_C_2 + Alt2_CL_2 * LV_metrohabit
  R[['Alt3']]  = TIME_2 * X_time3 + COST_2 * X_cost3 + TRANSFERDISTANCE_2 *  X_transferdistance3 + TRANSFERWAITTIME_2 * X_transferwaittime3 +
    X_Time_diff3 * BETA_Inertia_T_2 + X_Cost_diff3 * BETA_Inertia_C_2 + Alt3_CL_2 * LV_metrohabit
  R[['Alt4']]  = TIME_2 * X_time4 + COST_2 * X_cost4 + TRANSFERDISTANCE_2 *  X_transferdistance4 + TRANSFERWAITTIME_2 * X_transferwaittime4 +
    X_Time_diff4 * BETA_Inertia_T_2 + X_Cost_diff4 * BETA_Inertia_C_2
  R[['Alt5']]  = TIME_2 * X_time5 + COST_2 * X_cost5 + TRANSFERDISTANCE_2 *  X_transferdistance5 + TRANSFERWAITTIME_2 * X_transferwaittime5 +
    X_Time_diff5 * BETA_Inertia_T_2 + X_Cost_diff5 * BETA_Inertia_C_2
  R[['Alt6']]  = TIME_2 * X_time6 + COST_2 * X_cost6 + TRANSFERDISTANCE_2 *  X_transferdistance6 + TRANSFERWAITTIME_2 * X_transferwaittime6 +
    X_Time_diff6 * BETA_Inertia_T_2 + X_Cost_diff6 * BETA_Inertia_C_2 
  
  ###Calculating probabilities based on MNL function for class 2
  mnl_settings$V = lapply(R, "*", -1) ###the regrets must be negative (and used in MNL-settings as V)
  mnl_settings$V = R
  P[[2]] = apollo_mnl(mnl_settings, functionality)
  P[[2]] = apollo_panelProd(P[[2]], apollo_inputs ,functionality)

  ###Calculating choice probabilities using class membership and conditional probabilities 
  lc_settings   = list(inClassProb = P, classProb=pi_values)
  P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
  
  ##Average across inter-individual draws
  #P = apollo_avgInterDraws(P, apollo_inputs, functionality)
  
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}


### Optional: searching for starting value
apollo_beta = apollo_searchStart(apollo_beta,
                                 apollo_fixed,
                                 apollo_probabilities,
                                 apollo_inputs,
                                 searchStart_settings=list(nCandidates=20))
###Estimate the model  
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings=list(hessianRoutine="maxLik"))

###Display the output to console with p-values
apollo_modelOutput(model)

### Save output to file(s)
apollo_saveOutput(model)