I am trying to run some post-estimation analytics, where I would like to look into my classes in the latent class model to see what characterizes them. In this forum I found a very useful code that you wrote for a latent class mixed logit estimation, but some of the code is not applicable for me since I do not run the mixed logit part.
Below I attach my code:
Code: Select all
### Clear memory
rm(list = ls())
### Load libraries
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "LC_6CLASS_ASC",
modelDescr = "Latent Class model with SIX classes, removing socios and changing the specification of the ASCs.",
indivID = "ID",
nCores = 6,
outputDirectory = "output",
panelData = TRUE
)
#Använd olika startvärden samt olika algoritmer för att se vart jag kan hitta ett så nära globalt maximum som möjligt.
### Set working directory
setwd("C:/Users/gaso0009/OneDrive - Umeå universitet/RESEARCH/PAPER 1/R CODES")
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading choice data from excel
SOL_KOMBO <- read_excel("C:/Users/gaso0009/OneDrive - Umeå universitet/RESEARCH/PAPER 1/SOL_KOMBO.xlsx")
database <- SOL_KOMBO
### Create a choice variable, consisting of choice1, choice2, and choice3
database$choice <- ifelse(database$choice1 == 1, 1,
ifelse(database$choice2 == 1, 2,
ifelse(database$choice3 == 1, 3, NA)))
### Rescaling investment cost to thousands
database$invcost1 <- database$invcost1 / 1000
database$invcost2 <- database$invcost2 / 1000
database$invcost3 <- database$invcost3 / 1000
### Split up the gender parameter into female and male
database$female <- ifelse(database$gender == 2, 1, 0)
database$male <- ifelse(database$gender == 1, 1, 0)
database$nonbinary <- ifelse(database$gender == 3, 1, 0)
database$gender_other <- ifelse(database$gender == 4, 1, 0)
database$nogender <- ifelse(database$gender == 5, 1, 0)
###Mean income
#database$mean_income = mean(database$income_continuous)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(b_strg_no_a = -0.1,
b_strg_no_b = 0,
b_strg_no_c = 0,
b_strg_no_d = 0,
b_strg_no_e = 0,
b_strg_no_f = 0,
b_strg_1_a = -0.1,
b_strg_1_b = 0,
b_strg_1_c = 0,
b_strg_1_d = 0,
b_strg_1_e = 0,
b_strg_1_f = 0,
b_strg_10_a = -0.1,
b_strg_10_b = 0,
b_strg_10_c = 0,
b_strg_10_d = 0,
b_strg_10_e = 0,
b_strg_10_f = 0,
b_strg_24_a = -0.1,
b_strg_24_b = 0,
b_strg_24_c = 0,
b_strg_24_d = 0,
b_strg_24_e = 0,
b_strg_24_f = 0,
b_peer_low_a = -0.1,
b_peer_low_b = 0,
b_peer_low_c = 0,
b_peer_low_d = 0,
b_peer_low_e = 0,
b_peer_low_f = 0,
b_peer_med_a = -0.1,
b_peer_med_b = 0,
b_peer_med_c = 0,
b_peer_med_d = 0,
b_peer_med_e = 0,
b_peer_med_f = 0,
b_peer_high_a = -0.1,
b_peer_high_b = 0,
b_peer_high_c = 0,
b_peer_high_d = 0,
b_peer_high_e = 0,
b_peer_high_f = 0,
b_es_yes_a = -0.1,
b_es_yes_b = 0,
b_es_yes_c = 0,
b_es_yes_d = 0,
b_es_yes_e = 0,
b_es_yes_f = 0,
b_es_no_a = -0.1,
b_es_no_b = 0,
b_es_no_c = 0,
b_es_no_d = 0,
b_es_no_e = 0,
b_es_no_f = 0,
b_invcost_a = -0.1,
b_invcost_b = 0,
b_invcost_c = 0,
b_invcost_d = 0,
b_invcost_e = 0,
b_invcost_f = 0,
b_sq = 0, #Bara ASC för opt out?? Har opt out ett egenvärde?
delta_a = 0,
delta_b = 0,
delta_c = 0,
delta_d = 0,
delta_e = 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("b_strg_no_a", "b_strg_no_b", "b_strg_no_c", "b_strg_no_d", "b_strg_no_e", "b_strg_no_f", "b_peer_low_a", "b_peer_low_b", "b_peer_low_c", "b_peer_low_d", "b_peer_low_e", "b_peer_low_f", "b_es_no_a", "b_es_no_b", "b_es_no_c", "b_es_no_d", "b_es_no_e", "b_es_no_f")
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_strg_no"]] = list(b_strg_no_a, b_strg_no_b, b_strg_no_c, b_strg_no_d, b_strg_no_e, b_strg_no_f)
lcpars[["b_strg_1"]] = list(b_strg_1_a, b_strg_1_b, b_strg_1_c, b_strg_1_d, b_strg_1_e, b_strg_1_f)
lcpars[["b_strg_10"]] = list(b_strg_10_a, b_strg_10_b, b_strg_10_c, b_strg_10_d, b_strg_10_e, b_strg_10_f)
lcpars[["b_strg_24"]] = list(b_strg_24_a, b_strg_24_b, b_strg_24_c, b_strg_24_d, b_strg_24_e, b_strg_24_f)
lcpars[["b_peer_low"]] = list(b_peer_low_a, b_peer_low_b, b_peer_low_c, b_peer_low_d, b_peer_low_e, b_peer_low_f)
lcpars[["b_peer_med"]] = list(b_peer_med_a, b_peer_med_b, b_peer_med_c, b_peer_med_d, b_peer_med_e, b_peer_med_f)
lcpars[["b_peer_high"]] = list(b_peer_high_a, b_peer_high_b, b_peer_high_c, b_peer_high_d, b_peer_high_e, b_peer_high_f)
lcpars[["b_es_yes"]] = list(b_es_yes_a, b_es_yes_b, b_es_yes_c, b_es_yes_d, b_es_yes_e, b_es_yes_f)
lcpars[["b_es_no"]] = list(b_es_no_a, b_es_no_b, b_es_no_c, b_es_no_d, b_es_no_e, b_es_no_f)
lcpars[["b_invcost"]] = list(b_invcost_a, b_invcost_b, b_invcost_c, b_invcost_d, b_invcost_e, b_invcost_f)
V=list()
V[["class_a"]] = delta_a
V[["class_b"]] = delta_b
V[["class_c"]] = delta_c
V[["class_d"]] = delta_d
V[["class_e"]] = delta_e
V[["class_f"]] = 0
#Testa ta bort dessa
classAlloc_settings = list(
alternatives = c(class_a=1, class_b=2, class_c=3, class_d=4, class_e=5, class_f=6),
V = V
)
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
return(lcpars)
}
### Making sure that alternatives matches the choices
apollo_alternatives <- c(1, 2, 3) # Representing your alternatives
# ################################################################# #
#### 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 mnl_settings for apollo_mnl
mnl_settings = list(
alternatives = c(alt1 = 1, alt2 = 2, alt3 = 3), # Define each alternative
avail = 1, # Assume all alternatives are available
choiceVar = choice, # Variable for actual choice
panelData = TRUE
)
### Loop over classes
for(s in 1:6){
### List of utilities
V = list()
V[["alt1"]] = b_strg_no[[s]] * strg_no1 + b_strg_1[[s]] * strg_11 + b_strg_10[[s]] * strg_101 + b_strg_24[[s]] * strg_241 + b_peer_low[[s]] * peer_low1 + b_peer_med[[s]] * peer_med1 + b_peer_high[[s]] * peer_high1 + b_es_yes[[s]] * es_yes1 + b_es_no[[s]] * es_no1 + b_invcost[[s]] * invcost1
V[["alt2"]] = b_strg_no[[s]] * strg_no2 + b_strg_1[[s]] * strg_12 + b_strg_10[[s]] * strg_102 + b_strg_24[[s]] * strg_242 + b_peer_low[[s]] * peer_low2 + b_peer_med[[s]] * peer_med2 + b_peer_high[[s]] * peer_high2 + b_es_yes[[s]] * es_yes2 + b_es_no[[s]] * es_no2 + b_invcost[[s]] * invcost2
V[["alt3"]] = b_sq
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)
}
### 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)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
### Printing p-values
library(dplyr)
# Extract estimates and standard errors
estimates <- model$estimate
robust_se <- model$robse
# Calculate p-values
p_values <- 2 * (1 - pnorm(abs(model$estimate / model$robse)))
# Combine into a table with Estimates, Standard Errors, T-Ratios, and P-Values
results_table <- data.frame(
Parameter = names(estimates),
Estimate = estimates,
Rob.s.e. = robust_se,
T_Ratio = estimates / robust_se,
P_Value = p_values
)
# Print the table
print(results_table)
# ----------------------------------------------------------------- #
#---- LR TESTS & BAS TEST AGAINST EARLIER MODELS
# ----------------------------------------------------------------- #
### Comparing the models
#Likelihood ratio test
apollo_lrTest("MAIN", model)
apollo_lrTest("MAIN_ASC", model)
apollo_lrTest("MMNL_FULL_SPEC_NORMAL", model)
apollo_lrTest("MMNL_FULL_SPEC_ASC", model)
apollo_lrTest("MMNL_STORAGE_PEER_ES", model)
apollo_lrTest("MMNL_STORAGE_NORMAL", model)
apollo_lrTest("LC_2CLASS_2", model)
apollo_lrTest("LC_3CLASS_2", model)
apollo_lrTest("LC_4CLASS_2", model)
apollo_lrTest("LC_5CLASS_2", model)
apollo_lrTest("MMNL_STORAGE_INTER", model)
apollo_lrTest("MMNL_FULL_SPEC_INTER", model)
#Ben-Akiva & Swait test
apollo_basTest("MAIN", model)
apollo_basTest("MAIN_ASC", model)
apollo_basTest("MMNL_FULL_SPEC_NORMAL", model)
apollo_basTest("MMNL_STORAGE_PEER_ES", model)
apollo_basTest("MMNL_STORAGE_NORMAL", model)
apollo_basTest("LC_2CLASS_2", model)
apollo_basTest("LC_3CLASS_2", model)
apollo_basTest("LC_4CLASS_2", model)
apollo_basTest("LC_5CLASS_2", model)
apollo_basTest("LC_6CLASS_2", model)
apollo_basTest("MMNL_STORAGE_INTER", model)
apollo_basTest("MMNL_STORAGE_INTER_ASC", model)
apollo_basTest("MMNL_FULL_SPEC_INTER", model)
apollo_basTest("MMNL_FULL_SPEC_INTER_ASC", model)
# ----------------------------------------------------------------- #
#---- CALCULATE WTP ----
# ----------------------------------------------------------------- #
apollo_deltaMethod(model,
deltaMethod_settings = list(
expression=c(WTP_strg1_a = "(b_strg_1_a)/b_invcost_a",
WTP_strg10_a = "(b_strg_10_a)/b_invcost_a",
WTP_strg24_a = "(b_strg_24_a)/b_invcost_a",
WTP_es_yes_a = "(b_es_yes_a)/b_invcost_a",
WTP_peer_med_a = "(b_peer_med_a)/b_invcost_a",
WTP_peer_high_a = "(b_peer_high_a)/b_invcost_a",
WTP_strg1_b = "(b_strg_1_b)/b_invcost_b",
WTP_strg10_b = "(b_strg_10_b)/b_invcost_b",
WTP_strg24_b = "(b_strg_24_b)/b_invcost_b",
WTP_es_yes_b = "(b_es_yes_b)/b_invcost_b",
WTP_peer_med_b = "(b_peer_med_b)/b_invcost_b",
WTP_peer_high_b = "(b_peer_high_b)/b_invcost_b",
WTP_strg1_c = "(b_strg_1_c)/b_invcost_c",
WTP_strg10_c = "(b_strg_10_c)/b_invcost_c",
WTP_strg24_c = "(b_strg_24_c)/b_invcost_c",
WTP_es_yes_c = "(b_es_yes_c)/b_invcost_c",
WTP_peer_med_c = "(b_peer_med_c)/b_invcost_c",
WTP_peer_high_c = "(b_peer_high_c)/b_invcost_c",
WTP_strg1_d = "(b_strg_1_d)/b_invcost_d",
WTP_strg10_d = "(b_strg_10_d)/b_invcost_d",
WTP_strg24_d = "(b_strg_24_d)/b_invcost_d",
WTP_es_yes_d = "(b_es_yes_d)/b_invcost_d",
WTP_peer_med_d = "(b_peer_med_d)/b_invcost_d",
WTP_peer_high_d = "(b_peer_high_d)/b_invcost_d",
WTP_strg1_e = "(b_strg_1_e)/b_invcost_e",
WTP_strg10_e = "(b_strg_10_e)/b_invcost_e",
WTP_strg24_e = "(b_strg_24_e)/b_invcost_e",
WTP_es_yes_e = "(b_es_yes_e)/b_invcost_e",
WTP_peer_med_e = "(b_peer_med_e)/b_invcost_e",
WTP_peer_high_e = "(b_peer_high_e)/b_invcost_e",
WTP_strg1_f = "(b_strg_1_f)/b_invcost_f",
WTP_strg10_f = "(b_strg_10_f)/b_invcost_f",
WTP_strg24_f = "(b_strg_24_f)/b_invcost_f",
WTP_es_yes_f = "(b_es_yes_f)/b_invcost_f",
WTP_peer_med_f = "(b_peer_med_f)/b_invcost_f",
WTP_peer_high_f = "(b_peer_high_f)/b_invcost_f")))
apollo_saveOutput(model)
# ################################################################# #
##### POST-PROCESSING ####
# ################################################################# #
### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
apollo_sink()
# ----------------------------------------------------------------- #
#---- POSTERIOR ANALYSIS ----
# ----------------------------------------------------------------- #
# ################################################################# #
#### CONDITIONALS FOR CLASS ALLOCATION ####
# ################################################################# #
apollo_beta = model$estimate
class_prob = "pi_values"
L = apollo_probabilities(apollo_beta, apollo_inputs, functionality="output")
apollo_attach(apollo_beta, apollo_inputs)
lcpars = apollo_lcPars(apollo_beta, apollo_inputs)
classes = length(lcpars[[class_prob]])
post_pi = vector(mode="list", length=classes)
for(s in 1:classes) post_pi[[s]] = lcpars[[class_prob]][[s]]*L[[s]]/L[["model"]]
lc_conditionals = matrix(unlist(post_pi), ncol = length(post_pi), byrow = FALSE)
classnames = paste("Class ",seq(1:classes),sep="")
lc_conditionals=cbind(unique(database[,apollo_control$indivID]),lc_conditionals)
colnames(lc_conditionals) = c("ID",classnames)
apollo_detach(apollo_beta, apollo_inputs)
# ################################################################# #
#### CONDITIONALS FOR CONTINUOUS PARAMS ####
# ################################################################# #
apollo_probabilities_class_specific=function(apollo_beta, apollo_inputs, functionality="estimate"){
### Function initialisation: do not change the following three commands
### 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
mnl_settings = list(
alternatives = c(alt1 = 1, alt2 = 2, alt3 = 3),
avail = 1,
choiceVar = choice,
panelData = TRUE
)
s=apollo_inputs$s
### Compute class-specific utilities
### List of utilities
V = list()
V[["alt1"]] = b_strg_no[[s]] * strg_no1 + b_strg_1[[s]] * strg_11 + b_strg_10[[s]] * strg_101 + b_strg_24[[s]] * strg_241 + b_peer_low[[s]] * peer_low1 + b_peer_med[[s]] * peer_med1 + b_peer_high[[s]] * peer_high1 + b_es_yes[[s]] * es_yes1 + b_es_no[[s]] * es_no1 + b_invcost[[s]] * invcost1
V[["alt2"]] = b_strg_no[[s]] * strg_no2 + b_strg_1[[s]] * strg_12 + b_strg_10[[s]] * strg_102 + b_strg_24[[s]] * strg_242 + b_peer_low[[s]] * peer_low2 + b_peer_med[[s]] * peer_med2 + b_peer_high[[s]] * peer_high2 + b_es_yes[[s]] * es_yes2 + b_es_no[[s]] * es_no2 + b_invcost[[s]] * invcost2
V[["alt3"]] = b_sq
mnl_settings$V = V
### Compute within-class choice probabilities using MNL model
P[["model"]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs ,functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
apollo_beta=model$estimate
P=list()
for(s in 1:6){
apollo_inputs$s=s
P[[s]]=apollo_probabilities_class_specific(apollo_beta, apollo_inputs, functionality="conditionals")
}
apollo_attach(apollo_beta, apollo_inputs)
obsPerIndiv <- as.vector(table(database[,apollo_control$indivID]))
conditionals=list()
class1_pars=c("b_strg_no_a",
"b_strg_1_a",
"b_strg_10_a",
"b_strg_24_a",
"b_peer_low_a",
"b_peer_med_a",
"b_peer_high_a",
"b_es_no_a",
"b_es_yes_a",
"b_invcost_b")
class2_pars=c("b_strg_no_b",
"b_strg_1_b",
"b_strg_10_b",
"b_strg_24_b",
"b_peer_low_b",
"b_peer_med_b",
"b_peer_high_b",
"b_es_no_b",
"b_es_yes_b",
"b_invcost_b")
class3_pars=c("b_strg_no_c",
"b_strg_1_c",
"b_strg_10_c",
"b_strg_24_c",
"b_peer_low_c",
"b_peer_med_c",
"b_peer_high_c",
"b_es_no_c",
"b_es_yes_c",
"b_invcost_c")
class4_pars=c("b_strg_no_d",
"b_strg_1_d",
"b_strg_10_d",
"b_strg_24_d",
"b_peer_low_d",
"b_peer_med_d",
"b_peer_high_d",
"b_es_no_d",
"b_es_yes_d",
"b_invcost_d")
class5_pars=c("b_strg_no_e",
"b_strg_1_e",
"b_strg_10_e",
"b_strg_24_e",
"b_peer_low_e",
"b_peer_med_e",
"b_peer_high_e",
"b_es_no_e",
"b_es_yes_e",
"b_invcost_e")
class6_pars=c("b_strg_no_f",
"b_strg_1_f",
"b_strg_10_f",
"b_strg_24_f",
"b_peer_low_f",
"b_peer_med_f",
"b_peer_high_f",
"b_es_no_f",
"b_es_yes_f",
"b_invcost_f")
Until here, the code works perfectly. I just do not know how to alter the code to fit my data and how to move on from here.
Here is the continuation of the code I found:
Code: Select all
for(j in 1:length(randcoeff)){
b=randcoeff[[j]]
b <- rowsum(b, group=database[,apollo_control$indivID])
b=b/obsPerIndiv
P_use=(names(randcoeff)[j]%in%class1_pars)*P[[1]]+(names(randcoeff)[j]%in%class2_pars)*P[[2]]+(names(randcoeff)[j]%in%class3_pars)*P[[3]]
bn=(rowSums(b*P_use))/rowSums(P_use)
bns=sqrt(rowSums(P_use*(b-bn)^2)/(rowSums(P_use)))
conditionals[[names(randcoeff)[j]]]=cbind(unique(database[,apollo_control$indivID]),bn,bns)
colnames(conditionals[[names(randcoeff)[j]]])=c("ID","post. mean","post. sd")
rownames(conditionals[[names(randcoeff)[j]]])=c()
}
apollo_detach(apollo_beta, apollo_inputs)How do I look more into this?
Best,
Gabrielle Söderberg