Given the high number of possible classes, I am using the methodology proposed by Hole, 2011. This gave me plausible results.
From this model, I wanted to test the idea that attendance might vary across choice tasks. To get started, I wan to test if the probability of ANA for the first choice task is different from the probability of the other choice tasks (overall 6 choice tasks)
I have written the following code. It runs, but I get the following warning:
"The class allocation probabilities for model component "model" are calculated at the observation level in 'apollo_lcPars', but are used in
'apollo_probabilities' to multiply within class probabilities that are at the individual level. Apollo will average the class allocation
probabilities across observations for the same individual level before using them to multiply the within-class probabilities. If your class
allocation probabilities are constant across choice situations for the same individual, then this is of no concern. If your class allocation
probabilities however vary across choice tasks, then you should change your model specification in 'apollo_probabilities' to only call
'apollo_panelProd' after calling 'apollo_lc'. "
Since it is the case, that probabilities would vary across choice tasks, I tried to follow the advice, but I am unable to find the right way to write it.
I have tried to insert the two lines
Code: Select all
lc_settings = list(inClassProb = P, classProb=pi_values)
P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
before the class loop, but it throws an error.
Here is the code I have written so far.
Code: Select all
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
library(ggplot2)
library(dplyr)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "ANA_MXL_pref_space_HiProbs_Diff1",
modelDescr = "ANA + MXL HiProbs stability of probabilities? ",
indivID = "interview__key",
nCores = 2,
outputDirectory = "02_scripts/_pourPapier/output/",
mixing = TRUE,
panelData = TRUE
)
# # load data
load(file = "forForum.RData")
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(mu_asc = -2,
mu_work = -2,
sigma_work = 1,
mu_fodd = 1.6,
sigma_fodd = 1,
mu_legu = 4,
sigma_legu = 1,
mu_engr = -1,
sigma_engr = 1,
mu_subs = 0.3,
sigma_subs = 0.2,
delta_anafodd = 1,
delta_analegu = 1,
delta_anaengr = 0.1,
delta_anasubs = 0.1 ,
gamma_fodd = 0.1 ,
gamma_legu = 0.1 ,
gamma_engr = 0.1 ,
gamma_subs = 0.1
)
apollo_fixed = c()
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType="sobolOwenFaureTezuka",
interNDraws= 150,
interUnifDraws=c(),
interNormDraws=c("draws_work", "draws_fodd","draws_legu","draws_engr", "draws_subs"),
intraDrawsType="mlhs",
intraNDraws=0,
intraUnifDraws=c(),
intraNormDraws=c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_work"]] = (mu_work + sigma_work * draws_work)
randcoeff[["b_work0"]] = 0* (mu_work + sigma_work * draws_work)
randcoeff[["b_fodd"]] = (mu_fodd + sigma_fodd * draws_fodd)
randcoeff[["b_fodd0"]] = 0* (mu_fodd + sigma_fodd * draws_fodd)
randcoeff[["b_legu"]] = (mu_legu + sigma_legu * draws_legu)
randcoeff[["b_legu0"]] = 0* (mu_fodd + sigma_fodd * draws_fodd)
randcoeff[["b_engr"]] = (mu_engr + sigma_engr * draws_engr)
randcoeff[["b_engr0"]] = 0* (mu_fodd + sigma_fodd * draws_fodd)
randcoeff[["b_subs"]] = (mu_subs + sigma_subs * draws_subs)
randcoeff[["b_subs0"]] = 0* (mu_fodd + sigma_fodd * draws_fodd)
return(randcoeff)
}
#### DEFINE LATENT CLASS COMPONENTS ####
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["beta_work"]] = list(b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 , b_work0 )
lcpars[["beta_fodd"]] = list(b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 , b_fodd0 )
lcpars[["beta_legu"]] = list(b_legu , b_legu , b_legu , b_legu , b_legu0 , b_legu0 , b_legu0 , b_legu0 , b_legu , b_legu , b_legu , b_legu , b_legu0 , b_legu0 , b_legu0 , b_legu0 , b_legu , b_legu , b_legu , b_legu , b_legu0 , b_legu0 , b_legu0 , b_legu0 , b_legu , b_legu , b_legu , b_legu , b_legu0 , b_legu0 , b_legu0 , b_legu0 )
lcpars[["beta_engr"]] = list(b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 , b_engr , b_engr , b_engr0 , b_engr0 )
lcpars[["beta_subs"]] = list(b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 , b_subs , b_subs0 )
# do not take into account the last 16 classes
lcpars[["beta_work"]] = lcpars[["beta_work"]][-c(17:32)]
lcpars[["beta_fodd"]] = lcpars[["beta_fodd"]][-c(17:32)]
lcpars[["beta_legu"]] = lcpars[["beta_legu"]][-c(17:32)]
lcpars[["beta_engr"]] = lcpars[["beta_engr"]][-c(17:32)]
lcpars[["beta_subs"]] = lcpars[["beta_subs"]][-c(17:32)]
### Generic settings for class allocation models
classAlloc_settings = list(
classes = c(class_a = 1, class_b = 2, class_c = 3, class_d = 4,
class_e = 5, class_f = 6, class_g = 7, class_h = 8,
class_i = 9, class_j = 10, class_k = 11, class_l = 12,
class_m = 13, class_n = 14, class_o = 15, class_p = 16
),
avail = 1
)
# intermediate deltas , if flag (CTASK >1) then we add a shifter
delta_fodd = delta_anafodd + gamma_fodd * flag
delta_legu = delta_analegu + gamma_legu * flag
delta_engr = delta_anaengr + gamma_engr * flag
delta_subs = delta_anasubs + gamma_subs * flag
# Calculation of probabilities based on the idinvidual deltas
V = list()
V[[1]] = 0
V[[2]] = delta_subs
V[[3]] = delta_engr
V[[4]] = delta_engr + delta_subs
V[[5]] = delta_legu
V[[6]] = delta_legu + delta_subs
V[[7]] = delta_legu + delta_engr
V[[8]] = delta_legu + delta_engr + delta_subs
V[[9]] = delta_fodd
V[[10]] = delta_fodd + delta_subs
V[[11]] = delta_fodd + delta_engr
V[[12]] = delta_fodd + delta_engr + delta_subs
V[[13]] = delta_fodd + delta_legu
V[[14]] = delta_fodd + delta_legu + delta_subs
V[[15]] = delta_fodd + delta_legu + delta_engr
V[[16]] = delta_fodd + delta_legu + delta_engr + delta_subs
classAlloc_settings$utilities = V
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
return(lcpars)
}
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("1"=1, "2"=2, "SQ"=3),
avail = list("1"=1, "2"=1, "SQ"=1),
choiceVar = CHOICE)
### Loop over classes
for(s in 1:16){
### Compute class-specific utilities
V=list()
V[["SQ"]] = mu_asc + beta_work[[s]]*WORKSQ/100 + beta_fodd[[s]]*FODDSQ/100 + beta_legu[[s]]*LEGUSQ/100 + beta_engr[[s]]*ENGRSQ/100 + beta_subs[[s]]*SUBSSQ/100
V[["1"]] = beta_work[[s]]*WORK1/100 + beta_fodd[[s]]*FODD1/100 + beta_legu[[s]]*LEGU1/100 + beta_engr[[s]]*ENGR1/100 + beta_subs[[s]]*SUBS1/100
V[["2"]] = beta_work[[s]]*WORK2/100 + beta_fodd[[s]]*FODD2/100 + beta_legu[[s]]*LEGU2/100 + beta_engr[[s]]*ENGR2/100 + beta_subs[[s]]*SUBS2/100
mnl_settings$utilities = V
mnl_settings$componentName = paste0("Class_",s)
### Compute within-class choice probabilities using MNL model
P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs ,functionality)
# ### Average across inter-individual draws within classes
P[[paste0("Class_",s)]] = apollo_avgInterDraws(P[[paste0("Class_",s)]], apollo_inputs, functionality)
}
### 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 AND OUTPUT ####
# ################################################################# #
model = apollo_estimate(apollo_beta,
apollo_fixed, apollo_probabilities, apollo_inputs)
apollo_modelOutput(model)