Posterior probability more than 100%?
Posted: 05 Feb 2024, 21:44
Hi,
I am running a model using the latent class model to identify preferences for choosing a medical device. When I conducted the posterior analysis, no errors or warnings showed. However, I found that the maximum posterior conditional probability is over 100%. I checked over and over again but still have no clue what happened. The following are my codes:
Here is the output of unconditional/conditional probabilities:
Is there any error in my codes to make the probability over 100%?
Thanks in advance for answering my questions.
Tim
I am running a model using the latent class model to identify preferences for choosing a medical device. When I conducted the posterior analysis, no errors or warnings showed. However, I found that the maximum posterior conditional probability is over 100%. I checked over and over again but still have no clue what happened. The following are my codes:
Code: Select all
# ################################################################# #
#### 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 = "LC model",
modelDescr = "Latent class model (n=2)",
indivID = "ID",
nCores = 2,
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("cgm_N182d.csv",header=TRUE)
### for data dictionary, use ?apollo_swissRouteChoiceData
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc1 = 0,
asc2 = 0,
#Manual/Auto
#b_auto_a = 0,
#b_auto_b = 0,
b_manual_a = 0,
b_manual_b = -0.1382,
#Alarm
#b_noalm_a = 0,
#b_noalm_b = 0,
b_ncstalm_a = 0,
b_ncstalm_b = 0.05536,
b_cstalm_a = 0,
b_cstalm_b = 0.19059,
#Info
#b_5m_nop_a = 0,
#b_5m_nop_b = 0,
b_1m_nop_a = 0,
b_1m_nop_b = 0.06571,
b_5m_p_a = 0,
b_5m_p_b = 0.01974,
b_1m_p_a = 0,
b_1m_p_b = -0.09620,
#Acurcy
#b_acurcy0_a = 0,
#b_acurcy0_b = 0,
b_acurcy10_a = 0,
b_acurcy10_b = -0.01855,
b_acurcy15_a = 0,
b_acurcy15_b = -0.28813,
##Adjustment
#b_cali0_a = 0,
#b_cali0_b = 0,
b_cali2_a = 0,
b_cali2_b = 0.02384,
b_cali4_a = 0,
b_cali4_b = -0.25955,
#Battery
#b_life1_a =0,
b_life8_a =0,
b_life8_b =0.04006,
b_life24_a =0,
b_life24_b =0.18016,
#b_cost = 0 # When treat as continuous
#b_cost50_a = 0,
#b_cost50_b = 0,
b_cost150_a = 0,
b_cost150_b = 0.01128,
b_cost300_a = 0,
b_cost300_b = 0.09332,
b_cost500_a = 0,
b_cost500_b = -0.20953,
#Latent class specification
delta_a = 0,
delta_b = 0.1,
gamma_user_a = 0,
gamma_user_b = 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("asc1","delta_a","gamma_user_a")
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #
apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()
lcpars[["b_auto"]] = list(-b_manual_a, -b_manual_b) # Effects coding constraint
lcpars[["b_manual"]] = list(b_manual_a, b_manual_b)
lcpars[["b_noalm"]] = list(-b_ncstalm_a - b_cstalm_a, -b_ncstalm_b - b_cstalm_b)
lcpars[["b_ncstalm"]] = list(b_ncstalm_a, b_ncstalm_b)
lcpars[["b_cstalm"]] = list(b_cstalm_a, b_cstalm_b)
lcpars[["b_5m_nop"]] = list(-b_5m_p_a - b_1m_nop_a - b_1m_p_a, - b_5m_p_b - b_1m_nop_b - b_1m_p_b)
lcpars[["b_1m_nop"]] = list(b_1m_nop_a, b_1m_nop_b)
lcpars[["b_5m_p"]] = list(b_5m_p_a, b_5m_p_b)
lcpars[["b_1m_p"]] = list(b_1m_p_a, b_1m_p_b)
lcpars[["b_acurcy0"]] = list(-b_acurcy10_a - b_acurcy15_a, -b_acurcy10_b - b_acurcy15_b)
lcpars[["b_acurcy10"]] = list(b_acurcy10_a, b_acurcy10_b)
lcpars[["b_acurcy15"]] = list(b_acurcy15_a, b_acurcy15_b)
lcpars[["b_cali0"]] = list(-b_cali2_a - b_cali4_a, - b_cali2_b - b_cali4_b)
lcpars[["b_cali2"]] = list(b_cali2_a, b_cali2_b)
lcpars[["b_cali4"]] = list(b_cali4_a, b_cali4_b)
lcpars[["b_life1"]] = list(- b_life8_a - b_life24_a, - b_life8_b - b_life24_b)
lcpars[["b_life8"]] = list(b_life8_a, b_life8_b)
lcpars[["b_life24"]] = list(b_life24_a, b_life24_b)
lcpars[["b_cost50"]] = list(- b_cost150_a - b_cost300_a - b_cost500_a, - b_cost150_b - b_cost300_b - b_cost500_b)
lcpars[["b_cost150"]] = list(b_cost150_a, b_cost150_b)
lcpars[["b_cost300"]] = list(b_cost300_a, b_cost300_b)
lcpars[["b_cost500"]] = list(b_cost500_a, b_cost500_b)
### Utilities of class allocation model
#Without covariate
# V=list()
# V[["class_a"]] = delta_a
# V[["class_b"]] = delta_b
#With covariate
V=list()
V[["class_a"]] = delta_a + gamma_user_a*Rvl_device #Rvl_device <- current user?
V[["class_b"]] = delta_b + gamma_user_b*Rvl_device
classAlloc_settings = list(
classes = c(class_a=1, class_b=2),
utilities = V
)
lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
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),
avail = list(alt1=1, alt2=1),
choiceVar = Choice
)
### Loop over classes
for(s in 1:2){
### Compute class-specific utilities
V=list()
V[["alt1"]] = (asc1
+ b_auto[[s]]*(effort_1==0) + b_manual[[s]]*(effort_1==1)
+ b_noalm[[s]]*(alarm_1==0) + b_ncstalm[[s]]*(alarm_1==1) + b_cstalm[[s]]*(alarm_1==2)
+ b_5m_nop[[s]]*(info_1==0) + b_1m_nop[[s]]*(info_1==1) + b_5m_p[[s]]*(info_1==2) + b_1m_p[[s]]*(info_1==3)
+ b_acurcy0[[s]]*(accuracy_1==0)+b_acurcy10[[s]]*(accuracy_1==10)+b_acurcy15[[s]]*(accuracy_1==15)
+ b_cali0[[s]]*(cali_1==0)+b_cali2[[s]]*(cali_1==2)+b_cali4[[s]]*(cali_1==4)
+ b_life1[[s]]*(senslife_1==1)+b_life8[[s]]*(senslife_1==8)+b_life24[[s]]*(senslife_1==24)
+ b_cost50[[s]]*(cost_1==50)+b_cost150[[s]]*(cost_1==150)+b_cost300[[s]]*(cost_1==300)+b_cost500[[s]]*(cost_1==500)
)
V[["alt2"]] = (asc2
+ b_auto[[s]]*(effort_2==0) + b_manual[[s]]*(effort_2==1)
+ b_noalm[[s]]*(alarm_2==0) + b_ncstalm[[s]]*(alarm_2==1) + b_cstalm[[s]]*(alarm_2==2)
+ b_5m_nop[[s]]*(info_2==0) + b_1m_nop[[s]]*(info_2==1) + b_5m_p[[s]]*(info_2==2) + b_1m_p[[s]]*(info_2==3)
+ b_acurcy0[[s]]*(accuracy_2==0)+b_acurcy10[[s]]*(accuracy_2==10)+b_acurcy15[[s]]*(accuracy_2==15)
+ b_cali0[[s]]*(cali_2==0)+b_cali2[[s]]*(cali_2==2)+b_cali4[[s]]*(cali_2==4)
+ b_life1[[s]]*(senslife_2==1)+b_life8[[s]]*(senslife_2==8)+b_life24[[s]]*(senslife_2==24)
+ b_cost50[[s]]*(cost_2==50)+b_cost150[[s]]*(cost_2==150)+b_cost300[[s]]*(cost_2==300)+b_cost500[[s]]*(cost_2==500)
)
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 ####
# ################################################################ #
### Optional starting values search
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ################################################################# #
##### POST-PROCESSING ####
# ################################################################# #
### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
apollo_sink()
# ----------------------------------------------------------------- #
#---- ANALYSIS OF CLASS ALLOCATION ----
# ----------------------------------------------------------------- #
### Unconditional class allocation probabilities
unconditionals = apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
summary(as.data.frame(unconditionals[["pi_values"]]))
### Conditional (posterior) class allocation probabilities
conditionals = apollo_conditionals(model,apollo_probabilities, apollo_inputs)
summary(conditionals)
Code: Select all
> unconditionals = apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
Unconditional distributions computed
> summary(as.data.frame(unconditionals[["pi_values"]]))
class_a class_b
Min. :0.2906 Min. :0.5400
1st Qu.:0.2906 1st Qu.:0.5400
Median :0.2906 Median :0.7094
Mean :0.3734 Mean :0.6266
3rd Qu.:0.4600 3rd Qu.:0.7094
Max. :0.4600 Max. :0.7094
>
> ### Conditional (posterior) class allocation probabilities
> conditionals = apollo_conditionals(model,apollo_probabilities, apollo_inputs)
Calculating conditionals...
> summary(conditionals)
ID X1 X2
Min. : 1.00 Min. :0.0000000 Min. :0.001094
1st Qu.: 46.25 1st Qu.:0.0000165 1st Qu.:0.043277
Median : 91.50 Median :0.0080597 Median :0.775052
Mean : 91.50 Mean :0.4246550 Mean :0.660761
3rd Qu.:136.75 3rd Qu.:0.9462279 3rd Qu.:1.000000
Max. :182.00 Max. :1.5789786 Max. :1.313540
Thanks in advance for answering my questions.
Tim