Hi all
I am working on creating a general framework to work with mxiture models, specifically with Latent Class Analysis (LCA) and its extensions.
Now, I am working on a LCA model, using the bernoulli_logit_lpmf functions from Stan. This model so far works well when I estimate 2 clases, but when i estimate more than classes it generates bimodal posteriors.
I have looked at the post about Mixture Models Identification, and the solutions there dont work so far.
What I have done so far is to constraint the mixing proporcions to be ordered, as in this example ordered simplex. In addition to this, I have defined hyper priors for the parameters.
With these definition, the model provides good results for 2 classes, but cant figure out how to define it for larger number of classes.
Here is the Stan model so far an 9 item example
data{
int<lower=1> I; // number of items
int<lower=1> J; // number of respondents
int<lower=0, upper=1> y[J,I]; // score for obs n
int<lower=1> C; // # of attribute profiles (latent classes)
}
parameters {
matrix[I,C] pi; // log-odds of answering 1 for each item in each class
positive_ordered[C] lambda;
real mu_pi; // prior mean for the log-odds
real<lower=0> sigma_pi; // prior sd for the log-odds
real<lower=0> alpha_lambda; // prior alpha for lambda
real<lower=0> beta_lambda; // prior beta for lambda
}
transformed parameters {
vector[C] Ys[J]; // sum log for each subject for each class
matrix[I,C] pi_probs; // probability of answering 1 for each item on each class
simplex[C] nu = lambda / sum(lambda); // ordered mixing proportions
for(c in 1:C){
for(j in 1:J){
Ys[j,c] = bernoulli_logit_lpmf(y[j,1] | pi[1,c])+
bernoulli_logit_lpmf(y[j,2] | pi[2,c])+
bernoulli_logit_lpmf(y[j,3] | pi[3,c])+
bernoulli_logit_lpmf(y[j,4] | pi[4,c])+
bernoulli_logit_lpmf(y[j,5] | pi[5,c])+
bernoulli_logit_lpmf(y[j,6] | pi[6,c])+
bernoulli_logit_lpmf(y[j,7] | pi[7,c])+
bernoulli_logit_lpmf(y[j,8] | pi[8,c])+
bernoulli_logit_lpmf(y[j,9] | pi[9,c]);
}
}
pi_probs = inv_logit(pi);
}
model{
// Priors
for(i in 1:I){
pi[i] ~ normal(mu_pi, sigma_pi);
}
target += gamma_lpdf(lambda | alpha_lambda, beta_lambda);
mu_pi ~ normal(0,2);
sigma_pi ~ cauchy(0,2.5);
alpha_lambda ~ cauchy(1,2.5);
beta_lambda ~ cauchy(1,2.5);
target += log_mix(nu, Ys);
}
And here is the R code to run this example with a data set
library(poLCA)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
### https://stats.idre.ucla.edu/mplus/dae/latent-class-analysis/
dat <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/lca1.dat",
header=F)
colnames(dat) <- c("id", paste0("item",1:9))
head(dat)
apply(dat[,2:10], 2, table)
## 9 item data
dat2 <- dat[,2:10]
head(dat2)
apply(dat2, 2, table)
wide_data <- list(
I = ncol(dat2),
J = nrow(dat2),
y = dat2,
C = 3
)
pars <- c("nu","pi","pi_probs")
iter <- 6000
warmup <- iter - 1000
wide_model_bn <- stan(model_code=lca_wide_bn,
data = wide_data,
pars=pars,
chains = 3,
iter = iter,
warmup = warmup,
control = list(adapt_delta = 0.99))
wide_model_bn
sm_lca <- data.frame(summary(wide_model_bn)$summary)
sm_lca
bayesplot::mcmc_areas(wide_model_bn,
pars = rownames(sm_lca)[4:10])
### Here is the solution from the poLCA package
ch3 <- poLCA(cbind(item1, item2, item3, item4, item5, item6, item7, item8, item9)~1,
dat2+1,nclass=3, maxiter=20000)
ch3
Some of the bimodal posteriors look like this
One last thing I tried to identified the model, but made no difference, was to restric the log-odds for 1 item to be ordered, intending for it to define the direction of the classifications.
Any help would be appreciated!
