-
-
Save HyewonNamgung/c761f024f15290812440c81171f36d45 to your computer and use it in GitHub Desktop.
Latent Class_EMalgorithm with 3 classes.R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # ################################################################# # | |
| #### 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 = "EM_LC_with_covariates", | |
| modelDescr = "LC model with class allocation model on Swiss route choice data, EM algorithm", | |
| indivID = "ID", | |
| noValidation = TRUE, | |
| noDiagnostics = TRUE, | |
| outputDirectory = "output" | |
| ) | |
| # ################################################################# # | |
| #### LOAD DATA AND APPLY ANY TRANSFORMATIONS #### | |
| # ################################################################# # | |
| # variable | |
| database$job1 <-ifelse(database$job==1,1,0) #farmer-non | |
| database$job_no <-ifelse(database$job==6|database$job==12,1,0) #job-non | |
| database$young <- ifelse(database$age<3,1,0) | |
| database$middle <- ifelse(database$age==3|database$age==4|database$age==5 ,1,0) | |
| database$elderly <- ifelse(database$age>5,1,0) #elderly | |
| database$commuter <- ifelse(database$act1==2,1,0) # commuter | |
| database$iwami_commute <- ifelse(database$place1==1,1,0) # iwami inner commuter | |
| database$public_commute <- ifelse(database$to1==1|database$to1==2,1,0) # commute private car | |
| database$shop_no <- ifelse(database$act2==1,1,0) # shopping not going | |
| database$shop_1 <- ifelse(database$act2==4,1,0) # shopping less than 2 times | |
| database$shop_many <- ifelse(database$act2==2|database$act2==3,1,0) # shopping more than 3 times | |
| database$shop_tottori <- ifelse(database$place2==1,0,1) | |
| database$shop_car <- ifelse(database$to2==1|database$to2==2,1,0) | |
| ### create weights column for use in EM algorithm | |
| database$weights = 1 | |
| # ################################################################# # | |
| #### DEFINE MODEL PARAMETERS #### | |
| # ################################################################# # | |
| ### Vector of parameters, including any that are kept fixed in estimation | |
| apollo_beta = c(asc_1_a = 0, | |
| asc_1_b = 0, | |
| asc_1_c = 0, | |
| asc_2_a = 0, | |
| asc_2_b = 0, | |
| asc_2_c = 0, | |
| beta_wt_a = -0.4, | |
| beta_wt_b = -0.3, | |
| beta_wt_c = -0.5, | |
| beta_iv_a = -0.2, | |
| beta_iv_b = -0.13, | |
| beta_iv_c = -0.4, | |
| beta_md_a = -0.23, | |
| beta_md_b = -1, | |
| beta_md_c = -0.23, | |
| beta_bw_a = -0.3, | |
| beta_bw_b = -0.12, | |
| beta_bw_c = -0.145, | |
| beta_tc_a = -0.05, | |
| beta_tc_b = -0.02, | |
| beta_tc_c = -0.15, | |
| delta_a = 0, | |
| delta_c = 0, | |
| gamma_elderly_a = 0, | |
| gamma_commuter_a = 0, | |
| gamma_elderly_c = 0, | |
| gamma_commuter_c= 0, | |
| delta_b = 0, | |
| gamma_elderly_b = 0, | |
| gamma_commuter_b = 0, | |
| gamma_caruser_a = 0, | |
| gamma_caruser_b = 0, | |
| gamma_caruser_c = 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("asc_2_a", "asc_2_b", "asc_2_c","delta_c","gamma_commuter_c","gamma_elderly_c", "gamma_caruser_c") | |
| # ################################################################# # | |
| #### DEFINE LATENT CLASS COMPONENTS #### | |
| # ################################################################# # | |
| apollo_lcPars=function(apollo_beta, apollo_inputs){ | |
| lcpars = list() | |
| lcpars[["asc_1" ]] = list( asc_1_a, asc_1_b, asc_1_c) | |
| lcpars[["asc_2" ]] = list( asc_2_a, asc_2_b, asc_2_c) | |
| lcpars[["beta_wt"]] = list(beta_wt_a, beta_wt_b, beta_wt_c) | |
| lcpars[["beta_iv"]] = list(beta_iv_a, beta_iv_b, beta_iv_c) | |
| lcpars[["beta_md"]] = list(beta_md_a, beta_md_b, beta_md_c) | |
| lcpars[["beta_bw"]] = list(beta_bw_a, beta_bw_b, beta_bw_c) | |
| lcpars[["beta_tc"]] = list(beta_tc_a, beta_tc_b, beta_tc_c) | |
| V=list() | |
| V[["class_a"]] = delta_a + gamma_elderly_a*elderly + gamma_commuter_a*commuter + gamma_caruser_a * carowner | |
| V[["class_b"]] = delta_b + gamma_elderly_b*elderly + gamma_commuter_b*commuter + gamma_caruser_b * carowner | |
| V[["class_c"]] = delta_c + gamma_elderly_c*elderly + gamma_commuter_c*commuter + gamma_caruser_c * carowner | |
| ### Settings for class allocation models | |
| classAlloc_settings = list( | |
| classes = c(class_a=1, class_b=2, class_c=3), | |
| 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=2, alt2=1), | |
| avail = list(alt1=1, alt2=1), | |
| choiceVar = choice | |
| ) | |
| ### Loop over classes | |
| for(s in 1:length(pi_values)){ | |
| ### Compute class-specific utilities | |
| V=list() | |
| V[['alt1']] = asc_1[[s]] + beta_wt[[s]]*WT1 + beta_iv[[s]]*IV1 + beta_bw[[s]]*BW1 + beta_md[[s]]*DT1 + beta_tc[[s]]*TC1 | |
| V[['alt2']] = asc_2[[s]] + beta_wt[[s]]*WT2 + beta_iv[[s]]*IV2 + beta_bw[[s]]*BW2 + beta_md[[s]]*DT2 + beta_tc[[s]]*TC2 | |
| 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) | |
| } | |
| # ################################################################# # | |
| #### EM ESTIMATION FOR COVARIANCE MATRIX #### | |
| # ################################################################# # | |
| model=apollo_lcEM(apollo_beta, apollo_fixed, apollo_probabilities, | |
| apollo_inputs, lcEM_settings = list(EMmaxIterations=200)) | |
| # ----------------------------------------------------------------- # | |
| #---- FORMATTED OUTPUT (TO SCREEN) ---- | |
| # ----------------------------------------------------------------- # | |
| apollo_modelOutput(model) | |
| # ----------------------------------------------------------------- # | |
| #---- FORMATTED OUTPUT (TO FILE, using model name) ---- | |
| # ----------------------------------------------------------------- # | |
| #apollo_saveOutput(model) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Result:
Model run at : 2022-05-17 16:13:49
Estimation method : EM algorithm (bfgs) -> Maximum likelihood (bfgs)
Model diagnosis : successful convergence
Number of individuals : 1165
Number of rows in database : 4293
Number of modelled outcomes : 4293
Number of cores used : 1
Model without mixing
LL(start) : -2848.58
LL(0, whole model) : -2975.68
LL(C, whole model) : -2888.34
LL(final, whole model) : -1910.58
Rho-square (0) : 0.3579
Adj.Rho-square (0) : 0.3492
Rho-square (C) : 0.3385
Adj.Rho-square (C) : 0.3295
AIC : 3873.16
BIC : 4038.65
LL(0,Class_1) : -2975.68
LL(final,Class_1) : -5373.64
LL(0,Class_2) : -2975.68
LL(final,Class_2) : -16094.25
LL(0,Class_3) : -2975.68
LL(final,Class_3) : -7714.03
Estimated parameters : 26
Time taken (hh:mm:ss) : 00:02:13.96
pre-estimation : 00:00:4.48
estimation : 00:01:58.94
post-estimation : 00:00:10.55
Iterations : 90 (EM) & 35 (bfgs)
Min abs eigenvalue of Hessian : 0.043924
Some eigenvalues of Hessian are positive, indicating potential problems!
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_1_a 4.597119 1.276384 3.60167 1.509443 3.04557
asc_1_b -27.856248 NaN NaN 5.009692 -5.56047
asc_1_c -15.645011 NaN NaN 38.244792 -0.40908
asc_2_a 0.000000 NA NA NA NA
asc_2_b 0.000000 NA NA NA NA
asc_2_c 0.000000 NA NA NA NA
beta_wt_a -0.005569 0.082469 -0.06753 0.089096 -0.06251
beta_wt_b -2.070511 NaN NaN 0.421405 -4.91335
beta_wt_c -2.315043 NaN NaN 4.947627 -0.46791
beta_iv_a -0.051867 0.024606 -2.10786 0.025410 -2.04120
beta_iv_b -0.054558 0.011591 -4.70710 0.009526 -5.72718
beta_iv_c -1.031415 NaN NaN 1.986223 -0.51928
beta_md_a -0.214077 0.090520 -2.36496 0.103065 -2.07710
beta_md_b -1.174756 NaN NaN 0.278043 -4.22509
beta_md_c 0.608398 NaN NaN 0.456473 1.33282
beta_bw_a -0.182590 0.196766 -0.92795 0.207327 -0.88069
beta_bw_b -2.996197 NaN NaN 0.587952 -5.09599
beta_bw_c 1.723233 NaN NaN 1.046393 1.64683
beta_tc_a -0.020958 0.004238 -4.94579 0.004364 -4.80196
beta_tc_b -0.076817 NaN NaN 0.014767 -5.20179
beta_tc_c -0.071926 NaN NaN 0.025418 -2.82972
delta_a -0.238713 0.269398 -0.88610 0.288927 -0.82621
delta_c 0.000000 NA NA NA NA
gamma_elderly_a 1.158193 0.233494 4.96027 0.240400 4.81777
gamma_commuter_a 0.123817 0.247931 0.49940 0.265574 0.46622
gamma_elderly_c 0.000000 NA NA NA NA
gamma_commuter_c 0.000000 NA NA NA NA
delta_b -0.766038 0.314646 -2.43460 0.345228 -2.21893
gamma_elderly_b 1.375340 0.263549 5.21853 0.268600 5.12041
gamma_commuter_b -0.049522 0.289815 -0.17087 0.328972 -0.15054
gamma_caruser_a 0.252084 0.259699 0.97067 0.279820 0.90088
gamma_caruser_b 0.134804 0.292552 0.46078 0.317205 0.42497
gamma_caruser_c 0.000000 NA NA NA NA
Summary of class allocation for LC model component :
Mean prob.
Class_1 0.4629
Class_2 0.2642
Class_3 0.2729