I’m currently estimating a latent class model with one latent variable in Apollo, and I’ve encountered an issue during the estimation phase. I would really appreciate any insights or suggestions you might have.
Here’s the setup of my model:
1.It has two latent classes:
Class 1 follows the RUM (Random Utility Maximization) principle
Class 2 follows the RRM (Random Regret Minimization) principle
2.There is one latent variable: “Metro habit,” which is measured using five indicators representing the traveler’s historical use of the metro, all the indicators are binary variables (0 or 1).The latent variable is included in the utility functions of both classes.
However, when I run the estimation, I receive the following warning and error messages:
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)
I am a bit confused about this warning and error, especially regarding how to correctly use "apollo_randCoeff" within "apollo_probabilities". I tried to revise the code as suggested but the problem persists.
The full code for your reference.
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)
Thanks in advance for your help!
Best regards.