adamSales
commented
6 years ago Not sure if this is exactly what you're asking, or if it's even right, but I'm also trying to fit a multiple-group model with poLCA right now. Actually, I'm just trying to get a p-value testing the null hypothesis that the same LCA model holds for two groups (journals tend to like that sort of thing).
The approach I'm using is a likelihood ratio test, with Wilks theorem (e.g. https://en.wikipedia.org/wiki/Likelihood-ratio_test)
If some of the probabilities are zero or one, then this might not hold.
I wrote a function that takes 3 arguments: lca1 is a model fit in group 1, lca2 is fit in group 2, and lca12 is to pooled data from the two groups.
lrt <- function(lca1,lca2,lca12){
llikDif <- (lca1$llik+lca2$llik)-lca12$llik
nparDif <- lca1$npar+lca2$npar-lca12$npar
pchisq(llikDif*2,nparDif,lower.tail=FALSE)
}If anyone knows whether this is right/wrong, I'd be happy to hear
I want to make multi-group LCA with poLCA. however, I am new to the syntax and cannot figure out how to do it. based on a paper by your lanza 2013 (https://www.ncbi.nlm.nih.gov/pubmed/21318625) that is enclosed,
In this paper, the authors have proposed two methods to investigate the effect of treatment on outcome in subgroups identified by LCA. 1) classify-analyze approach and 2) model-based approach. I want to implement the two approach with R and created following code: set.seed(8) probs <- list(matrix(c(0.6,0.2,0.2, 0.6,0.3,0.1, 0.3,0.1,0.6 ),ncol=3,byrow=TRUE), # Y1 matrix(c(0.2,0.8, 0.7,0.3, 0.3,0.7 ),ncol=2,byrow=TRUE), # Y2 matrix(c(0.3,0.6,0.1, 0.1,0.3,0.6, 0.3,0.6,0.1 ),ncol=3,byrow=TRUE), # Y3 matrix(c(0.1,0.1,0.5,0.3, 0.5,0.3,0.1,0.1, 0.3,0.1,0.1,0.5),ncol=4,byrow=TRUE), # Y4 matrix(c(0.1,0.2,0.7, 0.1,0.8,0.1, 0.8,0.1,0.1 ),ncol=3,byrow=TRUE)) # Y5 simdat <- poLCA.simdata(N=1000,probs,P=c(0.2,0.3,0.5)) trt<-as.factor(sample(c("trt","ctrl"),replace=T,size=1000)) z <- 1 - as.numeric(trt)-2simdat$trueclass+0.5as.numeric(trt)simdat$trueclass pr <- 1/(1+exp(-z)) outcome <- rbinom(1000,1,pr) dat<-data.frame(simdat$dat,trt=trt,outcome=outcome)
classify-analyze approach
f1 <- cbind(Y1,Y2,Y3,Y4,Y5)~1 lc1 <- poLCA(f1,simdat$dat,nclass=3,nrep=5) mod<-glm(outcome~trt*as.factor(lc1$predclass), family="binomial") summary(mod)
model based approach
f2<-cbind(Y1,Y2,Y3,Y4,Y5)~outcome dat.trt<-dat[dat$trt=="trt",] dat.ctrl<-dat[dat$trt=="ctrl",] lc2.trt<-poLCA(f2,dat.trt,nclass=3,nrep=5) lc2.ctrl<-poLCA(f2,dat.ctrl,nclass=3,nrep=5) table(lc2.trt$predclass,dat.trt$outcome) prop<-rbind(ctrl=prop.table(table(lc2.ctrl$predclass,dat.ctrl$outcome),1)[4:6], trt=prop.table(table(lc2.trt$predclass,dat.trt$outcome),1)[4:6]) colnames(prop)<-c('class 1',"class 2","class 3") barplot(prop,beside =T, legend.text=c('ctrl',"trt"))
however, it appears that the model-based approach is wrong. How can I implement the model-based approach with poLCA?