Latent Class Model - Effect coded attributes

[catalinaposadab]

Hello,

I am running a LCM based on the example LC_no_covariates. Two of my variables are in effect coding, I estimated a MNL and MMNL model I used the "Effects coding constraint" that is used in the MNL_SP_effects example and this worked fine. However, with the LC I have not managed to integrate the constraints into my code, I have tried to do it in multiple ways, but I always get errors and have not been able to estimate the model. I was wondering if you know how I should include the Effect coding constraints or if I have to transform the data and do it in a different way. Thank you!

These are the constraints that I used in my other models:
### Effects coding constraint 1
b_sps = - b_land1 - b_land2
### Effects coding constraint 2
b_nobn = - b_col - b_ind
and I thought for the LCM I could include these :
beta_sps_j = - beta_land1_j - beta_land2_j
beta_sps_k = - beta_land1_k - beta_land2_k
beta_nobn_j = - beta_col_j - beta_ind_j
beta_nobn_k = - beta_col_k - beta_ind_k

Here is the rest of my code:
# \### LOAD LIBRARY AND DEFINE CORE SETTINGS
### Clear memory
rm(list = ls())

### Initialise code
library(apollo)
apollo_initialise()

### Set core controls
apollo_control = list( modelName = "LCM_simple_total",
modelDescr = "Simple LC Logit model on total data",
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)
load("~/R DCE/Data/rdata/anonymised_data.Rdata")
database = anonymised_data

#### DEFINE MODEL PARAMETERS

### Vector of parameters, including any that are kept fixed in estimation
### k:a and j:b

apollo_beta = c(asc_sq = 0.297499,
beta_price_k = 0.16677,
beta_ta_k = 0.004224,
beta_land1_k = 0.176853,
beta_land2_k = -0.407689,
beta_ind_k = 0.005,
beta_col_k = 0.0001,
beta_price_j = 0.16677,
beta_ta_j = 0.004224,
beta_land1_j = 0.176853,
beta_land2_j = -0.407689,
beta_ind_j = 0.005,
beta_col_j = 0.0001,
delta_j=0,
delta_k=0)

apollo_fixed = c()
apollo_beta_fixed = c()


# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()

lcpars[["beta_price"]] = list(beta_price_k, beta_price_j)
lcpars[["beta_ta"]] = list(beta_ta_k, beta_ta_j)
lcpars[["beta_sps"]] = list(-apollo_beta["beta_land1_k"] - apollo_beta["beta_land2_k"],
-apollo_beta["beta_land1_j"] - apollo_beta["beta_land2_j"])
lcpars[["beta_land1"]] = list(beta_land1_k, beta_land1_j)
lcpars[["beta_land2"]] = list(beta_land2_k, beta_land2_j)
lcpars[["beta_nobn"]] = list(-apollo_beta["beta_col_k"] - apollo_beta["beta_ind_k"],
-apollo_beta["beta_col_j"] - apollo_beta["beta_ind_j"])
lcpars[["beta_col"]] = list(beta_col_k, beta_col_j)
lcpars[["beta_land2"]] = list(beta_ind_k, beta_ind_j)


V=list()
V[["class_k"]] = delta_k
V[["class_j"]] = delta_j

classAlloc_settings = list(
classes = c(class_k=1, class_j=2),
utilities = V
)

lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)

return(lcpars)
}

# ################################################################# #
#### 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(ProgramA=1, ProgramB=2,
NoProgram=3),
choiceVar = choice
)
### Loop over classes
for(s in 1:2){

###utilities

V = list()
V[["ProgramA"]] = beta_sps* (A.landdist==0) + beta_land1*
(A.landdist==1) + beta_land2* (A.landdist==2) + beta_nobn* (A.bonus==0)+
beta_col* (A.bonus==1) + beta_ind* (A.bonus==2) + beta_ta* A.trainning +
beta_price* A.compensation
V[["ProgramB"]] = beta_sps* (B.landdist==0) +
beta_land1* (B.landdist==1) + beta_land2* (B.landdist==2) + beta_nobn*
(B.bonus==0) + beta_col* (B.bonus==1) + beta_ind* (B.bonus==2) + beta_ta*
B.trainning + beta_price* B.compensation
V[["NoProgram"]] = asc_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)

[stephanehess]

Hi

it seems like you are not using the index across classes for the parameters - can you first get a model set up with dummy coding that runs, and then we can help you further with that. Please have a look at the LC example online

Stephane

[catalinaposadab]

Hi,

Thank you very much for your answer, I did had some errors in the code I posted. Here is a similar code with 3 classes and dummy coding that runs, and that I would like to adapt for effects coding. Thank you!


# \### LOAD LIBRARY AND DEFINE CORE SETTINGS
### Clear memory

### Initialise code
library(apollo)
apollo_initialise()

### Set core controls
apollo_control = list( modelName = "LCM_simple_total",
modelDescr = "Simple LC Logit model on total data",
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)
database = dummy_anonymised


# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #



#### DEFINE MODEL PARAMETERS

### Vector of parameters, including any that are kept fixed in estimation
### k:a and j:b

apollo_beta = c(asc_sq_k = 0.297499,
asc_sq_j = 0.297499,
asc_sq_q = 0.297499,
beta_price_k = 0.16677,
beta_ta_k = 0.004224,
beta_land1_k = 0.176853,
beta_land2_k = -0.407689,
beta_ind_k = 0.005,
beta_col_k = 0.0001,
beta_price_j = 0.16677,
beta_ta_j = 0.004224,
beta_land1_j = 0.176853,
beta_land2_j = -0.407689,
beta_ind_j = 0.005,
beta_col_j = 0.0001,
beta_price_q = 0.16677,
beta_ta_q = 0.004224,
beta_land1_q = 0.176853,
beta_land2_q = -0.407689,
beta_ind_q = 0.005,
beta_col_q = 0.0001,
delta_j=0,
delta_q=0,
delta_k=0)

apollo_fixed = c("delta_q" , "delta_k")


# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
lcpars = list()

lcpars[["beta_price"]] = list(beta_price_k, beta_price_j, beta_price_q)
lcpars[["beta_ta"]] = list(beta_ta_k, beta_ta_j, beta_ta_q)
lcpars[["beta_land1"]] = list(beta_land1_k, beta_land1_j, beta_land1_q)
lcpars[["beta_land2"]] = list(beta_land2_k, beta_land2_j, beta_land2_q)
lcpars[["beta_ind"]] = list(beta_ind_k, beta_ind_j, beta_ind_q)
lcpars[["beta_col"]] = list(beta_col_k, beta_col_j, beta_col_q)
lcpars[["asc_sq"]] = list(asc_sq_k, asc_sq_j, asc_sq_q)



V=list()
V[["class_k"]] = delta_k
V[["class_j"]] = delta_j
V[["class_q"]] = delta_q


classAlloc_settings = list(
classes = c(class_k=1, class_j=2, class_q=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(ProgramA=1, ProgramB=2,
NoProgram=3),
choiceVar = choice
)
### Loop over classes
for(s in 1:3){

###utilities

V = list()
V[["ProgramA"]] =
beta_land1[[s]]* (A.land1==1) + beta_land2[[s]]* (A.land2==1)
+beta_col[[s]]* (A.collective==1) + beta_ind[[s]]* (A.individual==1) + beta_ta[[s]]* A.trainning +
beta_price[[s]]* A.compensation
V[["ProgramB"]] = beta_land1[[s]]* (B.land1==1) + beta_land2[[s]]* (B.land2==1) + beta_col[[s]]* (B.collective==1) + beta_ind[[s]]* (B.individual==1) + beta_ta[[s]]*
B.trainning + beta_price[[s]]* B.compensation
V[["NoProgram"]] = asc_sq[[s]]


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)

And the last part of the output is:
Estimated parameters with approximate standard errors from BHHH matrix:
Estimate BHHH se BHH t-ratio
asc_sq_k -2.361357 0.148232 -15.9301
asc_sq_j -0.460037 0.150163 -3.0636
asc_sq_q -5.820586 0.548437 -10.6130
beta_price_k 0.010792 0.016938 0.6371
beta_ta_k -0.026111 0.007755 -3.3671
beta_land1_k 0.334300 0.082724 4.0412
beta_land2_k 0.383821 0.090826 4.2259
beta_ind_k 0.761889 0.198213 3.8438
beta_col_k -0.053093 0.168126 -0.3158
beta_price_j 0.071689 0.039330 1.8228
beta_ta_j -0.007622 0.019113 -0.3988
beta_land1_j -1.390761 0.160255 -8.6784
beta_land2_j -3.661411 0.281364 -13.0131
beta_ind_j 0.327961 0.542066 0.6050
beta_col_j -2.365200 0.425853 -5.5540
beta_price_q 0.007343 0.044836 0.1638
beta_ta_q -0.047389 0.025519 -1.8570
beta_land1_q -0.766639 0.265867 -2.8835
beta_land2_q -3.774822 0.368112 -10.2545
beta_ind_q 1.159673 0.590248 1.9647
beta_col_q -0.198829 0.543497 -0.3658
delta_j -0.742564 0.181431 -4.0928
delta_q 0.000000 NA NA
delta_k 0.000000 NA NA

Final LL: -1732.1873


Summary of class allocation for model component :
Mean prob.
Class_1 0.4039
Class_2 0.1922
Class_3 0.4039

Calculating log-likelihood at equal shares (LL(0)) for applicable models...
Calculating log-likelihood at observed shares from estimation data (LL(c)) for applicable models...
Calculating LL of each model component...
Calculating other model fit measures
Computing covariance matrix using analytical gradient.
0%....25%....50%....75%....100%
Negative definite Hessian with maximum eigenvalue: -1.498396

[stephanehess]

Can you tell me one attribute where you want to use effects coding, and what the levels are? thanks

[catalinaposadab]

For instance, I have a "Bonus" attribute in effects coding from my Ngene design.
The reference level at 0 is "No bonus", when bonus equals 1 it is a "collective bonus", and when bonus is equal to 2 it is an "individual bonus".

In my MNL model, this is the constraint I use for that attribute: b_nobn = - b_col - b_ind

Thanks again!

[stephanehess]

so what you should be able to do is to include this in apollo_lcPars

Code: Select all

lcpars[["b_nobn"]] = list(b_col_k - b_ind_k,b_col_j - b_ind_j,b_col_q - b_ind_q)

or alternatively, inside the loop in apollo_probabilities, you could use

Code: Select all

b_nobn[[s]]=b_col_k[[s]] - b_ind_k[[s]]

[catalinaposadab]

Hi Stephane,

That makes sense and it works. Thank you very much !

Best,

Catalina