Errors when estimating a LC model with random components allowed for one of the classes
Posted: 19 Jan 2021, 20:07
Dear Stephane and David,
Thank you very much for inventing the Apollo package and creating the forum allowing our users to communicate.
I recently run a latent class model, but allow for inter-individual heterogeneity only within one of the classes (i.e., preferences are assumed to be non-random in the other classes). The model ends up with an error saying
"Error in apollo_avgInterDraws(P[[s]], apollo_inputs, functionality) :
No Inter-individuals draws to average over! "
I think the problem is the following: R tries to average across inter-individual draws for class, but the problem occurs when estimating the two non-random classes where no draws are generated.
My questions are:
(1) Is this model theoretically valid and/or practically sensible?
(2) If yes to (1), How can I code the model such as this (where let's say, attributes follow are random in one class, and non-random in the other two classes)?
I have attached my R code below.
Thank you very much.
Tim
-------------------------------------------
# ################################################################# #
#### 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 ="RPL_coop_LC3_ENVgroupRPL",
modelDescr ="RPL_coop_LC3_ENVgroupRPL",
indivID ="pool_ID",
nCores = 3,
mixing = TRUE,
noDiagnostics = TRUE
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
#load the data
database = read.csv("R_for_coop.csv",header=TRUE)
##renames##
#names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$dce_valid_CB1==1&dce_valid_CB2==1&
dce_valid_CB3==1&dce_valid_CB4==1&dce_valid_T2B1==1
&dce_valid_T2B2==1&dce_valid_T2B3==1&dce_valid_T2B4==1)
############################T1############################################
####T1 and control only
database= subset(database, database$T2==0)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_I = 0,
asc_III_1 = 0,
asc_III_2_mu = 0,
asc_III_2_sig = 0,
asc_III_3 = 0,
project_1 = 0,
project_2_mu = 0,
project_2_sig = 0,
project_3 = 0,
carbon_1 = 0,
carbon_2_mu = 0,
carbon_2_sig = 0,
carbon_3 = 0,
location_Wregion_1 = 0,
location_Wregion_2_mu = 0,
location_Wregion_2_sig = 0,
location_Wregion_3 = 0,
location_Wcountry_1 = 0,
location_Wcountry_2_mu = 0,
location_Wcountry_2_sig = 0,
location_Wcountry_3 = 0,
location_Ocountry_1 = 0,
location_Ocountry_2_mu = 0,
location_Ocountry_2_sig = 0,
location_Ocountry_3 = 0,
min_invest_1 = 0,
min_invest_2_mu = 0,
min_invest_2_sig = 0,
min_invest_3 = 0,
duration_I_1 = 0,
duration_I_2_mu = 0,
duration_I_2_sig = 0,
duration_I_3 = 0,
duration_II_1 = 0,
duration_II_2_mu = 0,
duration_II_2_sig = 0,
duration_II_3 = 0,
duration_V_1 = 0,
duration_V_2_mu = 0,
duration_V_2_sig = 0,
duration_V_3 = 0,
Q_participation_1 = 0,
Q_participation_2_mu = 0,
Q_participation_2_sig = 0,
Q_participation_3 = 0,
A_participation_1 = 0,
A_participation_2_mu = 0,
A_participation_2_sig = 0,
A_participation_3 = 0,
e_return_1 = 0,
e_return_2_mu = 0,
e_return_2_sig = 0,
e_return_3 = 0,
delta_a = 0,
delta_b = 0,
delta_c = 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("delta_b", "asc_I") #fix the scale of the control treatment
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("draws_project_2",
"draws_carbon_2",
"draws_location_Wregion_2",
"draws_location_Wcountry_2",
"draws_location_Ocountry_2",
"draws_min_invest_2",
"draws_duration_I_2",
"draws_duration_II_2",
"draws_duration_V_2",
"draws_Q_participation_2",
"draws_A_participation_2",
"draws_asc_2",
"draws_Ereturn_2"
),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["asc_III_2"]] = asc_III_2_mu + asc_III_2_sig*draws_asc_2
randcoeff[["project_2"]] = project_2_mu + project_2_sig*draws_project_2
randcoeff[["carbon_2"]] = carbon_2_mu + carbon_2_sig*draws_carbon_2
randcoeff[["location_Wregion_2"]] = location_Wregion_2_mu + location_Wregion_2_sig*draws_location_Wregion_2
randcoeff[["location_Wcountry_2"]] = location_Wcountry_2_mu + location_Wcountry_2_sig*draws_location_Wcountry_2
randcoeff[["location_Ocountry_2"]] = location_Ocountry_2_mu + location_Ocountry_2_sig*draws_location_Ocountry_2
randcoeff[["min_invest_2"]] = min_invest_2_mu + min_invest_2_sig*draws_min_invest_2
randcoeff[["duration_I_2"]] = duration_I_2_mu + duration_I_2_sig*draws_duration_I_2
randcoeff[["duration_II_2"]] = duration_II_2_mu + duration_II_2_sig*draws_duration_II_2
randcoeff[["duration_V_2"]] = duration_V_2_mu + duration_V_2_sig*draws_duration_V_2
randcoeff[["Q_participation_2"]] = A_participation_2_mu + A_participation_2_sig*draws_Q_participation_2
randcoeff[["A_participation_2"]] = A_participation_2_mu + A_participation_2_sig*draws_A_participation_2
randcoeff[["e_return_2"]] = e_return_2_mu + e_return_2_sig*draws_Ereturn_2
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["asc_III"]] = list(asc_III_1, asc_III_2,asc_III_3)
lcpars[["project"]] = list(project_1, project_2,project_3)
lcpars[["carbon"]] = list(carbon_1, carbon_2,carbon_3)
lcpars[["location_Wregion"]] = list(location_Wregion_1, location_Wregion_2,location_Wregion_3)
lcpars[["location_Wcountry"]] = list(location_Wcountry_1, location_Wcountry_2,location_Wcountry_3)
lcpars[["location_Ocountry"]] = list(location_Ocountry_1, location_Ocountry_2,location_Ocountry_3)
lcpars[["min_invest"]] = list(min_invest_1, min_invest_2,min_invest_3)
lcpars[["duration_I"]] = list(duration_I_1, duration_I_2,duration_I_3)
lcpars[["duration_II"]] = list(duration_II_1, duration_II_2,duration_II_3)
lcpars[["duration_V"]] = list(duration_V_1, duration_V_2,duration_V_3)
lcpars[["Q_participation"]] = list(Q_participation_1, Q_participation_2,Q_participation_3)
lcpars[["A_participation"]] = list(A_participation_1, A_participation_2,A_participation_3)
lcpars[["e_return"]] = list(e_return_1, e_return_2,e_return_3)
V=list()
V[["class_a"]] = delta_a
V[["class_b"]] = delta_b
V[["class_c"]] = delta_c #more classes need change
mnl_settings = list(
alternatives = c(class_a=1, class_b=2,class_c=3), #more classes need change
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()
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice
)
### Loop over classes
s=1
while(s<=3){ #change when there are more classes
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_I + project[[s]]*project1 + carbon[[s]]*rescaled_Land_co21 +
location_Wregion[[s]]*location_region_1 + location_Wcountry[[s]]*location_within_country_1 +
location_Ocountry[[s]]*location_out_country_1 + min_invest[[s]]*rescaled_min_invest1 +
duration_I[[s]]*duration_one_1 + duration_II[[s]]*duration_two_1 +
duration_V[[s]]*duration_five_1 + Q_participation[[s]]*participation_quarter_1 +
A_participation[[s]]*participation_annual_1 +
e_return[[s]]*e_return1
V[['alt2']] = asc_I + project[[s]]*project2 + carbon[[s]]*rescaled_Land_co22 +
location_Wregion[[s]]*location_region_2 + location_Wcountry[[s]]*location_within_country_2 +
location_Ocountry[[s]]*location_out_country_2 + min_invest[[s]]*rescaled_min_invest2 +
duration_I[[s]]*duration_one_2 + duration_II[[s]]*duration_two_2 +
duration_V[[s]]*duration_five_2 + Q_participation[[s]]*participation_quarter_2 +
A_participation[[s]]*participation_annual_2 +
e_return[[s]]*e_return2
V[['alt3']] = asc_III[[s]]
mnl_settings$V = V
### Compute within-class choice probabilities using MNL model
P[[s]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P[[s]] = apollo_panelProd(P[[s]], apollo_inputs ,functionality)
### Average across inter-individual draws within classes
P[[s]] = apollo_avgInterDraws(P[[s]], 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)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
Error in apollo_avgInterDraws(P[[s]], apollo_inputs, functionality) :
No Inter-individuals draws to average over!
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Thank you very much for inventing the Apollo package and creating the forum allowing our users to communicate.
I recently run a latent class model, but allow for inter-individual heterogeneity only within one of the classes (i.e., preferences are assumed to be non-random in the other classes). The model ends up with an error saying
"Error in apollo_avgInterDraws(P[[s]], apollo_inputs, functionality) :
No Inter-individuals draws to average over! "
I think the problem is the following: R tries to average across inter-individual draws for class, but the problem occurs when estimating the two non-random classes where no draws are generated.
My questions are:
(1) Is this model theoretically valid and/or practically sensible?
(2) If yes to (1), How can I code the model such as this (where let's say, attributes follow are random in one class, and non-random in the other two classes)?
I have attached my R code below.
Thank you very much.
Tim
-------------------------------------------
# ################################################################# #
#### 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 ="RPL_coop_LC3_ENVgroupRPL",
modelDescr ="RPL_coop_LC3_ENVgroupRPL",
indivID ="pool_ID",
nCores = 3,
mixing = TRUE,
noDiagnostics = TRUE
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
#load the data
database = read.csv("R_for_coop.csv",header=TRUE)
##renames##
#names(database)[names(database)=="id"] <- "ID"
database= subset(database, database$dce_valid_CB1==1&dce_valid_CB2==1&
dce_valid_CB3==1&dce_valid_CB4==1&dce_valid_T2B1==1
&dce_valid_T2B2==1&dce_valid_T2B3==1&dce_valid_T2B4==1)
############################T1############################################
####T1 and control only
database= subset(database, database$T2==0)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_I = 0,
asc_III_1 = 0,
asc_III_2_mu = 0,
asc_III_2_sig = 0,
asc_III_3 = 0,
project_1 = 0,
project_2_mu = 0,
project_2_sig = 0,
project_3 = 0,
carbon_1 = 0,
carbon_2_mu = 0,
carbon_2_sig = 0,
carbon_3 = 0,
location_Wregion_1 = 0,
location_Wregion_2_mu = 0,
location_Wregion_2_sig = 0,
location_Wregion_3 = 0,
location_Wcountry_1 = 0,
location_Wcountry_2_mu = 0,
location_Wcountry_2_sig = 0,
location_Wcountry_3 = 0,
location_Ocountry_1 = 0,
location_Ocountry_2_mu = 0,
location_Ocountry_2_sig = 0,
location_Ocountry_3 = 0,
min_invest_1 = 0,
min_invest_2_mu = 0,
min_invest_2_sig = 0,
min_invest_3 = 0,
duration_I_1 = 0,
duration_I_2_mu = 0,
duration_I_2_sig = 0,
duration_I_3 = 0,
duration_II_1 = 0,
duration_II_2_mu = 0,
duration_II_2_sig = 0,
duration_II_3 = 0,
duration_V_1 = 0,
duration_V_2_mu = 0,
duration_V_2_sig = 0,
duration_V_3 = 0,
Q_participation_1 = 0,
Q_participation_2_mu = 0,
Q_participation_2_sig = 0,
Q_participation_3 = 0,
A_participation_1 = 0,
A_participation_2_mu = 0,
A_participation_2_sig = 0,
A_participation_3 = 0,
e_return_1 = 0,
e_return_2_mu = 0,
e_return_2_sig = 0,
e_return_3 = 0,
delta_a = 0,
delta_b = 0,
delta_c = 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("delta_b", "asc_I") #fix the scale of the control treatment
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("draws_project_2",
"draws_carbon_2",
"draws_location_Wregion_2",
"draws_location_Wcountry_2",
"draws_location_Ocountry_2",
"draws_min_invest_2",
"draws_duration_I_2",
"draws_duration_II_2",
"draws_duration_V_2",
"draws_Q_participation_2",
"draws_A_participation_2",
"draws_asc_2",
"draws_Ereturn_2"
),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["asc_III_2"]] = asc_III_2_mu + asc_III_2_sig*draws_asc_2
randcoeff[["project_2"]] = project_2_mu + project_2_sig*draws_project_2
randcoeff[["carbon_2"]] = carbon_2_mu + carbon_2_sig*draws_carbon_2
randcoeff[["location_Wregion_2"]] = location_Wregion_2_mu + location_Wregion_2_sig*draws_location_Wregion_2
randcoeff[["location_Wcountry_2"]] = location_Wcountry_2_mu + location_Wcountry_2_sig*draws_location_Wcountry_2
randcoeff[["location_Ocountry_2"]] = location_Ocountry_2_mu + location_Ocountry_2_sig*draws_location_Ocountry_2
randcoeff[["min_invest_2"]] = min_invest_2_mu + min_invest_2_sig*draws_min_invest_2
randcoeff[["duration_I_2"]] = duration_I_2_mu + duration_I_2_sig*draws_duration_I_2
randcoeff[["duration_II_2"]] = duration_II_2_mu + duration_II_2_sig*draws_duration_II_2
randcoeff[["duration_V_2"]] = duration_V_2_mu + duration_V_2_sig*draws_duration_V_2
randcoeff[["Q_participation_2"]] = A_participation_2_mu + A_participation_2_sig*draws_Q_participation_2
randcoeff[["A_participation_2"]] = A_participation_2_mu + A_participation_2_sig*draws_A_participation_2
randcoeff[["e_return_2"]] = e_return_2_mu + e_return_2_sig*draws_Ereturn_2
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["asc_III"]] = list(asc_III_1, asc_III_2,asc_III_3)
lcpars[["project"]] = list(project_1, project_2,project_3)
lcpars[["carbon"]] = list(carbon_1, carbon_2,carbon_3)
lcpars[["location_Wregion"]] = list(location_Wregion_1, location_Wregion_2,location_Wregion_3)
lcpars[["location_Wcountry"]] = list(location_Wcountry_1, location_Wcountry_2,location_Wcountry_3)
lcpars[["location_Ocountry"]] = list(location_Ocountry_1, location_Ocountry_2,location_Ocountry_3)
lcpars[["min_invest"]] = list(min_invest_1, min_invest_2,min_invest_3)
lcpars[["duration_I"]] = list(duration_I_1, duration_I_2,duration_I_3)
lcpars[["duration_II"]] = list(duration_II_1, duration_II_2,duration_II_3)
lcpars[["duration_V"]] = list(duration_V_1, duration_V_2,duration_V_3)
lcpars[["Q_participation"]] = list(Q_participation_1, Q_participation_2,Q_participation_3)
lcpars[["A_participation"]] = list(A_participation_1, A_participation_2,A_participation_3)
lcpars[["e_return"]] = list(e_return_1, e_return_2,e_return_3)
V=list()
V[["class_a"]] = delta_a
V[["class_b"]] = delta_b
V[["class_c"]] = delta_c #more classes need change
mnl_settings = list(
alternatives = c(class_a=1, class_b=2,class_c=3), #more classes need change
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()
### Define settings for MNL model component that are generic across classes
mnl_settings = list(
alternatives = c(alt1=1, alt2=2,alt3=3),
avail = list(alt1=1, alt2=1,alt3=1),
choiceVar = choice
)
### Loop over classes
s=1
while(s<=3){ #change when there are more classes
### Compute class-specific utilities
V=list()
V[['alt1']] = asc_I + project[[s]]*project1 + carbon[[s]]*rescaled_Land_co21 +
location_Wregion[[s]]*location_region_1 + location_Wcountry[[s]]*location_within_country_1 +
location_Ocountry[[s]]*location_out_country_1 + min_invest[[s]]*rescaled_min_invest1 +
duration_I[[s]]*duration_one_1 + duration_II[[s]]*duration_two_1 +
duration_V[[s]]*duration_five_1 + Q_participation[[s]]*participation_quarter_1 +
A_participation[[s]]*participation_annual_1 +
e_return[[s]]*e_return1
V[['alt2']] = asc_I + project[[s]]*project2 + carbon[[s]]*rescaled_Land_co22 +
location_Wregion[[s]]*location_region_2 + location_Wcountry[[s]]*location_within_country_2 +
location_Ocountry[[s]]*location_out_country_2 + min_invest[[s]]*rescaled_min_invest2 +
duration_I[[s]]*duration_one_2 + duration_II[[s]]*duration_two_2 +
duration_V[[s]]*duration_five_2 + Q_participation[[s]]*participation_quarter_2 +
A_participation[[s]]*participation_annual_2 +
e_return[[s]]*e_return2
V[['alt3']] = asc_III[[s]]
mnl_settings$V = V
### Compute within-class choice probabilities using MNL model
P[[s]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P[[s]] = apollo_panelProd(P[[s]], apollo_inputs ,functionality)
### Average across inter-individual draws within classes
P[[s]] = apollo_avgInterDraws(P[[s]], 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)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
Error in apollo_avgInterDraws(P[[s]], apollo_inputs, functionality) :
No Inter-individuals draws to average over!
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)