sorry for my late response. I had not received an email that a reply to the forum entry had been created - or I deleted it unintentionally. Anyway: Yes, I would like the class allocation model to be a function of the LV.
But first a few words about the model: In addition to respondents' preferences (using the choice model), we would like to investigate why respondents ignored certain attributes or considered them less important in their choices (latent variable model). To do this, we queried three latent variables per attribute:
(1) The attribute is too complex to understand.
(2) It is not worth the effort (e.g., in terms of time) to think about the impact the attribute has on a respondent's daily life.
(3) The attribute is in fact irrelevant.
Factor analysis showed that the first two LVs cannot be separated, which is why we merged them. That leaves us with two latent variables:
We would like to specify the latent variable model with two stages, i.e. it should look like this:
Non-attendance is to be used for class allocation.
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
library(corrplot)
library(Hmisc)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "flexQgrid_ICLV_LC_3Classes_DR1_CutOff6_SO2000",
modelDescr = "Latent class model on flexQgrid panel data (classical estimation)",
indivID ="sys_RespNum",
nCores = 60, # parallel::detectCores()
HB = FALSE,
mixing = TRUE,
workInLogs = FALSE # TRUE for higher numeric stability at the expense of computational time; works with the sum of log-probabilities, rather than the product of probabilities
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
input.file <- "./input/CBC_results_merged.wide.rds"
database <- readRDS(input.file)
### Dual response
database[which(database$CBC_type == "RD" & database$Choice <6), "Choice"] <- 0
database[which(database$CBC_type == "RD" & database$Choice >=6), "Choice"] <- 1
### Dummy for first row of a respondent
database[, "sys_RespNum.first"] <- rep(0, times=nrow(database))
database[match( unique(database$sys_RespNum), database$sys_RespNum), "sys_RespNum.first"] <- 1
### Mean centering of indicators
NonAttendance_colIndices <- grep("NonAttendance", colnames(database))
database[,NonAttendance_colIndices] <- apply(database[,NonAttendance_colIndices], 2, scale, center=TRUE, scale=FALSE)
### Load priors from step 1
# priors.step1 = list()
# priors.step1[["FC"]] = readRDS(file = "./priors_step2/FC_step1.rds") # A vector of starting values for the fixed parameters
# priors.step1[["svN"]] = readRDS(file = "./priors_step2/svN_step1.rds") # A vector of starting values for the means of the underlying normals for the random parameters
# priors.step1[["pvMatrix"]] = readRDS(file = "./priors_step2/pvMatrix_step1.rds") # A custom prior variance-covariance matrix to be used in estimation
# ################################################################# #
#### PRE-PROCESSING DATA ####
# ################################################################# #
covariates.colnames <- c(
"COV_Gender",
"COV_Age_centered",
"COV_Education_log_centered",
"COV_IncomeClasses_log_centered",
"COV_HouseholdSize_centered",
"COV_ApartmentOwner",
"COV_ExistingApps",
"COV_CurrentMix",
"COV_EnergyDecisionMaker",
"COV_CheapTalk"
)
covariates <- database[,c("sys_RespNum",covariates.colnames)]
covariates <- covariates[!duplicated(covariates$sys_RespNum),] # only keep one row per respondent
covariates <- subset(covariates, select=-c(sys_RespNum))
## Corrplot
covariates.correlation = cor(x = covariates, method="pearson")
corrplot(covariates.correlation, method = "number")
covariates.correlation.p = rcorr(as.matrix(covariates), type="pearson") # "pearson" or "spearman"
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
# ----------------------------------------------------------------- #
#---- Generic parameters
# ----------------------------------------------------------------- #
### Alternative specific constants (ASC)
b_asc_1 = 0,
b_asc_2 = 0,
b_asc_3 = 0,
# ----------------------------------------------------------------- #
#---- Class Assignment (A)
# ----------------------------------------------------------------- #
delta_A = 0,
### Sociodemographics
# gamma_A_Gender = 0,
# gamma_A_Age = 0,
# gamma_A_Education = 0,
# gamma_A_IncomeClasses = 0,
#
# gamma_A_ExistingApps = 0,
# gamma_A_ApartmentOwner = 0,
### Influence of LV on class assignment
gamma_A_NACount = 0,
gamma_A_NADuration = 0,
gamma_A_NAPeriod = 0,
gamma_A_NACompensation = 0,
# ----------------------------------------------------------------- #
#---- Class Assignment (B)
# ----------------------------------------------------------------- #
delta_B = 0.1,
### Sociodemographics
# gamma_B_Gender = 0,
# gamma_B_Age = 0,
# gamma_B_Education = 0,
# gamma_B_IncomeClasses = 0,
#
# gamma_B_ExistingApps = 0,
# gamma_B_ApartmentOwner = 0,
### Influence of LV on class assignment
gamma_B_NACount = 0,
gamma_B_NADuration = 0,
gamma_B_NAPeriod = 0,
gamma_B_NACompensation = 0,
# ----------------------------------------------------------------- #
#---- Class Assignment (C)
# ----------------------------------------------------------------- #
delta_C = 0.1,
### Sociodemographics
# gamma_C_Gender = 0,
# gamma_C_Age = 0,
# gamma_C_Education = 0,
# gamma_C_IncomeClasses = 0,
#
# gamma_C_ExistingApps = 0,
# gamma_C_ApartmentOwner = 0,
### Influence of LV on class assignment
gamma_C_NACount = 0,
gamma_C_NADuration = 0,
gamma_C_NAPeriod = 0,
gamma_C_NACompensation = 0,
# ----------------------------------------------------------------- #
#---- Intra-Class-Parameters (A)
# ----------------------------------------------------------------- #
### DR_asc
A_b_DR_asc_0 = 0,
A_b_DR_asc_1 = 0,
### Quota_Count
A_b_Quota_Count = 0,
### Quota_Max_Duration
A_b_Quota_Max_Duration = 0,
### Quota_Period_midmorning
A_b_Quota_Period_midmorning = 0,
### Quota_Period_afternoon
A_b_Quota_Period_afternoon = 0,
### Quota_Period_evening
A_b_Quota_Period_evening = 0,
### Quota_Period_night
A_b_Quota_Period_night = 0,
### Quota_Compensation
A_b_Quota_Compensation = -0.1,
# ----------------------------------------------------------------- #
#---- Intra-Class-Parameters (B)
# ----------------------------------------------------------------- #
### DR_asc
B_b_DR_asc_0 = 0,
B_b_DR_asc_1 = 0,
### Quota_Count
B_b_Quota_Count = 0,
### Quota_Max_Duration
B_b_Quota_Max_Duration = 0,
### Quota_Period_midmorning
B_b_Quota_Period_midmorning = 0,
### Quota_Period_afternoon
B_b_Quota_Period_afternoon = 0,
### Quota_Period_evening
B_b_Quota_Period_evening = 0,
### Quota_Period_night
B_b_Quota_Period_night = 0,
### Quota_Compensation
B_b_Quota_Compensation = -0.1,
# ----------------------------------------------------------------- #
#---- Intra-Class-Parameters (C)
# ----------------------------------------------------------------- #
### DR_asc
C_b_DR_asc_0 = 0,
C_b_DR_asc_1 = 0,
### Quota_Count
C_b_Quota_Count = 0,
### Quota_Max_Duration
C_b_Quota_Max_Duration = 0,
### Quota_Period_midmorning
C_b_Quota_Period_midmorning = 0,
### Quota_Period_afternoon
C_b_Quota_Period_afternoon = 0,
### Quota_Period_evening
C_b_Quota_Period_evening = 0,
### Quota_Period_night
C_b_Quota_Period_night = 0,
### Quota_Compensation
C_b_Quota_Compensation = -0.2,
# ----------------------------------------------------------------- #
#---- Latent Variable Model
# ----------------------------------------------------------------- #
### Influence of LV on other LV
lambda_NACountComplexityCosts_NACount = 1,
lambda_NACountIrrelevance_NACount = 1,
lambda_NADurationComplexityCosts_NADuration = 1,
lambda_NADurationIrrelevance_NADuration = 1,
lambda_NAPeriodComplexityCosts_NAPeriod = 1,
lambda_NAPeriodIrrelevance_NAPeriod = 1,
lambda_NACompensationComplexityCosts_NACompensation = 1,
lambda_NACompensationIrrelevance_NACompensation = 1,
### Continuous measurement model
sigma_NACount_r1 = 1,
sigma_NACount_r2 = 1,
sigma_NACount_r3 = 1,
sigma_NACount_r4 = 1,
sigma_NACount_r5 = 1,
sigma_NACount_r6 = 1,
sigma_NADuration_r1 = 1,
sigma_NADuration_r2 = 1,
sigma_NADuration_r3 = 1,
sigma_NADuration_r4 = 1,
sigma_NADuration_r5 = 1,
sigma_NADuration_r6 = 1,
sigma_NAPeriod_r1 = 1,
sigma_NAPeriod_r2 = 1,
sigma_NAPeriod_r3 = 1,
sigma_NAPeriod_r4 = 1,
sigma_NAPeriod_r5 = 1,
sigma_NAPeriod_r6 = 1,
sigma_NACompensation_r1 = 1,
sigma_NACompensation_r2 = 1,
sigma_NACompensation_r3 = 1,
sigma_NACompensation_r4 = 1,
sigma_NACompensation_r5 = 1,
sigma_NACompensation_r6 = 1,
### Marketing scales
zeta_NACountComplexityCosts_r1_NACount= 1,
zeta_NACountComplexityCosts_r2_NACount = 1,
zeta_NACountComplexityCosts_r3_NACount = 1,
zeta_NACountComplexityCosts_r4_NACount = 1,
zeta_NACountIrrelevance_r5_NACount = 1,
zeta_NACountIrrelevance_r6_NACount = 1,
zeta_NADurationComplexityCosts_r1_NADuration= 1,
zeta_NADurationComplexityCosts_r2_NADuration = 1,
zeta_NADurationComplexityCosts_r3_NADuration = 1,
zeta_NADurationComplexityCosts_r4_NADuration = 1,
zeta_NADurationIrrelevance_r5_NADuration = 1,
zeta_NADurationIrrelevance_r6_NADuration = 1,
zeta_NAPeriodComplexityCosts_r1_NAPeriod= 1,
zeta_NAPeriodComplexityCosts_r2_NAPeriod = 1,
zeta_NAPeriodComplexityCosts_r3_NAPeriod = 1,
zeta_NAPeriodComplexityCosts_r4_NAPeriod = 1,
zeta_NAPeriodIrrelevance_r5_NAPeriod = 1,
zeta_NAPeriodIrrelevance_r6_NAPeriod = 1,
zeta_NACompensationComplexityCosts_r1_NACompensation= 1,
zeta_NACompensationComplexityCosts_r2_NACompensation = 1,
zeta_NACompensationComplexityCosts_r3_NACompensation = 1,
zeta_NACompensationComplexityCosts_r4_NACompensation = 1,
zeta_NACompensationIrrelevance_r5_NACompensation = 1,
zeta_NACompensationIrrelevance_r6_NACompensation = 1,
### Sociodemographics
gamma_Age_NACountComplexityCosts = 0,
gamma_Education_NACountComplexityCosts = 0,
gamma_ExistingApps_NACountComplexityCosts = 0,
sigma_NACountComplexityCosts = 1,
gamma_Age_NACountIrrelevance = 0,
gamma_Education_NACountIrrelevance = 0,
gamma_ExistingApps_NACountIrrelevance = 0,
sigma_NACountIrrelevance = 1,
gamma_Age_NADurationComplexityCosts = 0,
gamma_Education_NADurationComplexityCosts = 0,
gamma_ExistingApps_NADurationComplexityCosts = 0,
sigma_NADurationComplexityCosts = 1,
gamma_Age_NADurationIrrelevance = 0,
gamma_Education_NADurationIrrelevance = 0,
gamma_ExistingApps_NADurationIrrelevance = 0,
sigma_NADurationIrrelevance = 1,
gamma_Age_NAPeriodComplexityCosts = 0,
gamma_Education_NAPeriodComplexityCosts = 0,
gamma_ExistingApps_NAPeriodComplexityCosts = 0,
sigma_NAPeriodComplexityCosts = 1,
gamma_Age_NAPeriodIrrelevance = 0,
gamma_Education_NAPeriodIrrelevance = 0,
gamma_ExistingApps_NAPeriodIrrelevance = 0,
sigma_NAPeriodIrrelevance = 1,
gamma_Gender_NACompensationComplexityCosts = 0,
gamma_Age_NACompensationComplexityCosts = 0,
gamma_Education_NACompensationComplexityCosts = 0,
gamma_Income_NACompensationComplexityCosts = 0,
gamma_CheapTalk_NACompensationComplexityCosts = 0,
sigma_NACompensationComplexityCosts = 1,
gamma_Gender_NACompensationIrrelevance = 0,
gamma_Age_NACompensationIrrelevance = 0,
gamma_Education_NACompensationIrrelevance = 0,
gamma_Income_NACompensationIrrelevance = 0,
gamma_CheapTalk_NACompensationIrrelevance = 0,
sigma_NACompensationIrrelevance = 1
)
### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_bwtp_fixed = c() if none
apollo_fixed = c("b_asc_3",
"delta_A",
"gamma_A_NACount",
"gamma_A_NADuration",
"gamma_A_NAPeriod",
"gamma_A_NACompensation",
"lambda_NACountComplexityCosts_NACount",
"lambda_NADurationComplexityCosts_NADuration",
"lambda_NAPeriodComplexityCosts_NAPeriod",
"lambda_NACompensationComplexityCosts_NACompensation",
"zeta_NACountComplexityCosts_r1_NACount", # fix most reliable item if possible/known
"zeta_NACountIrrelevance_r5_NACount",
"zeta_NADurationComplexityCosts_r1_NADuration", # fix most reliable item if possible/known
"zeta_NADurationIrrelevance_r5_NADuration",
"zeta_NAPeriodComplexityCosts_r1_NAPeriod", # fix most reliable item if possible/known
"zeta_NAPeriodIrrelevance_r5_NAPeriod",
"zeta_NACompensationComplexityCosts_r1_NACompensation", # fix most reliable item if possible/known
"zeta_NACompensationIrrelevance_r5_NACompensation",
"A_b_DR_asc_0",
"B_b_DR_asc_0",
"C_b_DR_asc_0"
)
### Read estimates from other models
apollo_beta = apollo_readBeta(apollo_beta, apollo_fixed, "flexQgrid_LC_3Classes_DR1_CutOff6", overwriteFixed=FALSE)
# apollo_beta_new = replace(apollo_beta, names(priors.step1$FC), priors.step1$FC)
# apollo_beta_new = replace(apollo_beta_new, names(priors.step1$svN), priors.step1$svN)
# apollo_beta = apollo_beta_new[names(apollo_beta)]
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "sobolOwen",
interNDraws = 2000,
interUnifDraws = c(),
interNormDraws = c("eta_NACountComplexityCosts", "eta_NACountIrrelevance",
"eta_NADurationComplexityCosts", "eta_NADurationIrrelevance",
"eta_NAPeriodComplexityCosts", "eta_NAPeriodIrrelevance",
"eta_NACompensationComplexityCosts", "eta_NACompensationIrrelevance"),
intraDrawsType = "sobolOwen",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs) {
randcoeff = list()
# ----------------------------------------------------------------- #
#---- Latent Variables
# ----------------------------------------------------------------- #
randcoeff[["LV_NACountComplexityCosts"]] = gamma_Age_NACountComplexityCosts * COV_Age_centered +
gamma_Education_NACountComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NACountComplexityCosts * COV_ExistingApps +
sigma_NACountComplexityCosts * eta_NACountComplexityCosts
randcoeff[["LV_NACountIrrelevance"]] = gamma_Age_NACountIrrelevance * COV_Age_centered +
gamma_Education_NACountIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NACountIrrelevance * COV_ExistingApps +
sigma_NACountIrrelevance * eta_NACountIrrelevance
randcoeff[["LV_NACount"]] = lambda_NACountComplexityCosts_NACount *
( # LV_NACountComplexityCosts
gamma_Age_NACountComplexityCosts * COV_Age_centered +
gamma_Education_NACountComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NACountComplexityCosts * COV_ExistingApps +
sigma_NACountComplexityCosts * eta_NACountComplexityCosts
) +
lambda_NACountIrrelevance_NACount *
( # LV_NACountIrrelevance
gamma_Age_NACountIrrelevance * COV_Age_centered +
gamma_Education_NACountIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NACountIrrelevance * COV_ExistingApps +
sigma_NACountIrrelevance * eta_NACountIrrelevance
)
randcoeff[["LV_NADurationComplexityCosts"]] = gamma_Age_NADurationComplexityCosts * COV_Age_centered +
gamma_Education_NADurationComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NADurationComplexityCosts * COV_ExistingApps +
sigma_NADurationComplexityCosts * eta_NADurationComplexityCosts
randcoeff[["LV_NADurationIrrelevance"]] = gamma_Age_NADurationIrrelevance * COV_Age_centered +
gamma_Education_NADurationIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NADurationIrrelevance * COV_ExistingApps +
sigma_NADurationIrrelevance * eta_NADurationIrrelevance
randcoeff[["LV_NADuration"]] = lambda_NADurationComplexityCosts_NADuration *
( # LV_NADurationComplexityCosts
gamma_Age_NADurationComplexityCosts * COV_Age_centered +
gamma_Education_NADurationComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NADurationComplexityCosts * COV_ExistingApps +
sigma_NADurationComplexityCosts * eta_NADurationComplexityCosts
) +
lambda_NADurationIrrelevance_NADuration *
( # LV_NADurationIrrelevance
gamma_Age_NADurationIrrelevance * COV_Age_centered +
gamma_Education_NADurationIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NADurationIrrelevance * COV_ExistingApps +
sigma_NADurationIrrelevance * eta_NADurationIrrelevance
)
randcoeff[["LV_NAPeriodComplexityCosts"]] = gamma_Age_NAPeriodComplexityCosts * COV_Age_centered +
gamma_Education_NAPeriodComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NAPeriodComplexityCosts * COV_ExistingApps +
sigma_NAPeriodComplexityCosts * eta_NAPeriodComplexityCosts
randcoeff[["LV_NAPeriodIrrelevance"]] = gamma_Age_NAPeriodIrrelevance * COV_Age_centered +
gamma_Education_NAPeriodIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NAPeriodIrrelevance * COV_ExistingApps +
sigma_NAPeriodIrrelevance * eta_NAPeriodIrrelevance
randcoeff[["LV_NAPeriod"]] = lambda_NAPeriodComplexityCosts_NAPeriod *
( # LV_NAPeriodComplexityCosts
gamma_Age_NAPeriodComplexityCosts * COV_Age_centered +
gamma_Education_NAPeriodComplexityCosts * COV_Education_log_centered +
gamma_ExistingApps_NAPeriodComplexityCosts * COV_ExistingApps +
sigma_NAPeriodComplexityCosts * eta_NAPeriodComplexityCosts
) +
lambda_NAPeriodIrrelevance_NAPeriod *
( # LV_NAPeriodIrrelevance
gamma_Age_NAPeriodIrrelevance * COV_Age_centered +
gamma_Education_NAPeriodIrrelevance * COV_Education_log_centered +
gamma_ExistingApps_NAPeriodIrrelevance * COV_ExistingApps +
sigma_NAPeriodIrrelevance * eta_NAPeriodIrrelevance
)
randcoeff[["LV_NACompensationComplexityCosts"]] = gamma_Gender_NACompensationComplexityCosts * COV_Gender +
gamma_Age_NACompensationComplexityCosts * COV_Age_centered +
gamma_Education_NACompensationComplexityCosts * COV_Education_log_centered +
(1-IncomeQuan_miss) * gamma_Income_NACompensationComplexityCosts * COV_IncomeClasses_log_centered +
gamma_CheapTalk_NACompensationComplexityCosts * COV_CheapTalk +
sigma_NACompensationComplexityCosts * eta_NACompensationComplexityCosts
randcoeff[["LV_NACompensationIrrelevance"]] = gamma_Gender_NACompensationIrrelevance * COV_Gender +
gamma_Age_NACompensationIrrelevance * COV_Age_centered +
gamma_Education_NACompensationIrrelevance * COV_Education_log_centered +
(1-IncomeQuan_miss) * gamma_Income_NACompensationIrrelevance * COV_IncomeClasses_log_centered +
gamma_CheapTalk_NACompensationIrrelevance * COV_CheapTalk +
sigma_NACompensationIrrelevance * eta_NACompensationIrrelevance
randcoeff[["LV_NACompensation"]] = lambda_NACompensationComplexityCosts_NACompensation *
( # LV_NACompensationComplexityCosts
gamma_Gender_NACompensationComplexityCosts * COV_Gender +
gamma_Age_NACompensationComplexityCosts * COV_Age_centered +
gamma_Education_NACompensationComplexityCosts * COV_Education_log_centered +
(1-IncomeQuan_miss) * gamma_Income_NACompensationComplexityCosts * COV_IncomeClasses_log_centered +
gamma_CheapTalk_NACompensationComplexityCosts * COV_CheapTalk +
sigma_NACompensationComplexityCosts * eta_NACompensationComplexityCosts
) +
lambda_NACompensationIrrelevance_NACompensation *
( # LV_NACompensationIrrelevance
gamma_Gender_NACompensationIrrelevance * COV_Gender +
gamma_Age_NACompensationIrrelevance * COV_Age_centered +
gamma_Education_NACompensationIrrelevance * COV_Education_log_centered +
(1-IncomeQuan_miss) * gamma_Income_NACompensationIrrelevance * COV_IncomeClasses_log_centered +
gamma_CheapTalk_NACompensationIrrelevance * COV_CheapTalk +
sigma_NACompensationIrrelevance * eta_NACompensationIrrelevance
)
return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars = function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_DR_asc_0"]] = list(A_b_DR_asc_0, B_b_DR_asc_0, C_b_DR_asc_0)
lcpars[["b_DR_asc_1"]] = list(A_b_DR_asc_1, B_b_DR_asc_1, C_b_DR_asc_1)
lcpars[["b_Quota_Count"]] = list(A_b_Quota_Count, B_b_Quota_Count, C_b_Quota_Count)
lcpars[["b_Quota_Max_Duration"]] = list(A_b_Quota_Max_Duration, B_b_Quota_Max_Duration, C_b_Quota_Max_Duration)
lcpars[["b_Quota_Period_midmorning"]] = list(A_b_Quota_Period_midmorning, B_b_Quota_Period_midmorning, C_b_Quota_Period_midmorning)
lcpars[["b_Quota_Period_afternoon"]] = list(A_b_Quota_Period_afternoon, B_b_Quota_Period_afternoon, C_b_Quota_Period_afternoon)
lcpars[["b_Quota_Period_evening"]] = list(A_b_Quota_Period_evening, B_b_Quota_Period_evening, C_b_Quota_Period_evening)
lcpars[["b_Quota_Period_night"]] = list(A_b_Quota_Period_night, B_b_Quota_Period_night, C_b_Quota_Period_night)
lcpars[["b_Quota_Compensation"]] = list(A_b_Quota_Compensation, B_b_Quota_Compensation, C_b_Quota_Compensation)
V=list()
V[["class_A"]] = delta_A +
gamma_A_NACount * LV_NACount +
gamma_A_NADuration * LV_NADuration +
gamma_A_NAPeriod * LV_NAPeriod +
gamma_A_NACompensation * LV_NACompensation
V[["class_B"]] = delta_B +
gamma_B_NACount * LV_NACount +
gamma_B_NADuration * LV_NADuration +
gamma_B_NAPeriod * LV_NAPeriod +
gamma_B_NACompensation * LV_NACompensation
V[["class_C"]] = delta_C +
gamma_C_NACount * LV_NACount +
gamma_C_NADuration * LV_NADuration +
gamma_C_NAPeriod * LV_NAPeriod +
gamma_C_NACompensation * LV_NACompensation
mnl_settings = list(
alternatives = c(class_A=1, class_B=2, class_C=3),
avail = 1,
choiceVar = NA,
V = V
)
lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
lcpars[["pi_values"]] = apollo_firstRow(lcpars[["pi_values"]], apollo_inputs)
return(lcpars)
}
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate") {
### Attach inputs and detach after function exit
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
### Create list of probabilities P
P <- list()
# ----------------------------------------------------------------- #
#---- Likelihood of indicators
# ----------------------------------------------------------------- #
### NACount
normalDensity_settings1 = list(outcomeNormal = NonAttendance_r1,
xNormal = zeta_NACountComplexityCosts_r1_NACount * LV_NACountComplexityCosts,
mu = 0,
sigma = sigma_NACount_r1,
rows=(sys_RespNum.first==1),
componentName = "NACount_r1_NACount")
normalDensity_settings2 = list(outcomeNormal = NonAttendance_r2,
xNormal = zeta_NACountComplexityCosts_r2_NACount * LV_NACountComplexityCosts,
mu = 0,
sigma = sigma_NACount_r2,
rows=(sys_RespNum.first==1),
componentName = "NACount_r2_NACount")
normalDensity_settings3 = list(outcomeNormal = NonAttendance_r3,
xNormal = zeta_NACountComplexityCosts_r3_NACount * LV_NACountComplexityCosts,
mu = 0,
sigma = sigma_NACount_r3,
rows=(sys_RespNum.first==1),
componentName = "NACount_r3_NACount")
normalDensity_settings4 = list(outcomeNormal = NonAttendance_r4,
xNormal = zeta_NACountComplexityCosts_r4_NACount * LV_NACountComplexityCosts,
mu = 0,
sigma = sigma_NACount_r4,
rows=(sys_RespNum.first==1),
componentName = "NACount_r4_NACount")
normalDensity_settings5 = list(outcomeNormal = NonAttendance_r5,
xNormal = zeta_NACountIrrelevance_r5_NACount * LV_NACountIrrelevance,
mu = 0,
sigma = sigma_NACount_r5,
rows=(sys_RespNum.first==1),
componentName = "NACount_r5_NACount")
normalDensity_settings6 = list(outcomeNormal = NonAttendance_r6,
xNormal = zeta_NACountIrrelevance_r6_NACount * LV_NACountIrrelevance,
mu = 0,
sigma = sigma_NACount_r6,
rows=(sys_RespNum.first==1),
componentName = "NACount_r6_NACount")
### NADuration
normalDensity_settings7 = list(outcomeNormal = NonAttendance2_r1,
xNormal = zeta_NADurationComplexityCosts_r1_NADuration * LV_NADurationComplexityCosts,
mu = 0,
sigma = sigma_NADuration_r1,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r1_NADuration")
normalDensity_settings8 = list(outcomeNormal = NonAttendance2_r2,
xNormal = zeta_NADurationComplexityCosts_r2_NADuration * LV_NADurationComplexityCosts,
mu = 0,
sigma = sigma_NADuration_r2,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r2_NADuration")
normalDensity_settings9 = list(outcomeNormal = NonAttendance2_r3,
xNormal = zeta_NADurationComplexityCosts_r3_NADuration * LV_NADurationComplexityCosts,
mu = 0,
sigma = sigma_NADuration_r3,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r3_NADuration")
normalDensity_settings10 = list(outcomeNormal = NonAttendance2_r4,
xNormal = zeta_NADurationComplexityCosts_r4_NADuration * LV_NADurationComplexityCosts,
mu = 0,
sigma = sigma_NADuration_r4,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r4_NADuration")
normalDensity_settings11 = list(outcomeNormal = NonAttendance2_r5,
xNormal = zeta_NADurationIrrelevance_r5_NADuration * LV_NADurationIrrelevance,
mu = 0,
sigma = sigma_NADuration_r5,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r5_NADuration")
normalDensity_settings12 = list(outcomeNormal = NonAttendance2_r6,
xNormal = zeta_NADurationIrrelevance_r6_NADuration * LV_NADurationIrrelevance,
mu = 0,
sigma = sigma_NADuration_r6,
rows=(sys_RespNum.first==1),
componentName = "NADuration_r6_NADuration")
### NAPeriod
normalDensity_settings13 = list(outcomeNormal = NonAttendance3_r1,
xNormal = zeta_NAPeriodComplexityCosts_r1_NAPeriod * LV_NAPeriodComplexityCosts,
mu = 0,
sigma = sigma_NAPeriod_r1,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r1_NAPeriod")
normalDensity_settings14 = list(outcomeNormal = NonAttendance3_r2,
xNormal = zeta_NAPeriodComplexityCosts_r2_NAPeriod * LV_NAPeriodComplexityCosts,
mu = 0,
sigma = sigma_NAPeriod_r2,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r2_NAPeriod")
normalDensity_settings15 = list(outcomeNormal = NonAttendance3_r3,
xNormal = zeta_NAPeriodComplexityCosts_r3_NAPeriod * LV_NAPeriodComplexityCosts,
mu = 0,
sigma = sigma_NAPeriod_r3,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r3_NAPeriod")
normalDensity_settings16 = list(outcomeNormal = NonAttendance3_r4,
xNormal = zeta_NAPeriodComplexityCosts_r4_NAPeriod * LV_NAPeriodComplexityCosts,
mu = 0,
sigma = sigma_NAPeriod_r4,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r4_NAPeriod")
normalDensity_settings17 = list(outcomeNormal = NonAttendance3_r5,
xNormal = zeta_NAPeriodIrrelevance_r5_NAPeriod * LV_NAPeriodIrrelevance,
mu = 0,
sigma = sigma_NAPeriod_r5,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r5_NAPeriod")
normalDensity_settings18 = list(outcomeNormal = NonAttendance3_r6,
xNormal = zeta_NAPeriodIrrelevance_r6_NAPeriod * LV_NAPeriodIrrelevance,
mu = 0,
sigma = sigma_NAPeriod_r6,
rows=(sys_RespNum.first==1),
componentName = "NAPeriod_r6_NAPeriod")
### NACompensation
normalDensity_settings19 = list(outcomeNormal = NonAttendance4_r1,
xNormal = zeta_NACompensationComplexityCosts_r1_NACompensation * LV_NACompensationComplexityCosts,
mu = 0,
sigma = sigma_NACompensation_r1,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r1_NACompensation")
normalDensity_settings20 = list(outcomeNormal = NonAttendance4_r2,
xNormal = zeta_NACompensationComplexityCosts_r2_NACompensation * LV_NACompensationComplexityCosts,
mu = 0,
sigma = sigma_NACompensation_r2,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r2_NACompensation")
normalDensity_settings21 = list(outcomeNormal = NonAttendance4_r3,
xNormal = zeta_NACompensationComplexityCosts_r3_NACompensation * LV_NACompensationComplexityCosts,
mu = 0,
sigma = sigma_NACompensation_r3,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r3_NACompensation")
normalDensity_settings22 = list(outcomeNormal = NonAttendance4_r4,
xNormal = zeta_NACompensationComplexityCosts_r4_NACompensation * LV_NACompensationComplexityCosts,
mu = 0,
sigma = sigma_NACompensation_r4,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r4_NACompensation")
normalDensity_settings23 = list(outcomeNormal = NonAttendance4_r5,
xNormal = zeta_NACompensationIrrelevance_r5_NACompensation * LV_NACompensationIrrelevance,
mu = 0,
sigma = sigma_NACompensation_r5,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r5_NACompensation")
normalDensity_settings24 = list(outcomeNormal = NonAttendance4_r6,
xNormal = zeta_NACompensationIrrelevance_r6_NACompensation * LV_NACompensationIrrelevance,
mu = 0,
sigma = sigma_NACompensation_r6,
rows=(sys_RespNum.first==1),
componentName = "NACompensation_r6_NACompensation")
### NACount
P[["indic_NACount_NACount_r1"]] = apollo_normalDensity(normalDensity_settings1, functionality)
P[["indic_NACount_NACount_r2"]] = apollo_normalDensity(normalDensity_settings2, functionality)
P[["indic_NACount_NACount_r3"]] = apollo_normalDensity(normalDensity_settings3, functionality)
P[["indic_NACount_NACount_r4"]] = apollo_normalDensity(normalDensity_settings4, functionality)
P[["indic_NACount_NACount_r5"]] = apollo_normalDensity(normalDensity_settings5, functionality)
P[["indic_NACount_NACount_r6"]] = apollo_normalDensity(normalDensity_settings6, functionality)
### NADuration
P[["indic_NADuration_NADuration_r1"]] = apollo_normalDensity(normalDensity_settings7, functionality)
P[["indic_NADuration_NADuration_r2"]] = apollo_normalDensity(normalDensity_settings8, functionality)
P[["indic_NADuration_NADuration_r3"]] = apollo_normalDensity(normalDensity_settings9, functionality)
P[["indic_NADuration_NADuration_r4"]] = apollo_normalDensity(normalDensity_settings10, functionality)
P[["indic_NADuration_NADuration_r5"]] = apollo_normalDensity(normalDensity_settings11, functionality)
P[["indic_NADuration_NADuration_r6"]] = apollo_normalDensity(normalDensity_settings12, functionality)
### NAPeriod
P[["indic_NAPeriod_NAPeriod_r1"]] = apollo_normalDensity(normalDensity_settings13, functionality)
P[["indic_NAPeriod_NAPeriod_r2"]] = apollo_normalDensity(normalDensity_settings14, functionality)
P[["indic_NAPeriod_NAPeriod_r3"]] = apollo_normalDensity(normalDensity_settings15, functionality)
P[["indic_NAPeriod_NAPeriod_r4"]] = apollo_normalDensity(normalDensity_settings16, functionality)
P[["indic_NAPeriod_NAPeriod_r5"]] = apollo_normalDensity(normalDensity_settings17, functionality)
P[["indic_NAPeriod_NAPeriod_r6"]] = apollo_normalDensity(normalDensity_settings18, functionality)
### NACompensation
P[["indic_NACompensation_NACompensation_r1"]] = apollo_normalDensity(normalDensity_settings19, functionality)
P[["indic_NACompensation_NACompensation_r2"]] = apollo_normalDensity(normalDensity_settings20, functionality)
P[["indic_NACompensation_NACompensation_r3"]] = apollo_normalDensity(normalDensity_settings21, functionality)
P[["indic_NACompensation_NACompensation_r4"]] = apollo_normalDensity(normalDensity_settings22, functionality)
P[["indic_NACompensation_NACompensation_r5"]] = apollo_normalDensity(normalDensity_settings23, functionality)
P[["indic_NACompensation_NACompensation_r6"]] = apollo_normalDensity(normalDensity_settings24, functionality)
### Combined model
P = apollo_combineModels(P, apollo_inputs, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
# ### Average across inter-individual draws
# P = apollo_avgInterDraws(P, apollo_inputs, functionality)
### Rename model
names(P)[which(names(P)=="model")] <- "Measurement_model"
# ----------------------------------------------------------------- #
#---- LC model
# ----------------------------------------------------------------- #
### Create list of probabilities P
P_within <- list() # list for within class models
### Loop over classes
for (s in 1:3) {
### Create list of probabilities P.componentName
P.componentName <- list()
### Compute class-specific utilities
### CBC (Forced Choices)
V = list()
V[['alt1']] = b_asc_1 +
b_Quota_Count[[s]] * Quota_Count.1 +
b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
b_Quota_Compensation[[s]] * Quota_Compensation.1
V[['alt2']] = b_asc_2 +
b_Quota_Count[[s]] * Quota_Count.2 +
b_Quota_Max_Duration[[s]] * Quota_Max_Duration.2 +
b_Quota_Period_midmorning[[s]] * Quota_Period2.2 + b_Quota_Period_afternoon[[s]] * Quota_Period3.2 + b_Quota_Period_evening[[s]] * Quota_Period4.2 + b_Quota_Period_night[[s]] * Quota_Period5.2 +
b_Quota_Compensation[[s]] * Quota_Compensation.2
V[['alt3']] = b_asc_3 +
b_Quota_Count[[s]] * Quota_Count.3 +
b_Quota_Max_Duration[[s]] * Quota_Max_Duration.3 +
b_Quota_Period_midmorning[[s]] * Quota_Period2.3 + b_Quota_Period_afternoon[[s]] * Quota_Period3.3 + b_Quota_Period_evening[[s]] * Quota_Period4.3 + b_Quota_Period_night[[s]] * Quota_Period5.3 +
b_Quota_Compensation[[s]] * Quota_Compensation.3
componentName1 = paste("CBC_", toString(s), sep="")
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3),
avail = list(alt1=1, alt2=1, alt3=1),
choiceVar = Choice,
V = V,
rows = (CBC_type=="R"),
componentName = componentName1
)
### Compute within-class choice probabilities using MNL model
P.componentName[[componentName1]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P.componentName[[componentName1]] = apollo_panelProd(P.componentName[[componentName1]], apollo_inputs, functionality)
### CBC (Free Choices)
V = list()
V[['alt0']] = b_DR_asc_0[[s]]
V[['alt1']] = b_DR_asc_1[[s]] +
b_Quota_Count[[s]] * Quota_Count.1 +
b_Quota_Max_Duration[[s]] * Quota_Max_Duration.1 +
b_Quota_Period_midmorning[[s]] * Quota_Period2.1 + b_Quota_Period_afternoon[[s]] * Quota_Period3.1 + b_Quota_Period_evening[[s]] * Quota_Period4.1 + b_Quota_Period_night[[s]] * Quota_Period5.1 +
b_Quota_Compensation[[s]] * Quota_Compensation.1
componentName2 = paste("CBC_DR_", toString(s), sep="")
mnl_settings = list(
alternatives = c(alt0=0, alt1=1),
avail = list(alt0=1, alt1=1),
choiceVar = Choice,
V = V,
rows = (CBC_type=="RD"),
componentName = componentName2
)
### Compute within-class choice probabilities using MNL model
P.componentName[[componentName2]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P.componentName[[componentName2]] = apollo_panelProd(P.componentName[[componentName2]], apollo_inputs, functionality)
### Combined model
P_within[[s]] = apollo_combineModels(P.componentName, apollo_inputs, functionality, asList=FALSE)
}
### Compute latent class model probabilities
lc_settings = list(inClassProb = P_within, classProb=pi_values)
P[["LC_model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
# ----------------------------------------------------------------- #
#---- Combine model components
# ----------------------------------------------------------------- #
### Combined model
P = apollo_combineModels(P, apollo_inputs, functionality)
# ### Take product across observation for same individual
# P = apollo_panelProd(P, apollo_inputs, functionality)
### Average across inter-individual draws
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
estimate_settings = list(hessianRoutine = "maxLik",
estimationRoutine = "bfgs",
maxIterations = 300,
printLevel = 3)
model = apollo_estimate(apollo_beta,
apollo_fixed,
apollo_probabilities,
apollo_inputs,
estimate_settings = estimate_settings)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
saveOutput_settings = list(
printClassical = TRUE,
printPVal = TRUE,
printDiagnostics = TRUE,
printCovar = TRUE,
printCorr = TRUE,
printOutliers = 100,
printChange = TRUE,
saveEst = TRUE,
saveCov = TRUE,
saveCorr = TRUE,
saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)