ApolloChoiceModelling forum

Forum for users of the Apollo choice modelling software
http://www.apollochoicemodelling.com/forum/

Posterior estimation latent class model with continuous random heterogeneity

http://www.apollochoicemodelling.com/forum/viewtopic.php?t=470

Page 1 of 1

Posterior estimation latent class model with continuous random heterogeneity

Posted: 12 Jun 2022, 04:13
by fayyazkhan
Hi Stephane and David,

I have estimated a latent class mixed logit with normal distribution for alternative specific constants (two classes, with only constant in the class allocation model). I would like to compute posterior values for different covariates (to profile individuals by figuring out which individuals belong to which class). Because of the random parameter, I cannot use apollo_lcConditionals (as this can only be used for latent class models without continuous random heterogeneity). I came across the below question on this forum, where Stephane has written a code for estimating posterior expected values of random coefficients.

http://www.apollochoicemodelling.com/fo ... ?f=16&t=50

I'm trying to implement the code in the above question for covariates in the class allocation model (e.g., gender), how should I change the last part of the above code. Many thanks.

Best regards,
Fayyaz



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()
}

Re: Posterior estimation latent class model with continuous random heterogeneity

Posted: 27 Jun 2022, 11:41
by stephanehess
Hi

apologies for the slow reply.

Could you share your model code, please?

Thanks

Stephane

Re: Posterior estimation latent class model with continuous random heterogeneity

Posted: 27 Jun 2022, 23:20
by fayyazkhan
Hi Stephane,

Thanks for your reply. Please see the below code.

Best,
Fayyaz

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  ="EC_LC_2_v1_0",
  modelDescr ="EC Latent class model with two classes no covariates",
  mixing     = TRUE,
  indivID    ="id",
  nCores = 3
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database = read.csv("data_2.csv",header=TRUE) # reads csv files


# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(mu_POSTIE_1             = 0.58729,
              mu_POSTIE_2             = 1.01714,
              POSTIE_nosafe_1         =-1.45606,
              POSTIE_nosafe_2         =-0.94492,
              POSTIE_regional_1       = 0.34415,
              POSTIE_regional_2       = 1.59364,
              mu_DRONE_1              =-0.48321,
              mu_DRONE_2              = 1.59027,
              DRONE_value100_1        = 0.13095,
              DRONE_value100_2        =-0.05695,
              DRONE_nosafe_1          =-2.37235,
              DRONE_nosafe_2          =-1.20607,
              LOCKER                  = 0.00000,
              LOCKER_value100_1       = 0.11362,
              LOCKER_value100_2       = 1.36205,
              SPEED5bd                = 0.00000,
              SPEED3bd_1              = 0.40906,
              SPEED3bd_2              = 0.26957,
              SPEED2bd_1              = 0.71360,
              SPEED2bd_2              = 0.39687,
              SPEEDnextbd_1           = 0.84788,
              SPEEDnextbd_2           = 0.93944,
              SPEEDsamed_1            = 1.80863,
              SPEEDsamed_2            = 0.91754,
              SPEED2hr_1              = 2.02390,
              SPEED2hr_2              = 0.82560,
              METHODfd                = 0.00000,
              METHODsafe_1            = 0.31718,
              METHODsafe_2            = 0.01002,
              METHODsig_1             = 0.15320,
              METHODsig_2             = 0.50737,
              METHODsig_nosafe_1      = 0.98429,
              METHODsig_nosafe_2      = 0.03496,
              METHODsig_lockdown_1    =-0.44788,
              METHODsig_lockdown_2    =-1.19845,
              WINDOWno                = 0.00000,
              WINDOW2hr_1             =-0.13254,
              WINDOW2hr_2             = 0.02532,
              WINDOWday_nosafe_1      = 0.16347,
              WINDOWday_nosafe_2      = 0.61101,
              WINDOW1hr_1             = 0.08986,
              WINDOW1hr_2             = 0.18075,
              WINDOW30min_1           = 0.04719,
              WINDOW30min_2           =-0.14218,
              WINDOWeven_1            =-0.19406,
              WINDOWeven_2            =-0.11872,
              COST_1                  =-0.22698,
              COST_2                  =-1.20887,
              COST_value50_1          = 0.07659,
              COST_value50_2          = 0.46906,
              COST_value100_1         = 0.20893,
              COST_value100_2         = 0.73285,
              Income_elas_1           =-0.19298,
              Income_elas_2           =-0.11279,
              COST_incomeNR_1         =-0.17783,
              COST_incomeNR_2         =-0.87834,
              delta_1                 = 0.46851,
              delta_2                 = 0.00000,
              sig_POSTIE_1            = 2.28863,
              sig_POSTIE_2            = 1.49950,
              sig_DRONE_1             =-3.04564,
              sig_DRONE_2             =-2.00225,
              ORDER1_1                = 0.180242,
              ORDER1_2                = 0.075402,
              ORDER2_1                = 0.094992,
              ORDER2_2                = 0.047412,
              ORDER3_1                = 0.000000,
              ORDER3_2                = 0.000000,
              LAMBDA_regional_1       = 0.947236,
              LAMBDA_regional_2       = 1.135653
)

### 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("LOCKER","SPEED5bd","METHODfd","WINDOWno","ORDER3_1","ORDER3_2",
                 "delta_2")

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 500,
  interUnifDraws = c(),
  interNormDraws = c("postie","drone"),
  intraDrawsType = "halton",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  
  randcoeff[["POSTIE_1"]]        =  mu_POSTIE_1  +  sig_POSTIE_1  *   postie 
  randcoeff[["POSTIE_2"]]        =  mu_POSTIE_2  +  sig_POSTIE_2  *   postie 
  randcoeff[["DRONE_1"]]         =  mu_DRONE_1   +  sig_DRONE_1   *   drone 
  randcoeff[["DRONE_2"]]         =  mu_DRONE_2   +  sig_DRONE_2   *   drone 
  
  
  
  return(randcoeff)
}
# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS                              ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
  lcpars = list()
  lcpars[["POSTIE"]]               = list(POSTIE_1            ,    POSTIE_2)
  lcpars[["POSTIE_nosafe"]]        = list(POSTIE_nosafe_1     ,    POSTIE_nosafe_2)
  lcpars[["POSTIE_regional"]]      = list(POSTIE_regional_1   ,    POSTIE_regional_2)
  lcpars[["DRONE"]]                = list(DRONE_1             ,    DRONE_2)
  lcpars[["DRONE_value100"]]       = list(DRONE_value100_1    ,    DRONE_value100_2)
  lcpars[["DRONE_nosafe"]]         = list(DRONE_nosafe_1      ,    DRONE_nosafe_2)
  lcpars[["LOCKER_value100"]]      = list(LOCKER_value100_1   ,    LOCKER_value100_2)
  lcpars[["SPEED3bd"]]             = list(SPEED3bd_1          ,    SPEED3bd_2)
  lcpars[["SPEED2bd"]]             = list(SPEED2bd_1          ,    SPEED2bd_2)
  lcpars[["SPEEDnextbd"]]          = list(SPEEDnextbd_1       ,    SPEEDnextbd_2)
  lcpars[["SPEEDsamed"]]           = list(SPEEDsamed_1        ,    SPEEDsamed_2)
  lcpars[["SPEED2hr"]]             = list(SPEED2hr_1          ,    SPEED2hr_2)
  lcpars[["METHODsafe"]]           = list(METHODsafe_1        ,    METHODsafe_2)
  lcpars[["METHODsig"]]            = list(METHODsig_1         ,    METHODsig_2)
  lcpars[["METHODsig_nosafe"]]     = list(METHODsig_nosafe_1  ,    METHODsig_nosafe_2)
  lcpars[["METHODsig_lockdown"]]   = list(METHODsig_lockdown_1,    METHODsig_lockdown_2)  
  lcpars[["WINDOW2hr"]]            = list(WINDOW2hr_1         ,    WINDOW2hr_2)
  lcpars[["WINDOWday_nosafe"]]     = list(WINDOWday_nosafe_1  ,    WINDOWday_nosafe_2)
  lcpars[["WINDOW1hr"]]            = list(WINDOW1hr_1         ,    WINDOW1hr_2)
  lcpars[["WINDOW30min"]]          = list(WINDOW30min_1       ,    WINDOW30min_2)
  lcpars[["WINDOWeven"]]           = list(WINDOWeven_1        ,    WINDOWeven_2)
  lcpars[["COST"]]                 = list(COST_1              ,    COST_2)
  lcpars[["COST_value50"]]         = list(COST_value50_1      ,    COST_value50_2)
  lcpars[["COST_value100"]]        = list(COST_value100_1     ,    COST_value100_2)
  lcpars[["Income_elas"]]          = list(Income_elas_1       ,    Income_elas_2)
  lcpars[["COST_incomeNR"]]        = list(COST_incomeNR_1     ,    COST_incomeNR_2)
  lcpars[["ORDER1"]]               = list(ORDER1_1            ,        ORDER1_2         )
  lcpars[["ORDER2"]]               = list(ORDER2_1            ,        ORDER2_2         )
  lcpars[["ORDER3"]]               = list(ORDER3_1            ,        ORDER3_2         )
  lcpars[["LAMBDA_regional"]]      = list(LAMBDA_regional_1   ,        LAMBDA_regional_2)
  
  
  
  
  V=list()
  V[["class_1"]] = delta_1 
  V[["class_2"]] = delta_2
  
  classAlloc_settings = list(
    classes      = c(class_1=1, class_2=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
  mnl_settings = list(
    alternatives  = c(postie=1, drone=2, locker=3), 
    avail         = 1, 
    choiceVar     = choice
  )
  
  
  ### Loop over classes
  for(s in 1:2){
    
    ### Compute class-specific utilities
    V = list()
    V[['postie']]  = LAMBDA_regional[[s]]^(metro==0)  * ( ORDER1[[s]] * (porder==1) + ORDER2[[s]] * (porder==2) + ORDER3[[s]] * (porder==3) + POSTIE[[s]] + POSTIE_nosafe[[s]] * (safepostie==0) + POSTIE_regional[[s]] * (metro==0)
                                                          + SPEEDnextbd[[s]] * (pspeed==4) + SPEED2bd[[s]] * (pspeed==5)+ SPEED3bd[[s]] * (pspeed==6) + SPEED5bd * (pspeed==1) 
                                                          + SPEED2hr[[s]] * (pspeed==2) + SPEEDsamed[[s]] * (pspeed==3) 
                                                          + METHODsafe[[s]] * (pmethod==1) + METHODsig[[s]] * (pmethod==2) + METHODfd * (pmethod==3) + METHODsig[[s]] * (pmethod==2) + METHODsig_nosafe[[s]] * (pmethod==2) * (safepostie==0) + METHODsig_lockdown[[s]] * (pmethod==2) * (pilot==0)
                                                          + WINDOWno * (pwindow==1) + WINDOW30min[[s]] * (pwindow==2) + WINDOW1hr[[s]] * (pwindow==3) + WINDOW2hr[[s]] * (pwindow==4) + WINDOWday_nosafe[[s]] * (pwindow==2) * (safepostie==0) + WINDOWday_nosafe[[s]] * (pwindow==3) * (safepostie==0) + WINDOWday_nosafe[[s]] * (pwindow==4) * (safepostie==0) + WINDOWeven[[s]] * (pwindow==5)
                                                          + (householdincome<12)  * (COST[[s]]          *  pcost + COST_value50[[s]]          * pcost * pvalue50 + COST_value100[[s]]          * pcost * pvalue100 ) * (con_hhincome/mincome)^Income_elas[[s]] 
                                                          + (householdincome==12) * COST_incomeNR[[s]] *  pcost ) 
    V[['drone']]   = LAMBDA_regional[[s]]^(metro==0)   * ( ORDER1[[s]] * (dorder==1) + ORDER2[[s]] * (dorder==2) + ORDER3[[s]] * (dorder==3) + DRONE[[s]] + DRONE_nosafe[[s]] * (safedrone==0) + DRONE_value100[[s]] * pvalue100 
                                                           + SPEEDnextbd[[s]] * (dspeed==4) + SPEED2bd[[s]] * (dspeed==5) + SPEED3bd[[s]] * (dspeed==6) + SPEED5bd * (dspeed==1) 
                                                           + SPEED2hr[[s]] * (dspeed==2) + SPEEDsamed[[s]] * (dspeed==3)
                                                           + METHODsafe[[s]] * (dmethod==1) + METHODsig[[s]] * (dmethod==2) + METHODfd * (dmethod==3) 
                                                           + WINDOWno * (dwindow==1) + WINDOW30min[[s]] * (dwindow==2) + WINDOW1hr[[s]] * (dwindow==3) + WINDOW2hr[[s]] * (dwindow==4) + WINDOWday_nosafe[[s]] * (dwindow==2) * (safedrone==0) + WINDOWday_nosafe[[s]] * (dwindow==3) * (safedrone==0) + WINDOWday_nosafe[[s]] * (dwindow==4) * (safedrone==0) + WINDOWeven[[s]] * (dwindow==5)
                                                           + (householdincome<12)  * (COST[[s]]          *  dcost + COST_value50[[s]]          * dcost * pvalue50 + COST_value100[[s]]          * dcost * pvalue100 ) * (con_hhincome/mincome)^Income_elas[[s]] 
                                                           + (householdincome==12) * COST_incomeNR[[s]] * dcost )                                   
    V[['locker']]  = LAMBDA_regional[[s]]^(metro==0)   * ( ORDER1[[s]] * (lorder==1) + ORDER2[[s]] * (lorder==2) + ORDER3[[s]] * (lorder==3) + LOCKER + LOCKER_value100[[s]] * pvalue100
                                                           + SPEEDnextbd[[s]] * (lspeed==4) + SPEED2bd[[s]] * (lspeed==5) + SPEED3bd[[s]] * (lspeed==6) + SPEED5bd * (lspeed==1) + SPEED2hr[[s]] * (lspeed==2) +SPEEDsamed[[s]] * (lspeed==3) 
                                                           + (householdincome<12)  * (COST[[s]]          *  lcost + COST_value50[[s]]          * lcost * pvalue50 + COST_value100[[s]]          * lcost * pvalue100 ) * (con_hhincome/mincome)^Income_elas[[s]] 
                                                           + (householdincome==12) * COST_incomeNR[[s]] * lcost ) 
    
    
    
    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)
}


# ################################################################# #
#### EM ESTIMATION FOR COVARIANCE MATRIX                         ####
# ################################################################# #

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model,modelOutput_settings=list(printClassical=TRUE,printT1=TRUE))

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model,saveOutput_settings=list(printClassical=TRUE,printT1=TRUE))

Re: Posterior estimation latent class model with continuous random heterogeneity

Posted: 10 Oct 2022, 17:04
by stephanehess
Fayyaz

apologies for the slow reply. If you can share the data with me, then I'll try to prepare the code for you

Stephane
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited