Bvr

[daob]

Hi Drew,

As requested, a pull request including documentation. Function seems to work fine with covariates also.

I removed the rounding from expected counts because this can mess with downstream statistics, hope that's OK. Also I added unit tests - feel free to ignore those if you don't want them.

Best,
Daniel

[zh-zhang1984]

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 Lanza 2013 (https://www.ncbi.nlm.nih.gov/pubmed/21318625) that is enclosed,
In this paper, the author 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(outcometrtas.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?