Nothing
R/MLCAfunctions.R
In multilevLCA: Estimates and Plots Single-Level and Multilevel Latent Class Models
Defines functions
plot.multiLCA
print.multiLCA
check_inputs2
check_inputs1
update_YmY_nonnum
update_YmY
clean_cov
clean_output5
clean_output4
clean_output3
clean_output2
clean_output1
sel_other_poly
lukosel_fun_poly
simultsel_fun_poly
sel_other
lukosel_fun
simultsel_fun
meas_Init_poly
LCA_fast_init_whigh_poly
LCA_fast_init_wcov_poly
LCA_fast_init_poly
meas_Init
LCA_fast_init_whigh
LCA_fast_init_wcov
LCA_fast_init
multiLCA
.onAttach
Documented in
multiLCA
plot.multiLCA
.onAttach <- function(libname, pkgname){
packageStartupMessage("multilevLCA: Estimates and plots single- and multilevel latent class models.")
}
multiLCA = function(data, Y, iT, id_high = NULL, iM = NULL, Z = NULL, Zh = NULL,
extout = FALSE, dataout = TRUE, kmea = TRUE, sequential = TRUE,
numFreeCores = 2, maxIter = 1e3, tol = 1e-8,
reord = TRUE, fixedpars = 1,
NRmaxit = 100, NRtol = 1e-6, verbose = TRUE){
#
if(verbose){
pb=txtProgressBar(char="Cleaning data...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
check_inputs1(data,Y,iT,id_high,iM,Z,Zh)
nrow_data = nrow(data)
data = data[complete.cases(data[,Y]),c(Y,id_high,Z,Zh)]
approach = check_inputs2(data,Y,iT,id_high,iM,Z,Zh)
reord = as.numeric(reord)
fixed = fixedpars
if((fixed==2)&is.null(iM)) fixed=fixedpars=1
if(!is.null(id_high)){
id_high_levs = sort(unique(data[,id_high]))
data[,id_high] = as.numeric(factor(data[,id_high]))
data = data[order(data[,id_high]),]
id_high_name = id_high
id_high = data[,id_high]
vNj = table(id_high)
}
mY = as.matrix(data[,Y])
ivItemcat = apply(mY,2,function(x){length(unique(x))})
itemchar = any(apply(data[,Y],2,function(x){!is.numeric(x)}))
if(any(ivItemcat>2)&!itemchar){
mY = update_YmY(mY,Y,ivItemcat)
Y = mY$Y
mY = mY$mY
} else if(itemchar){
mY = update_YmY_nonnum(mY,Y,ivItemcat)
Y = mY$Y
mY = mY$mY
}
if(!is.null(Z)){
mZ = clean_cov(data,Z)
Z = mZ$covnam
mZ = mZ$cov
}
if(!is.null(Zh)){
mZh = clean_cov(data,Zh)
Zh = mZh$covnam
mZh = mZh$cov
}
if(approach=="direct"){
#
if(verbose){
if(fixed==0){
pb=txtProgressBar(char="Fitting LC model...",width=1)
} else{
pb=txtProgressBar(char="Fitting measurement model...",width=1)
}
setTxtProgressBar(pb,1)
close(pb)
}
#
if(any(ivItemcat>2)){
#
if(is.null(id_high)&is.null(Z)&is.null(Zh)){
out = LCA_fast_init_poly(mY,iT,ivItemcat,kmea,maxIter,tol,reord)
out = clean_output1(out,Y,iT,length(Y),extout,dataout,ivItemcat)
} else if(is.null(id_high)&!is.null(Z)&is.null(Zh)){
out = LCA_fast_init_wcov_poly(mY,mZ,iT,ivItemcat,kmea,maxIter,tol,fixed,reord,
NRtol,NRmaxit,verbose)
out = clean_output2(out,Y,iT,c("Intercept",Z),length(Y),length(Z)+1,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&is.null(Z)&is.null(Zh)){
init = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
out = MLTLCA_poly(mY,vNj,init$vOmega_start,init$mPi_start,init$mPhi_start,
ivItemcat,maxIter,tol,reord)
out = clean_output3(out,Y,iT,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&is.null(Zh)){
init1 = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
init2 = MLTLCA_poly(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
vOmega_start = init2$vOmega
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
#
mY = mY[complete.cases(mZ),]
id_high = id_high[complete.cases(mZ)]
vNj = table(id_high)
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_cov_poly(mY,mZ,vNj,vOmega_start,cGamma_start,mPhi_start,mStep1Var,
ivItemcat,maxIter,tol,fixedpars,1,NRtol,NRmaxit)
out = clean_output4(out,Y,iT,iM,c("Intercept",Z),mY,mZ,id_high,length(Y),P,id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&!is.null(Zh)){
init1 = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
init2 = MLTLCA_poly(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
P_high = ncol(mZh)
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
mDelta_start = matrix(0,iM-1,P_high)
mDelta_start[,1] = init2$vDelta
#
nomissing = complete.cases(mZ)&complete.cases(mZh)
mY = mY[nomissing,]
id_high = id_high[nomissing]
vNj = table(id_high)
mZ = mZ[nomissing,]
mZh = mZh[nomissing,]
mZh = mZh[!duplicated(id_high),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_covlowhigh_poly(mY,mZ,mZh,vNj,mDelta_start,cGamma_start,mPhi_start,
mStep1Var,ivItemcat,maxIter,tol,fixedpars,1,NRtol,NRmaxit)
out = clean_output5(out,Y,iT,iM,c("Intercept",Z),c("Intercept",Zh),mY,mZ,mZh,id_high,length(Y),P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat)
}
} else{
if(is.null(id_high)&is.null(Z)&is.null(Zh)){
out = LCA_fast_init(mY,iT,kmea,maxIter,tol,reord)
out = clean_output1(out,Y,iT,length(Y),extout,dataout,ivItemcat)
} else if(is.null(id_high)&!is.null(Z)&is.null(Zh)){
out = LCA_fast_init_wcov(mY,mZ,iT,kmea,maxIter,tol,fixed,reord,
NRtol,NRmaxit,verbose)
out = clean_output2(out,Y,iT,c("Intercept",Z),length(Y),length(Z)+1,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&is.null(Z)&is.null(Zh)){
init = meas_Init(mY,id_high,vNj,iM,iT,kmea)
out = MLTLCA(mY,vNj,init$vOmega_start,init$mPi_start,init$mPhi_start,
maxIter,tol,reord)
out = clean_output3(out,Y,iT,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&is.null(Zh)){
init1 = meas_Init(mY,id_high,vNj,iM,iT,kmea)
init2 = MLTLCA(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
maxIter,tol,reord)
P = ncol(mZ)
vOmega_start = init2$vOmega
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
#
mY = mY[complete.cases(mZ),]
id_high = id_high[complete.cases(mZ)]
vNj = table(id_high)
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_cov(mY,mZ,vNj,vOmega_start,cGamma_start,mPhi_start,mStep1Var,
maxIter,tol,fixedpars,NRtol,NRmaxit)
out = clean_output4(out,Y,iT,iM,c("Intercept",Z),mY,mZ,id_high,length(Y),P,id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&!is.null(Zh)){
init1 = meas_Init(mY,id_high,vNj,iM,iT,kmea)
init2 = MLTLCA(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
maxIter,tol,reord)
P = ncol(mZ)
P_high = ncol(mZh)
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
mDelta_start = matrix(0,iM-1,P_high)
mDelta_start[,1] = init2$vDelta
#
nomissing = complete.cases(mZ)&complete.cases(mZh)
mY = mY[nomissing,]
id_high = id_high[nomissing]
vNj = table(id_high)
mZ = mZ[nomissing,]
mZh = mZh[nomissing,]
mZh = mZh[!duplicated(id_high),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_covlowhigh(mY,mZ,mZh,vNj,mDelta_start,cGamma_start,mPhi_start,
mStep1Var,maxIter,tol,fixedpars,NRtol,NRmaxit)
out = clean_output5(out,Y,iT,iM,c("Intercept",Z),c("Intercept",Zh),mY,mZ,mZh,id_high,length(Y),P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat)
}
}
} else if(approach=="model selection on low"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output1(out$outmuLCA,Y,out$iT_best,length(Y),extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output1(out$outmuLCA,Y,out$iT_best,length(Y),extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on low with high"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,out$iT_best,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,out$iT_best,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on high"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iM),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iM_best,dimnames=list("iM=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,iT,out$iM_best,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iM),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iM_best,dimnames=list("iM=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,iT,out$iM_best,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on low and high"){
if(sequential){
if(any(ivItemcat>2)){
out = lukosel_fun_poly(mY,id_high,iT,iM,ivItemcat,verbose)
step1mat = matrix(as.character(round(unlist(out[4:8]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step1mat[step1mat=="Inf"|step1mat=="-Inf"|is.na(step1mat)] = "-"
step2mat = matrix(as.character(round(unlist(out[9:13]),2)),length(iM),5,
dimnames=list(paste0("iT*,","iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step2mat[step2mat=="Inf"|step2mat=="-Inf"|is.na(step2mat)] = "-"
step3mat = matrix(as.character(round(unlist(out[14:18]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM*"),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step3mat[step3mat=="Inf"|step3mat=="-Inf"|is.na(step3mat)] = "-"
optimal = matrix(c(out$iT_opt,out$iM_opt),dimnames=list(c("iT=","iM="),""))
list_luk = list(step1 = step1mat,step2 = step2mat,step3 = step3mat, optimal = optimal)
if(out$iM_opt>1){
out = clean_output3(out$outmuLCA_step3,Y,out$iT_opt,out$iM_opt,mY,id_high,length(Y),
id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(out$iM_opt==1){
out = clean_output1(out$outmuLCA_step3,Y,out$iT_opt,length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_luk
if(verbose)print(noquote(list_luk))
} else{
out = lukosel_fun(mY,id_high,iT,iM,verbose)
step1mat = matrix(as.character(round(unlist(out[4:8]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step1mat[step1mat=="Inf"|step1mat=="-Inf"|is.na(step1mat)] = "-"
step2mat = matrix(as.character(round(unlist(out[9:13]),2)),length(iM),5,
dimnames=list(paste0("iT*,","iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step2mat[step2mat=="Inf"|step2mat=="-Inf"|is.na(step2mat)] = "-"
step3mat = matrix(as.character(round(unlist(out[14:18]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM*"),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step3mat[step3mat=="Inf"|step3mat=="-Inf"|is.na(step3mat)] = "-"
optimal = matrix(c(out$iT_opt,out$iM_opt),dimnames=list(c("iT=","iM="),""))
list_luk = list(step1 = step1mat,step2 = step2mat,step3 = step3mat, optimal = optimal)
if(out$iM_opt>1){
out = clean_output3(out$outmuLCA_step3,Y,out$iT_opt,out$iM_opt,mY,id_high,length(Y),
id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(out$iM_opt==1){
out = clean_output1(out$outmuLCA_step3,Y,out$iT_opt,length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_luk
if(verbose)print(noquote(list_luk))
}
} else{
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement models simultaneously...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(any(ivItemcat>2)){
mYfoo = mY
id_highfoo = id_high
ivItemcatfoo = ivItemcat
select_mat = as.matrix(tidyr::expand_grid(iT,iM))
nmod = nrow(select_mat)
iV = parallel::detectCores()-numFreeCores
cluster = parallel::makeCluster(iV)
parallel::clusterExport(cluster,c("mYfoo","id_highfoo","ivItemcatfoo","select_mat","nmod"),
envir = environment())
parallel::clusterExport(cluster,
unclass(lsf.str(envir = asNamespace("multilevLCA"),
all = T)),
envir = as.environment(asNamespace("multilevLCA"))
)
parallel::clusterEvalQ(cluster, {
library(clustMixType)
library(multilevLCA)
library(mclust)
library(tictoc)
library(numDeriv)
library(klaR)
library(tidyverse)
library(MASS)})
simultaneous_out = parallel::parLapply(cluster,1:nmod,
function(x){
iFoo = select_mat[x,]
iT_curr = iFoo[1]
iM_curr = iFoo[2]
out_simult = simultsel_fun_poly(mYfoo,id_highfoo,iT_curr,iM_curr,ivItemcatfoo)
BIClow = out_simult$BIClow
BIChigh = out_simult$BIChigh
AIC = out_simult$AIC
ICLlow = out_simult$ICL_BIClow
ICLhigh = out_simult$ICL_BIChigh
nout = c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")
out = c(BIClow,BIChigh,AIC,ICLlow,ICLhigh)
names(out) = nout
return(list(out=out,fit=out_simult$modFIT))
}
)
parallel::stopCluster(cluster)
mod_sel = apply(t(sapply(simultaneous_out,function(x){x$out})),
2,
function(x){round(as.numeric(x),2)})
rownames(mod_sel) = paste0("iT=",select_mat[,1],",iM=",select_mat[,2])
r_optimal = which.min(mod_sel[,1])
optimal = matrix(c(select_mat[r_optimal,1],select_mat[r_optimal,2]),dimnames=list(c("iT=","iM="),""))
mod_sel[mod_sel==Inf|mod_sel==-Inf|is.na(mod_sel)] = "-"
list_simul = list(model_selection=mod_sel,optimal=optimal)
if(select_mat[r_optimal,2]>1){
out = clean_output3(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],select_mat[r_optimal,2],
mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(select_mat[r_optimal,2]==1){
out = clean_output1(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_simul
if(verbose)print(noquote(list_simul))
} else{
mYfoo = mY
id_highfoo = id_high
select_mat = as.matrix(tidyr::expand_grid(iT,iM))
nmod = nrow(select_mat)
iV = parallel::detectCores()-numFreeCores
cluster = parallel::makeCluster(iV)
parallel::clusterExport(cluster,c("mYfoo","id_highfoo","select_mat","nmod"), envir = environment())
parallel::clusterExport(cluster,
unclass(lsf.str(envir = asNamespace("multilevLCA"),
all = T)),
envir = as.environment(asNamespace("multilevLCA"))
)
parallel::clusterEvalQ(cluster, {
library(clustMixType)
library(multilevLCA)
library(mclust)
library(tictoc)
library(numDeriv)
library(klaR)
library(tidyverse)
library(MASS)})
simultaneous_out = parallel::parLapply(cluster,1:nmod,
function(x){
iFoo = select_mat[x,]
iT_curr = iFoo[1]
iM_curr = iFoo[2]
out_simult = simultsel_fun(mYfoo,id_highfoo,iT_curr,iM_curr)
BIClow = out_simult$BIClow
BIChigh = out_simult$BIChigh
AIC = out_simult$AIC
ICLlow = out_simult$ICL_BIClow
ICLhigh = out_simult$ICL_BIChigh
nout = c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")
out = c(BIClow,BIChigh,AIC,ICLlow,ICLhigh)
names(out) = nout
return(list(out=out,fit=out_simult$modFIT))
}
)
parallel::stopCluster(cluster)
mod_sel = apply(t(sapply(simultaneous_out,function(x){x$out})),
2,
function(x){round(as.numeric(x),2)})
rownames(mod_sel) = paste0("iT=",select_mat[,1],",iM=",select_mat[,2])
r_optimal = which.min(mod_sel[,1])
optimal = matrix(c(select_mat[r_optimal,1],select_mat[r_optimal,2]),dimnames=list(c("iT=","iM="),""))
mod_sel[mod_sel==Inf|mod_sel==-Inf|is.na(mod_sel)] = "-"
list_simul = list(model_selection=mod_sel,optimal=optimal)
if(select_mat[r_optimal,2]>1){
out = clean_output3(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],select_mat[r_optimal,2],
mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(select_mat[r_optimal,2]==1){
out = clean_output1(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_simul
if(verbose)print(noquote(list_simul))
}
}
}
if(nrow(data)<nrow_data){
if(fixed==1|fixed==2){
warning(paste("Missing values in columns for indicators,",
"sample size for estimation of measurement model:",
nrow(data)),call.=FALSE)
} else if(fixed==0&is.null(Z)){
warning(paste("Missing values in columns for indicators,",
"sample size for estimation of measurement model:",
nrow(data)),call.=FALSE)
} else if(fixed==0&!is.null(Z)){
if(nrow(na.omit(mZ))<nrow(data)){
warning(paste("Missing values in columns for indicators",
"and columns for covariates,",
"sample size for estimation:",
nrow(na.omit(mZ))),call.=FALSE)
} else{
warning(paste("Missing values in columns for indicators,",
"sample size for estimation:",
nrow(data)),call.=FALSE)
}
}
}
if(approach=="direct"&!is.null(Z)){
if(nrow(na.omit(mZ))<nrow(data)){
if(fixed==1|fixed==2){
warning(paste("Missing values in columns for covariates,",
"sample size for estimation of structural model:",
nrow(na.omit(mZ))),call.=FALSE)
} else if(fixed==0){
if(nrow(data)==nrow_data){
warning(paste("Missing values in columns for covariates,",
"sample size for estimation of structural model:",
nrow(na.omit(mZ))),call.=FALSE)
}
}
}
}
out$call = match.call()
class(out) = "multiLCA"
return(out)
}
#
LCA_fast_init = function(mY, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
out$mY_unique = mY_unique
out$freq = freq
return(out)
}
#
LCA_fast_init_wcov = function(mY, mZ, iT, kmea = T, maxIter = 1e3, tol = 1e-8, fixed = 0, reord = 1,
NRtol = 1e-6, NRmaxit = 100, verbose){
# mZ must include a column of ones!
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
P = ncol(mZ)
mBeta_init = matrix(0,P,iT-1)
mBeta_init[1,] = out$alphas
Step1Var = out$Varmat
mY = mY[complete.cases(mZ),]
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
outcov = LCAcov(mY,mZ,iT,out$mPhi,mBeta_init,Step1Var,fixed,maxIter,
tol,NRtol,NRmaxit)
return(list(out=out,outcov=outcov,mY=mY,mZ=mZ))
}
#
LCA_fast_init_whigh = function(mY, id_high, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
group_by_all = NULL
mY_df = data.frame(mY,id_high)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
mY_aggr = mY_aggr[order(mY_aggr[,iHf-1]),]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-c(iHf-1,iHf)]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
mU = out$mU[rep(1:nrow(mY_unique),freq),]
return(list(out=out,mU=mU))
}
#
meas_Init = function(mY, id_high, vNj, iM, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# Fixed number of low- (iT) and high- (iM) level classes
iJ = length(vNj)
iN = dim(mY)[1]
iK = dim(mY)[2]
#
# Working out starting values at the low level first
#
out_LCA = LCA_fast_init_whigh(mY,id_high,iT,kmea,maxIter,tol,reord)
#
# Now turning to the higher level
#
lowlev_relclassprop = matrix(0,iJ,iT)
foo = 0
for(j in 1:iJ){
lowlev_relclassprop[j,] = colMeans(out_LCA$mU[foo+(1:vNj[j]),,drop=FALSE])
foo = foo +vNj[j]
}
high_out = kmeans(lowlev_relclassprop,centers = iM,iter.max=100,nstart=100)
vOmega_start = high_out$size/iJ
Wmodal_mat = vecTomatClass(high_out$cluster)
# Reordering in decreasing order
highreord = order(vOmega_start,decreasing = T)
Wmodal_mat_reord = Wmodal_mat[,highreord]
Wmodal = apply(Wmodal_mat_reord,1,which.max)
vOmega_start = vOmega_start[highreord]
index_indiv_high = rep(Wmodal,times=vNj)
#
mPhi_start = out_LCA$out$mPhi
mPi_fast = table(index_indiv_high,rep(out_LCA$out$vModalAssnm+1,times=out_LCA$out$freq))
mPi_start = t(mPi_fast/rowSums(mPi_fast))
if(nrow(mPi_start)!=iT){
warning(paste0("Initialization suggests less than iT=",iT," lower-level classes."),call.=FALSE)
mPi_start_replace = matrix(NA,nrow=iT,ncol=iM)
mPi_start_replace[as.numeric(rownames(mPi_start)),] = mPi_start
mPi_start_replace[is.na(mPi_start_replace)] = 1e-4
mPi_start_replace = apply(mPi_start_replace,2,function(x){x/sum(x)})
mPi_start = mPi_start_replace
}
#
return(list(vOmega_start=vOmega_start, mPi_start = mPi_start, mPhi_start = mPhi_start))
}
#
LCA_fast_init_poly = function(mY, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
out$mY_unique = mY_unique
out$freq = freq
return(out)
}
#
LCA_fast_init_wcov_poly = function(mY, mZ, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, fixed = 0, reord = 1,
NRtol = 1e-6, NRmaxit = 100, verbose){
# mZ must include a column of ones!
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
mBeta_init = matrix(0,P,iT-1)
mBeta_init[1,] = out$alphas
Step1Var = out$Varmat
mY = mY[complete.cases(mZ),]
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
outcov = LCAcov_poly(mY,mZ,iT,out$mPhi,mBeta_init,Step1Var,ivItemcat,fixed,maxIter,
tol,NRtol,NRmaxit)
return(list(out=out,outcov=outcov,mY=mY,mZ=mZ))
}
#
LCA_fast_init_whigh_poly = function(mY, id_high, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY,id_high)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
mY_aggr = mY_aggr[order(mY_aggr[,iHf-1]),]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-c(iHf-1,iHf)]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
mU = out$mU[rep(1:nrow(mY_unique),freq),]
return(list(out=out,mU=mU))
}
#
meas_Init_poly = function(mY, id_high, vNj, iM, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# Fixed number of low- (iT) and high- (iM) level classes
# ivItemcat is the vector of number of categories for each item
iJ = length(vNj)
iN = dim(mY)[1]
iK = dim(mY)[2]
#
# Working out starting values at the low level first
#
out_LCA = LCA_fast_init_whigh_poly(mY,id_high,iT,ivItemcat,kmea,maxIter,tol,reord)
#
# now turning to higher level
#
lowlev_relclassprop = matrix(0,iJ,iT)
foo = 0
for(j in 1:iJ){
lowlev_relclassprop[j,] = colMeans(out_LCA$mU[foo+(1:vNj[j]),,drop=FALSE])
foo = foo +vNj[j]
}
high_out = kmeans(lowlev_relclassprop,centers = iM,iter.max=100,nstart=100)
vOmega_start = high_out$size/iJ
Wmodal_mat = vecTomatClass(high_out$cluster)
# Reordering in decreasing order
highreord = order(vOmega_start,decreasing = T)
Wmodal_mat_reord = Wmodal_mat[,highreord]
Wmodal = apply(Wmodal_mat_reord,1,which.max)
vOmega_start = vOmega_start[highreord]
index_indiv_high = rep(Wmodal,times=vNj)
#
mPhi_start = out_LCA$out$mPhi
mPi_fast = table(index_indiv_high,rep(out_LCA$out$vModalAssnm+1,times=out_LCA$out$freq))
mPi_start = t(mPi_fast/rowSums(mPi_fast))
if(nrow(mPi_start)!=iT){
warning(paste0("Initialization suggests less than iT=",iT," lower-level classes."),call.=FALSE)
mPi_start_replace = matrix(NA,nrow=iT,ncol=iM)
mPi_start_replace[as.numeric(rownames(mPi_start)),] = mPi_start
mPi_start_replace[is.na(mPi_start_replace)] = 1e-4
mPi_start_replace = apply(mPi_start_replace,2,function(x){x/sum(x)})
mPi_start = mPi_start_replace
}
#
return(list(vOmega_start=vOmega_start, mPi_start = mPi_start, mPhi_start = mPhi_start))
}
#
#
#
simultsel_fun = function(mY,id_high,iT,iM){
LCAout = NULL
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
if(iM ==1 & iT ==1){
ll = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow = -2*ll + K*log(sum(vNj))
BIChigh = -2*ll + K*log(iN)
AIC = -2*ll + 2*K
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM > 1 & iT == 1){
ll = -Inf
BIClow = Inf
BIChigh = Inf
AIC = Inf
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM==1 & iT > 1){
LCAout = LCA_fast_init(mY,iT)
ll = tail(LCAout$LLKSeries,1)
npar = iT*iH + iT - 1
BIClow = LCAout$BIC
BIChigh = -2*ll + npar*log(iN)
AIC = -2*ll + 2*npar
ICL_BIClow = BIClow + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh = Inf
outmuLCA = LCAout
}else{
# Actual multilevel LC fit
start = meas_Init(mY,id_high,vNj,iM,iT)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA = MLTLCA(mY, vNj, vOmegast, mPi, start$mPhi_start)
ll = tail(outmuLCA$LLKSeries,1)
BIClow = outmuLCA$BIClow
BIChigh = outmuLCA$BIChigh
AIC = outmuLCA$AIC
ICL_BIClow = outmuLCA$ICL_BIClow
ICL_BIChigh = outmuLCA$ICL_BIChigh
}
return(list(modFIT = outmuLCA,ll=ll, BIClow = BIClow, BIChigh = BIChigh, AIC = AIC,
ICL_BIClow = ICL_BIClow, ICL_BIChigh = ICL_BIChigh))
}
#
lukosel_fun = function(mY,id_high,iT_range,iM_range,verbose){
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
# Step 1 - lower level
LCAout = list()
BIClow_step1 = rep(NA,num_iT)
BIChigh_step1 = rep(NA,num_iT)
AIC_step1 = rep(NA,num_iT)
ICL_BIClow_step1 = rep(NA,num_iT)
ICL_BIChigh_step1 = rep(NA,num_iT)
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting single-level measurement model (step 1)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow_step1[1] = -2*llfoo + K*log(sum(vNj))
BIChigh_step1[1] = -2*llfoo + K*log(iN)
AIC_step1[1] = -2*llfoo + 2*K
ICL_BIClow_step1[1] = Inf
ICL_BIChigh_step1[1] = Inf
} else{
LCAout = LCA_fast_init(mY,i)
ll = tail(LCAout$LLKSeries,1)
npar = i*iH + i - 1
BIClow_step1[1+i-iT_min] = LCAout$BIC
BIChigh_step1[1+i-iT_min] = -2*ll + npar*log(iN)
AIC_step1[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow_step1[1+i-iT_min] = BIClow_step1[i] + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh_step1[1+i-iT_min] = Inf
}
}
iT_currbest = which.min(BIClow_step1)+iT_min-1
if(iT_currbest==1) stop("iT=1 optimal.",call.=FALSE)
# Step 2 - higher level
outmuLCA2 = list()
BIClow_step2 = rep(NA,num_iM)
BIChigh_step2 = rep(NA,num_iM)
AIC_step2 = rep(NA,num_iM)
ICL_BIClow_step2 = rep(NA,num_iM)
ICL_BIChigh_step2 = rep(NA,num_iM)
if(iT_currbest > 1){
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 2)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
BIClow_step2[1] = BIClow_step1[1+iT_currbest-iT_min]
BIChigh_step2[1] = BIChigh_step1[1+iT_currbest-iT_min]
AIC_step2[1] = AIC_step1[1+iT_currbest-iT_min]
ICL_BIClow_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
ICL_BIChigh_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
} else{
start = meas_Init(mY,id_high,vNj,i,iT_currbest)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA2[[1+i-iM_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIClow
BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIChigh
AIC_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$AIC
ICL_BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_currbest = which.min(BIChigh_step2)+iM_min-1
} else{
iM_currbest = 1
}
# Step 3 - revisiting lower level
outmuLCA3 = list()
BIClow_step3 = rep(NA,num_iT)
BIChigh_step3 = rep(NA,num_iT)
AIC_step3 = rep(NA,num_iT)
ICL_BIClow_step3 = rep(NA,num_iT)
ICL_BIChigh_step3 = rep(NA,num_iT)
if(iM_currbest > 1){
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 3)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init(mY,id_high,vNj,iM_currbest,i)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA3[[1+i-iT_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIClow
BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIChigh
AIC_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$AIC
ICL_BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_currbest = which(BIClow_step3==min(BIClow_step3,na.rm=T))+iT_min-1
outmuLCA_step3 = outmuLCA3[[1+iT_currbest-iT_min]]
} else{
outmuLCA_step3 = LCA_fast_init(mY,iT_currbest)
}
return(list(outmuLCA_step3=outmuLCA_step3,iT_opt=iT_currbest, iM_opt=iM_currbest,
BIClow_step1 = BIClow_step1, BIChigh_step1 = BIChigh_step1, AIC_step1 = AIC_step1,
ICL_BIClow_step1 = ICL_BIClow_step1, ICL_BIChigh_step1 = ICL_BIChigh_step1,
BIClow_step2 = BIClow_step2, BIChigh_step2 = BIChigh_step2, AIC_step2 = AIC_step2,
ICL_BIClow_step2 = ICL_BIClow_step2, ICL_BIChigh_step2 = ICL_BIChigh_step2,
BIClow_step3 = BIClow_step3, BIChigh_step3 = BIChigh_step3, AIC_step3 = AIC_step3,
ICL_BIClow_step3 = ICL_BIClow_step3, ICL_BIChigh_step3 = ICL_BIChigh_step3))
}
#
sel_other = function(mY,id_high,iT_range,iM_range,approach,verbose){
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
if(!is.null(iM_range)){
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
if(num_iM==1){
if(iM_range==1){
iM_range=NULL
}
}
}
#
if(approach=="model selection on low"|approach=="model selection on low with high"){
BIClow = rep(NA,num_iT)
BIChigh = rep(NA,num_iT)
AIC = rep(NA,num_iT)
ICL_BIClow = rep(NA,num_iT)
ICL_BIChigh = rep(NA,num_iT)
if(is.null(iM_range)){
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = NA
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow[1] = -2*llfoo + K*log(sum(vNj))
if(is.null(iM_range)) BIClow[1] = -2*llfoo + K*log(nrow(mY))
if(!is.null(iM_range)) BIChigh[1] = -2*llfoo + K*log(iN)
AIC[1] = -2*llfoo + 2*K
ICL_BIClow[1] = Inf
ICL_BIChigh[1] = Inf
} else{
outmuLCA[[1+i-iT_min]] = LCA_fast_init(mY,i)
ll = tail(outmuLCA[[1+i-iT_min]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIC
BIChigh[1+i-iT_min] = -2*ll + npar*log(iN)
AIC[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow[1+i-iT_min] = BIClow[i] + 2*sum(outmuLCA[[1+i-iT_min]]$freq*apply(-outmuLCA[[1+i-iT_min]]$mU*log(outmuLCA[[1+i-iT_min]]$mU),1,sum))
ICL_BIChigh[1+i-iT_min] = Inf
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
} else{
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init(mY,id_high,vNj,iM_range,i)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iT_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIClow
BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIChigh
AIC[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$AIC
ICL_BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
}
if(iT_best==1) stop("iT=1 optimal.",call.=FALSE)
return(list(outmuLCA=outmuLCA,iT_best=iT_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
} else if(approach=="model selection on high"){
BIClow = rep(NA,num_iM)
BIChigh = rep(NA,num_iM)
AIC = rep(NA,num_iM)
ICL_BIClow = rep(NA,num_iM)
ICL_BIChigh = rep(NA,num_iM)
if(iT_range > 1){
outmuLCA = list()
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = LCA_fast_init(mY,iT_range)
ll = tail(outmuLCA[[1]]$LLKSeries,1)
npar = iT_range*iH + iT_range - 1
BIClow[1] = outmuLCA[[1]]$BIC
BIChigh[1] = -2*ll + npar*log(iN)
AIC[1] = -2*ll + 2*npar
ICL_BIClow[1] = BIClow[1] + 2*sum(outmuLCA[[1]]$freq*apply(-outmuLCA[[1]]$mU*log(outmuLCA[[1]]$mU),1,sum))
ICL_BIChigh[1] = Inf
} else{
start = meas_Init(mY,id_high,vNj,i,iT_range)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iM_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIClow
BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIChigh
AIC[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$AIC
ICL_BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_best = which.min(BIChigh)+iM_min-1
outmuLCA = outmuLCA[[1+iM_best-iM_min]]
} else{
stop("iM=1 optimal.",call.=FALSE)
}
return(list(outmuLCA=outmuLCA,iM_best=iM_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
}
}
#
simultsel_fun_poly = function(mY,id_high,iT,iM,ivItemcat){
LCAout = NULL
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
if(iM ==1 & iT ==1){
ll = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow = -2*ll + K*log(sum(vNj))
BIChigh = -2*ll + K*log(iN)
AIC = -2*ll + 2*K
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM > 1 & iT == 1){
ll = -Inf
BIClow = Inf
BIChigh = Inf
AIC = Inf
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM==1 & iT > 1){
LCAout = LCA_fast_init_poly(mY,iT,ivItemcat)
ll = tail(LCAout$LLKSeries,1)
npar = iT*iH + iT - 1
BIClow = LCAout$BIC
BIChigh = -2*ll + npar*log(iN)
AIC = -2*ll + 2*npar
ICL_BIClow = BIClow + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh = Inf
outmuLCA = LCAout
}else{
# Actual multilevel LC fit
start = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA = MLTLCA_poly(mY,vNj,vOmegast,mPi,
start$mPhi_start,ivItemcat)
ll = tail(outmuLCA$LLKSeries,1)
BIClow = outmuLCA$BIClow
BIChigh = outmuLCA$BIChigh
AIC = outmuLCA$AIC
ICL_BIClow = outmuLCA$ICL_BIClow
ICL_BIChigh = outmuLCA$ICL_BIChigh
}
return(list(modFIT = outmuLCA,ll=ll, BIClow = BIClow, BIChigh = BIChigh, AIC = AIC,
ICL_BIClow = ICL_BIClow, ICL_BIChigh = ICL_BIChigh))
}
#
lukosel_fun_poly = function(mY,id_high,iT_range,iM_range,ivItemcat,verbose){
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
# Step 1 - lower level
LCAout = list()
BIClow_step1 = rep(NA,num_iT)
BIChigh_step1 = rep(NA,num_iT)
AIC_step1 = rep(NA,num_iT)
ICL_BIClow_step1 = rep(NA,num_iT)
ICL_BIChigh_step1 = rep(NA,num_iT)
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting single-level measurement model (step 1)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow_step1[1] = -2*llfoo + K*log(sum(vNj))
BIChigh_step1[1] = -2*llfoo + K*log(iN)
AIC_step1[1] = -2*llfoo + 2*K
ICL_BIClow_step1[1] = Inf
ICL_BIChigh_step1[1] = Inf
} else{
LCAout = LCA_fast_init_poly(mY,i,ivItemcat)
ll = tail(LCAout$LLKSeries,1)
npar = i*iH + i - 1
BIClow_step1[1+i-iT_min] = LCAout$BIC
BIChigh_step1[1+i-iT_min] = -2*ll + npar*log(iN)
AIC_step1[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow_step1[1+i-iT_min] = BIClow_step1[i] + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh_step1[1+i-iT_min] = Inf
}
}
iT_currbest = which.min(BIClow_step1)+iT_min-1
if(iT_currbest==1) stop("iT=1 optimal.",call.=FALSE)
# Step 2 - higher level
outmuLCA2 = list()
BIClow_step2 = rep(NA,num_iM)
BIChigh_step2 = rep(NA,num_iM)
AIC_step2 = rep(NA,num_iM)
ICL_BIClow_step2 = rep(NA,num_iM)
ICL_BIChigh_step2 = rep(NA,num_iM)
if(iT_currbest > 1){
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 2)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
BIClow_step2[1] = BIClow_step1[1+iT_currbest-iT_min]
BIChigh_step2[1] = BIChigh_step1[1+iT_currbest-iT_min]
AIC_step2[1] = AIC_step1[1+iT_currbest-iT_min]
ICL_BIClow_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
ICL_BIChigh_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
} else{
start = meas_Init_poly(mY,id_high,vNj,i,iT_currbest,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA2[[1+i-iM_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIClow
BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIChigh
AIC_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$AIC
ICL_BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_currbest = which.min(BIChigh_step2)+iM_min-1
} else{
iM_currbest = 1
}
# Step 3 - revisiting lower level
outmuLCA3 = list()
BIClow_step3 = rep(NA,num_iT)
BIChigh_step3 = rep(NA,num_iT)
AIC_step3 = rep(NA,num_iT)
ICL_BIClow_step3 = rep(NA,num_iT)
ICL_BIChigh_step3 = rep(NA,num_iT)
if(iM_currbest > 1){
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 3)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>2){
start = meas_Init_poly(mY,id_high,vNj,iM_currbest,i,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA3[[1+i-iT_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIClow
BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIChigh
AIC_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$AIC
ICL_BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_currbest = which(BIClow_step3==min(BIClow_step3,na.rm=T))+iT_min-1
outmuLCA_step3 = outmuLCA3[[1+iT_currbest-iT_min]]
} else{
outmuLCA_step3 = LCAout
}
return(list(outmuLCA_step3=outmuLCA_step3,iT_opt=iT_currbest, iM_opt=iM_currbest,
BIClow_step1 = BIClow_step1, BIChigh_step1 = BIChigh_step1, AIC_step1 = AIC_step1,
ICL_BIClow_step1 = ICL_BIClow_step1, ICL_BIChigh_step1 = ICL_BIChigh_step1,
BIClow_step2 = BIClow_step2, BIChigh_step2 = BIChigh_step2, AIC_step2 = AIC_step2,
ICL_BIClow_step2 = ICL_BIClow_step2, ICL_BIChigh_step2 = ICL_BIChigh_step2,
BIClow_step3 = BIClow_step3, BIChigh_step3 = BIChigh_step3, AIC_step3 = AIC_step3,
ICL_BIClow_step3 = ICL_BIClow_step3, ICL_BIChigh_step3 = ICL_BIChigh_step3))
}
#
sel_other_poly = function(mY,id_high,iT_range,iM_range,ivItemcat,approach,verbose){
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
if(!is.null(iM_range)){
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
if(num_iM==1){
if(iM_range==1){
iM_range=NULL
}
}
}
#
if(approach=="model selection on low"|approach=="model selection on low with high"){
BIClow = rep(NA,num_iT)
BIChigh = rep(NA,num_iT)
AIC = rep(NA,num_iT)
ICL_BIClow = rep(NA,num_iT)
ICL_BIChigh = rep(NA,num_iT)
if(is.null(iM_range)){
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = NA
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow[1] = -2*llfoo + K*log(sum(vNj))
if(is.null(iM_range)) BIClow[1] = -2*llfoo + K*log(nrow(mY))
if(!is.null(iM_range)) BIChigh[1] = -2*llfoo + K*log(iN)
AIC[1] = -2*llfoo + 2*K
ICL_BIClow[1] = Inf
ICL_BIChigh[1] = Inf
} else{
outmuLCA[[1+i-iT_min]] = LCA_fast_init_poly(mY,i,ivItemcat)
ll = tail(outmuLCA[[1+i-iT_min]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIC
BIChigh[1+i-iT_min] = -2*ll + npar*log(iN)
AIC[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow[1+i-iT_min] = BIClow[i] + 2*sum(outmuLCA[[1+i-iT_min]]$freq*apply(-outmuLCA[[1+i-iT_min]]$mU*log(outmuLCA[[1+i-iT_min]]$mU),1,sum))
ICL_BIChigh[1+i-iT_min] = Inf
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
} else{
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init_poly(mY,id_high,vNj,iM_range,i,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iT_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIClow
BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIChigh
AIC[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$AIC
ICL_BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
}
if(iT_best==1) stop("iT=1 optimal.",call.=FALSE)
return(list(outmuLCA=outmuLCA,iT_best=iT_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
} else if(approach=="model selection on high"){
BIClow = rep(NA,num_iM)
BIChigh = rep(NA,num_iM)
AIC = rep(NA,num_iM)
ICL_BIClow = rep(NA,num_iM)
ICL_BIChigh = rep(NA,num_iM)
if(iT_range > 1){
outmuLCA = list()
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = LCA_fast_init_poly(mY,iT_range,ivItemcat)
ll = tail(outmuLCA[[1]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1] = outmuLCA[[1]]$BIC
BIChigh[1] = -2*ll + npar*log(iN)
AIC[1] = -2*ll + 2*npar
ICL_BIClow[1] = BIClow[i] + 2*sum(outmuLCA[[1]]$freq*apply(-outmuLCA[[1]]$mU*log(outmuLCA[[1]]$mU),1,sum))
ICL_BIChigh[1] = Inf
} else{
start = meas_Init_poly(mY,id_high,vNj,i,iT_range,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iM_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIClow
BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIChigh
AIC[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$AIC
ICL_BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_best = which.min(BIChigh)+iM_min-1
outmuLCA = outmuLCA[[1+iM_best-iM_min]]
} else{
stop("iM=1 optimal.",call.=FALSE)
}
return(list(outmuLCA=outmuLCA,iM_best=iM_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
}
}
#
#
#
clean_output1 = function(output,Y,iT,iH,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add vector of class proportions
out$vPi = output$pg
rownames(out$vPi) = paste0("P(C",1:iT,")")
colnames(out$vPi) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add vector of class proportions
out$vPi = output$pg
rownames(out$vPi) = paste0("P(C",1:iT,")")
colnames(out$vPi) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add matrix of posterior class membership probabilities
out$mU = output$mU[rep(1:dim(output$mU)[1],output$freq),]
colnames(out$mU) = colnames(out$mPhi)
if(dataout){
out$mU = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$mU)
rownames(out$mU) = NULL
}
# Add matrix of modal class assignment
out$mU_modal = output$mModalAssnm[rep(1:dim(output$mModalAssnm)[1],output$freq),]
colnames(out$mU_modal) = colnames(out$mPhi)
if(dataout){
out$mU_modal = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$mU_modal)
rownames(out$mU_modal) = NULL
}
# Add vector of modal class assignment
out$vU_modal = output$vModalAssnm[rep(1:dim(output$vModalAssnm)[1],output$freq),,drop=FALSE]+1
colnames(out$vU_modal) = "C"
if(dataout){
out$vU_modal = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$vU_modal)
rownames(out$vU_modal) = NULL
}
# Add matrix of classification errors
out$mClassErr = output$mClassErr
rownames(out$mClassErr) = paste0("C",1:iT,"_true")
colnames(out$mClassErr) = paste0("C",1:iT,"_pred")
# Add matrix of proportion of classification errors
out$mClassErrProb = output$mClassErrProb
rownames(out$mClassErrProb) = rownames(out$mClassErr)
colnames(out$mClassErrProb) = colnames(out$mClassErr)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$vGamma = output$alphas
rownames(out$vGamma) = paste0("gamma(C",2:iT,")")
colnames(out$vGamma) = ""
# Add betas
out$mBeta = output$gamma
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vGamma),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of standard errors
out$SEs = output$SEs
rownames(out$SEs) = rownames(out$parvec)
colnames(out$SEs) = ""
# Add variance-covariance matrix
out$Varmat = output$Varmat
rownames(out$Varmat) = colnames(out$Varmat) = rownames(out$parvec)
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add specification
out$spec = as.matrix("Single-level LC model")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output2 = function(output,Y,iT,Z,iH,P,extout,dataout,ivItemcat){
mY = output$mY
mZ = output$mZ
output = output$outcov
if(!extout){
# Create empty list for output
out = list()
# Add vector of sample means of class proportions
out$vPi_avg = as.matrix(apply(output$mPg,2,mean))
rownames(out$vPi_avg) = paste0("P(C",1:iT,")")
colnames(out$vPi_avg) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = output$beta
rownames(out$cGamma) = paste0("gamma(",Z,"|C)")
colnames(out$cGamma) = paste0("C",2:iT)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = matrix(output$SEs_cor[1:((iT-1)*P),],P,iT-1)
rownames(out$SEs_cor_gamma) = rownames(out$cGamma)
colnames(out$SEs_cor_gamma) = colnames(out$cGamma)
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model with covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of class proportions
out$mPi = output$mPg
colnames(out$mPi) = paste0("P(C",1:iT,")")
# Add vector of sample means of class proportions
out$vPi_avg = as.matrix(apply(output$mPg,2,mean))
rownames(out$vPi_avg) = colnames(out$mPi)
colnames(out$vPi_avg) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add matrix of posterior class membership probabilities
out$mU = output$mU
colnames(out$mU) = colnames(out$mPhi)
if(dataout){
data = cbind(mY,mZ)
rownames(data) = NULL
colnames(data) = c(Y,Z)
data = data[,-which(colnames(data)=="Intercept"),drop=FALSE]
out$mU = cbind(data,out$mU)
rownames(out$mU) = NULL
}
# Add matrix of classification errors
out$mClassErr = output$mClassErr
rownames(out$mClassErr) = paste0("C",1:iT,"_true")
colnames(out$mClassErr) = paste0("C",1:iT,"_pred")
# Add matrix of proportion of classification errors
out$mClassErrProb = output$mClassErrProb
rownames(out$mClassErrProb) = paste0("C",1:iT,"_true")
colnames(out$mClassErrProb) = paste0("C",1:iT,"_pred")
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = output$beta
rownames(out$cGamma) = paste0("gamma(",Z,"|C)")
colnames(out$cGamma) = paste0("C",2:iT)
# Add betas
out$mBeta = output$gamma
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(paste0(rep(substr(rownames(out$cGamma),1,nchar(rownames(out$cGamma))-2),iT-1),rep(colnames(out$cGamma),rep(P,iT-1)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:((iT-1)*P),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:((iT-1)*P),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for gamma
out$SEs_cor_gamma = matrix(output$SEs_cor[1:((iT-1)*P),],P,iT-1)
rownames(out$SEs_cor_gamma) = rownames(out$cGamma)
colnames(out$SEs_cor_gamma) = colnames(out$cGamma)
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mQ) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat_unc[1:((iT-1)*P),1:((iT-1)*P)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$Varmat_cor
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$mV2)
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = paste0(rep(substr(rownames(out$cGamma),1,nchar(rownames(out$cGamma))-2),iT-1),rep(colnames(out$cGamma),rep(P,iT-1)),")")
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model with covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output3 = function(output,Y,iT,iM,mY,id_high,iH,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = output$mPi
rownames(out$mPi) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = output$mPi
rownames(out$mPi) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi)))
if(dataout){
data = cbind(mY,id_high)
rownames(data) = NULL
colnames(data) = c(Y,id_high_name)
for(i in 1:length(id_high_levs)){
data[data[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX,dim(output$cLogPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX,dim(output$cPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi)
if(dataout){
out$mSumPX = cbind(data,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi)
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi)
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi)
if(dataout){
out$mPW_N = cbind(data,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$vAlpha = output$vDelta
rownames(out$vAlpha) = paste0("alpha(G",2:iM,")")
colnames(out$vAlpha) = ""
# Add gammas
out$mGamma = output$mGamma
rownames(out$mGamma) = paste0("gamma(C",2:iT,"|G)")
colnames(out$mGamma) = colnames(out$mPi)
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vAlpha),paste0(rep(substr(rownames(out$mGamma),1,nchar(rownames(out$mGamma))-2),iM),rep(colnames(out$mGamma),rep(iT-1,iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of standard errors
out$SEs = output$SEs
rownames(out$SEs) = rownames(out$parvec)
colnames(out$SEs) = ""
# Add variance-covariance matrix
out$Varmat = output$Varmat
rownames(out$Varmat) = colnames(out$Varmat) = rownames(out$parvec)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = colnames(out$Infomat) = rownames(out$parvec)
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add specification
out$spec = as.matrix("Multilevel LC model")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output4 = function(output,Y,iT,iM,Z,mY,mZ,id_high,iH,P,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add vector of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[iM:(iM-1+(P*(iT-1)*iM))],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model with lower-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = matrix(output$cPi,nrow(output$cPi),iT*iM)
colnames(out$mPi) = paste0(rep(paste0("P(C",1:iT,"|"), iM),rep(paste0("G",1:iM,")"),rep(iT,iM)))
# Add vector of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi_avg)))
if(dataout){
data = cbind(mY,id_high,mZ)
rownames(data) = NULL
colnames(data) = c(Y,id_high_name,Z)
for(i in 1:length(id_high_levs)){
data[data[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
data = data[,-which(colnames(data)=="Intercept"),drop=FALSE]
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX, dim(output$cLogPMX), dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi_avg)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX, dim(output$cPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi_avg)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi_avg)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi_avg)
if(dataout){
out$mSumPX = cbind(data,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi_avg)
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi_avg)
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi_avg)
if(dataout){
out$mPW_N = cbind(data,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$vAlpha = output$vDelta
rownames(out$vAlpha) = paste0("alpha(G",2:iM,")")
colnames(out$vAlpha) = ""
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),paste0("G",1:iM)))
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vAlpha),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:(iM-1+P*(iT-1)*iM),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:(iM-1+P*(iT-1)*iM),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[iM:(iM-1+(P*(iT-1)*iM))],c(iT-1,P,iM)),3,t),dim(out$cGamma), dimnames=dimnames(out$cGamma))
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mQ) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat[1:(iM-1+P*(iT-1)*iM),1:(iM-1+P*(iT-1)*iM)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$mVar_corr
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$mV2)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = rownames(out$parvec)
# Add cube of Fisher information matrix for gamma
out$cGamma_Info = array(output$cGamma_Info,dim(output$cGamma_Info),dimnames=list(paste0("gamma(",Z,"|C,G)"),paste0("gamma(",Z,"|C,G)"),paste0("C", 2:iT,rep(paste0(",G", 1:iM),rep(iT-1,iM)))))
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = c(rownames(out$vAlpha),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"))
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add subset of matrix of model parameter contributions to log-likelihood score for gamma
out$mGamma_Score = output$mGamma_Score
colnames(out$mGamma_Score) = paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")")
# Add specification
out$spec = as.matrix("Multilevel LC model with lower-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output5 = function(output,Y,iT,iM,Z,Zh,mY,mZ,mZh,id_high,iH,P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of sample means of high-level class proportions
out$vOmega_avg = as.matrix(apply(output$mOmega,2,mean))
rownames(out$vOmega_avg) = paste0("P(G",1:iM,")")
colnames(out$vOmega_avg) = ""
# Add matrix of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$mAlpha = t(output$mDelta)
rownames(out$mAlpha) = paste0("alpha(",Zh,"|G)")
colnames(out$mAlpha) = paste0("G",2:iM)
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add corrected standard errors for alpha
out$SEs_cor_alpha = t(matrix(output$SEs_cor[1:(P_high*(iM-1))],iM-1,P_high))
rownames(out$SEs_cor_alpha) = rownames(out$mAlpha)
colnames(out$SEs_cor_alpha) = colnames(out$mAlpha)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[(1+P_high*(iM-1)):(P_high*(iM-1)+P*(iT-1)*iM)],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model with lower- and higher-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$mOmega = output$mOmega
colnames(out$mOmega) = paste0("P(G",1:iM,")")
# Add matrix of sample means of high-level class proportions
out$vOmega_avg = as.matrix(apply(output$mOmega,2,mean))
rownames(out$vOmega_avg) = colnames(out$mOmega)
colnames(out$vOmega_avg) = ""
# Add matrix of conditional low-level class proportions
out$mPi = matrix(output$cPi,nrow(output$cPi),iT*iM)
colnames(out$mPi) = paste0(rep(paste0("P(C",1:iT,"|"),iM),rep(paste0("G",1:iM,")"),rep(iT,iM)))
# Add matrix of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),paste0("G", 1:iM)))
if(dataout){
data_low = cbind(mY,id_high,mZ)
rownames(data_low) = NULL
colnames(data_low) = c(Y,id_high_name,Z)
for(i in 1:length(id_high_levs)){
data_low[data_low[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
data_low = data_low[,-which(colnames(data_low)=="Intercept"),drop=FALSE]
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data_low,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data_low)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX,dim(output$cLogPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),paste0("G", 1:iM)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data_low,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data_low)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX,dim(output$cPX),dimnames=list(NULL,paste0("C", 1:iT,"|G"),paste0("G", 1:iM)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data_low,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data_low)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C", 1:iT,"|G"),paste0("G", 1:iM)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data_low,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data_low)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi_avg)
if(dataout){
out$mSumPX = cbind(data_low,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi_avg)
if(dataout){
data_high = mZh
rownames(data_high) = id_high_levs
colnames(data_high) = Zh
data_high = data_high[,-which(colnames(data_high)=="Intercept"),drop=FALSE]
out$mPW = cbind(data_high,out$mPW)
}
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi_avg)
if(dataout){
out$mlogPW = cbind(data_high,out$mlogPW)
}
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi_avg)
if(dataout){
out$mPW_N = cbind(data_low,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data_low,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$mAlpha = t(output$mDelta)
rownames(out$mAlpha) = paste0("alpha(",Zh,"|G)")
colnames(out$mAlpha) = paste0("G",2:iM)
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")"),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:(P_high*(iM-1)+P*(iT-1)*iM),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:(P_high*(iM-1)+P*(iT-1)*iM),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for alpha
out$SEs_cor_alpha = t(matrix(output$SEs_cor[1:(P_high*(iM-1))],iM-1,P_high))
rownames(out$SEs_cor_alpha) = rownames(out$mAlpha)
colnames(out$SEs_cor_alpha) = colnames(out$mAlpha)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[(1+P_high*(iM-1)):(P_high*(iM-1)+P*(iT-1)*iM)],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mV2) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat[1:(P_high*(iM-1)+P*(iT-1)*iM),1:(P_high*(iM-1)+P*(iT-1)*iM)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$mVar_corr
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$Varmat_unc)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = colnames(out$Infomat) = rownames(out$parvec)
# Add cube of Fisher information matrix for alpha
out$cAlpha_Info = array(output$cDelta_Info, dim(output$cDelta_Info), dimnames=list(paste0("alpha(",Zh,"|G)"),paste0("alpha(",Zh,"|G)"),colnames(out$mAlpha)))
# Add cube of Fisher information matrix for gamma
out$cGamma_Info = array(output$cGamma_Info, dim(output$cGamma_Info), dimnames=list(paste0("gamma(",Z,"|C,G)"),paste0("gamma(",Z,"|C,G)"),paste0("C",2:iT,rep(paste0(",G",1:iM),rep(iT-1,iM)))))
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = c(paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")"),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"))
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add subset of matrix of model parameter contributions to log-likelihood score for alpha
out$mAlpha_Score = output$mDelta_Score
rownames(out$mAlpha_Score) = id_high_levs
colnames(out$mAlpha_Score) = paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")")
# Add subset of matrix of model parameter contributions to log-likelihood score for gamma
out$mGamma_Score = output$mGamma_Score
colnames(out$mGamma_Score) = paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")")
# Add specification
out$spec = as.matrix("Multilevel LC model with lower- and higher-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
#
clean_cov = function(data,covnam){
cov = data[,covnam,drop=FALSE]
if(any(!sapply(cov,is.numeric))){
covnam_upd = c()
for(i in covnam){
if(!sapply(cov,is.numeric)[which(colnames(cov)==i)]){
covnam_lev = sort(unique(cov[,i]))
cov_missing = is.na(as.numeric(factor(cov[,i])))
cov_upd = vecTomatClass(as.numeric(factor(cov[,i])))
cov_upd[cov_missing,] = NA
colnames(cov_upd) = paste0(i,".",covnam_lev)
cov_upd = cov_upd[,-1,drop=FALSE]
covnam_upd = c(covnam_upd,colnames(cov_upd))
} else{
cov_upd = cov[,i,drop=FALSE]
covnam_upd = c(covnam_upd,i)
}
cov = cov[,-which(colnames(cov)==i),drop=FALSE]
cov = cbind(cov,cov_upd)
}
covnam = covnam_upd
colnames(cov) = covnam
}
cov = as.matrix(cbind(1,cov))
return(list(cov=cov,covnam=covnam))
}
#
update_YmY = function(mY,Y,ivItemcat){
Y_upd = c()
for(i in Y){
if(ivItemcat[i]>2){
mY_upd = vecTomatClass(mY[,i]+1)
colnames(mY_upd) = paste0(i,".",1:ncol(mY_upd)-1)
Y_upd = c(Y_upd,colnames(mY_upd))
} else{
mY_upd = mY[,i,drop=FALSE]
Y_upd = c(Y_upd,i)
}
mY = mY[,-which(colnames(mY)==i)]
mY = cbind(mY,mY_upd)
}
Y = Y_upd
colnames(mY) = Y
return(list(mY=mY,Y=Y))
}
#
update_YmY_nonnum = function(mY,Y,ivItemcat){
itemnam = lapply(data.frame(mY),function(x){levels(factor(x))})
mY = apply(mY,2,function(x){as.numeric(factor(x))-1})
Y_upd = c()
for(i in Y){
if(ivItemcat[i]>2){
mY_upd = vecTomatClass(mY[,i]+1)
colnames(mY_upd) = paste0(i,".",itemnam[[i]])
Y_upd = c(Y_upd,colnames(mY_upd))
} else{
mY_upd = mY[,i,drop=FALSE]
Y_upd = c(Y_upd,paste0(i,".",itemnam[[i]][2]))
}
mY = mY[,-which(colnames(mY)==i)]
mY = cbind(mY,mY_upd)
}
Y = Y_upd
colnames(mY) = Y
return(list(mY=mY,Y=Y))
}
#
check_inputs1 = function(data,Y,iT,id_high,iM,Z,Zh){
# <data> is matrix or dataframe
if(!(is.matrix(data)|is.data.frame(data))){
stop("data must be matrix or dataframe.",call.=FALSE)
}
# <Y>, <id_high>, <Z> and <Zh> are column names in <data>
if(any(!c(Y,id_high,Z,Zh)%in%colnames(data))){
if(any(!Y%in%colnames(data))){
stop("Not all Y in data.",call.=FALSE)
}
if(!is.null(id_high)&any(!id_high%in%colnames(data))){
stop("id_high not in data.",call.=FALSE)
}
if(!is.null(Z)&any(!Z%in%colnames(data))){
stop("Not all Z in data.",call.=FALSE)
}
if(!is.null(Zh)&any(!Zh%in%colnames(data))){
stop("Not all Zh in data.",call.=FALSE)
}
}
}
#
check_inputs2 = function(data,Y,iT,id_high,iM,Z,Zh){
# Single <id_high>
if(!is.null(id_high)){
if(length(id_high)!=1){
stop("Invalid id_high.",call.=FALSE)
}
}
# <iT> positive integer or sequential positive integers
if(!length(iT)>=1){
stop("Invalid iT.",call.=FALSE)
} else{
if(length(iT)==1){
if(!is.numeric(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(iT!=round(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(iT==1){
stop("Invalid iT: not a mixture model.",call.=FALSE)
} else if(iT<1){
stop("Invalid iT.",call.=FALSE)
}
} else{
if(!is.vector(iT)|is.list(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(!is.numeric(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(any(!iT==round(iT))){
stop("Invalid iT.",call.=FALSE)
} else if(min(iT)<1){
stop("Invalid iT.",call.=FALSE)
} else if(length(iT)!=length(min(iT):max(iT))){
stop("Invalid iT.",call.=FALSE)
} else if(any(!iT==min(iT):max(iT))){
stop("Invalid iT.",call.=FALSE)
}
}
}
# <iM> positive integer or sequential positive integers
if(!is.null(iM)&!length(iM)>=1){
stop("Invalid iM.",call.=FALSE)
} else if(!is.null(iM)){
if(length(iM)==1){
if(!is.numeric(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(iM!=round(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(iM<1){
stop("Invalid iM.",call.=FALSE)
}
} else{
if(!is.vector(iM)|is.list(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(!is.numeric(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(any(!iM==round(iM))){
stop("Invalid iM.",call.=FALSE)
} else if(min(iM)<1){
stop("Invalid iM.",call.=FALSE)
} else if(length(iM)!=length(min(iM):max(iM))){
stop("Invalid iM.",call.=FALSE)
} else if(any(!iM==min(iM):max(iM))){
stop("Invalid iM.",call.=FALSE)
}
}
}
# <Y> consecutive integers from 0 or be non-numeric character
if(any(sapply(as.data.frame(data[,Y]),function(x){(any(sort(unique(x))!=0:(length(unique(x))-1)))&!is.character(x)}))){
stop("Items must contain consecutive integers from 0 or be non-numeric character.",call.=FALSE)
}
# No missings in <id_high>
if(any(is.na(data[,id_high]))){
stop("id_high cannot contain missing values.",call.=FALSE)
}
# <id_high> not labelled "Intercept"
if(!is.null(id_high)){
if(id_high=="Intercept"){
stop("id_high may not be labelled 'Intercept'.",call.=FALSE)
}
}
# No duplicate <Y>, <Z> or >Zh>
if(any(duplicated(as.list(as.data.frame(data[,Y]))))){
stop("Duplicate items.",call.=FALSE)
}
if(!is.null(Z)){
if(length(Z)>2){
if(any(duplicated(as.list(data[,Z])))){
stop("Duplicate lower-level covariates.",call.=FALSE)
}
}
}
if(!is.null(Zh)){
if(length(Zh)>2){
if(any(duplicated(as.list(data[,Zh])))){
stop("Duplicate higher-level covariates.",call.=FALSE)
}
}
}
# No constant <Y>, <Z> or <Zh>
if(any(apply(data[,Y],2,function(x){length(unique(x))==1}))){
stop("Constant in indicators.",call.=FALSE)
}
if(!is.null(Z)){
if(any(apply(data[,Z,drop=FALSE],2,function(x){length(unique(x))==1}))){
stop("Constant in lower-level covariates.",call.=FALSE)
}
}
if(!is.null(Zh)){
if(any(apply(data[,Zh,drop=FALSE],2,function(x){length(unique(x))==1}))){
stop("Constant in higher-level covariates.",call.=FALSE)
}
}
# Identify specification
if(is.null(id_high)&is.null(iM)){
specification = "single-level"
} else if(!is.null(id_high)&!is.null(iM)){
specification = "multilevel"
} else{
if(is.null(id_high)){
stop("iM specified, id_high not.",call.=FALSE)
}
if(is.null(iM)){
stop("id_high specified, iM not.",call.=FALSE)
}
}
if(is.null(Z)){
specification = paste(specification,"measurement")
if(!is.null(Zh)){
stop("Zh specified, Z not.",call.=FALSE)
}
} else if(!is.null(Z)){
specification = paste(specification,"structural")
if(specification=="single-level structural"&!is.null(Zh)){
stop("Zh specified, iM and id_high not.",call.=FALSE)
}
}
# Identify approach
if(length(iT)==1&is.null(iM)){
approach = "direct"
} else if(length(iT)==1&length(iM)==1){
approach = "direct"
} else if(length(iT)>1&is.null(iM)){
approach = "model selection on low"
} else if(length(iT)>1&length(iM)==1){
if(iM==1){
approach = "model selection on low"
} else{
approach = "model selection on low with high"
}
} else if(length(iT)==1&length(iM)>1){
approach = "model selection on high"
} else if(length(iT)>1&length(iM)>1){
approach = "model selection on low and high"
}
return(approach)
}
#
#
#
print.multiLCA = function(x,...){
out = NULL
if(x$spec == "Single-level LC model"){
######################################
### Single-level measurement model ###
######################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nCLASS PROPORTIONS:\n")
print(round(x$vPi, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$AvgClassErrProb, 4)), as.character(round(x$R2entr, 4)), as.character(round(x$BIC, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("ClassErr", "EntR-sqr", "BIC", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
} else if(x$spec == "Single-level LC model with covariates"){
#####################################
### Single-level structural model ###
#####################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n")
print(round(x$vPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$AvgClassErrProb, 4)), as.character(round(x$R2entr, 4)), as.character(round(x$BIC, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("ClassErr", "EntR-sqr", "BIC", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR CLASS MEMBERSHIP:\n\n")
gammas = format(round(as.matrix(x$cGamma), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma)/as.matrix(x$SEs_cor_gamma), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma)/as.matrix(x$SEs_cor_gamma)))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (i in 1:ncol(gammas)){
C = paste0("C", 1 + i)
logit_params = noquote(cbind(gammas[,i], SE[,i], Zscore[,i], matrix(paste0(pval[,i], psignf[,i]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma)), 1, nchar(rownames(as.matrix(x$cGamma))) - 1), 1 + i, ")")
cat("\nMODEL FOR", C, "(BASE C1)\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
} else if(x$spec == "Multilevel LC model"){
####################################
### Multilevel measurement model ###
####################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS:\n")
print(round(x$vOmega, 4))
cat("\nCLASS PROPORTIONS:\n")
print(round(x$mPi, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
} else if(x$spec == "Multilevel LC model with lower-level covariates"){
###############################################################
### Multilevel structural model with lower-level covariates ###
###############################################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS:\n")
print(round(x$vOmega, 4))
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n\n")
print(round(x$mPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR LOWER-LEVEL CLASS MEMBERSHIP:\n\n")
for (i in 1:dim(x$cGamma)[3]){
gammas = format(round(as.matrix(x$cGamma[,,i]), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i])))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (j in 1:ncol(gammas)){
G = paste0("G", i)
C = paste0("C", 1 + j)
logit_params = noquote(cbind(gammas[,j], SE[,j], Zscore[,j], matrix(paste0(pval[,j], psignf[,j]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma[,,i])), 1, nchar(rownames(as.matrix(x$cGamma[,,i]))) - 1), j + 1, ",G", i, ")")
cat("\nMODEL FOR", C, "(BASE C1) GIVEN", G, "\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
}
} else if(x$spec == "Multilevel LC model with lower- and higher-level covariates"){
#########################################################################
### Multilevel structural model with low- and higher-level covariates ###
#########################################################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS (SAMPLE MEAN):\n")
print(round(x$vOmega_avg, 4))
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n\n")
print(round(x$mPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR HIGHER-LEVEL CLASS MEMBERSHIP:\n\n")
alphas = format(round(as.matrix(x$mAlpha), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_alpha), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$mAlpha)/as.matrix(x$SEs_cor_alpha), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$mAlpha)/as.matrix(x$SEs_cor_alpha)))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(alphas), ncol(alphas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (i in 1:ncol(alphas)){
G = paste0("G", 1 + i)
logit_params = noquote(cbind(alphas[,i], SE[,i], Zscore[,i], matrix(paste0(pval[,i], psignf[,i]))))
colnames(logit_params) = c("Alpha", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$mAlpha)), 1, nchar(rownames(as.matrix(x$mAlpha))) - 1), 1 + i, ")")
cat("\nMODEL FOR", G, "(BASE G1)\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR LOWER-LEVEL CLASS MEMBERSHIP:\n\n")
for (i in 1:dim(x$cGamma)[3]){
gammas = format(round(as.matrix(x$cGamma[,,i]), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i])))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (j in 1:ncol(gammas)){
G = paste0("G", i)
C = paste0("C", 1 + j)
logit_params = noquote(cbind(gammas[,j], SE[,j], Zscore[,j], matrix(paste0(pval[,j], psignf[,j]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma[,,i])), 1, nchar(rownames(as.matrix(x$cGamma[,,i]))) - 1), j + 1, ",G", i, ")")
cat("\nMODEL FOR", C, "(BASE C1) GIVEN", G, "\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
}
}
}
#
plot.multiLCA = function(x, horiz = FALSE, clab = NULL, ...){
out = NULL
iT = ncol(x$mPhi)
items = 1:nrow(x$mPhi)
itemnames = substr(rownames(x$mPhi), 3, nchar(rownames(x$mPhi)) - 3)
if(horiz == TRUE){
las = 1
} else{
las = 3
}
if(is.null(clab)){
legend = paste0("Class ", 1:iT)
} else{
if(!is.vector(clab) | is.list(clab)){
stop("Invalid clab.",call.=FALSE)
} else if(length(clab) != iT | !is.character(clab)){
stop("Invalid clab.",call.=FALSE)
}
legend = clab
}
oldpar = par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar = c(5, 5, 3, 8), xpd = TRUE)
plot(items, x$mPhi[,1], type = "b", xaxt = "n",yaxt = "n", pch = 1,
xlab = "", ylab = "Response probability", frame.plot = FALSE, ylim = c(0,1), lty = 1, ...)
axis(1, at = items, labels = itemnames, las = las)
axis(2, las = 1)
for(h in 2:iT){
lines(items, x$mPhi[,h], pch = h, type = "b", lty = h)
}
legend("topright", legend = legend, inset = c(-0.4,0),
lty = (1:iT), pch = 1:iT, cex = 0.8)
}
Try the multilevLCA package in your browser
Any scripts or data that you put into this service are public.
multilevLCA documentation built on Oct. 13, 2024, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.multiLCA print.multiLCA check_inputs2 check_inputs1 update_YmY_nonnum update_YmY clean_cov clean_output5 clean_output4 clean_output3 clean_output2 clean_output1 sel_other_poly lukosel_fun_poly simultsel_fun_poly sel_other lukosel_fun simultsel_fun meas_Init_poly LCA_fast_init_whigh_poly LCA_fast_init_wcov_poly LCA_fast_init_poly meas_Init LCA_fast_init_whigh LCA_fast_init_wcov LCA_fast_init multiLCA .onAttach
Documented in multiLCA plot.multiLCA
.onAttach <- function(libname, pkgname){
packageStartupMessage("multilevLCA: Estimates and plots single- and multilevel latent class models.")
}
multiLCA = function(data, Y, iT, id_high = NULL, iM = NULL, Z = NULL, Zh = NULL,
extout = FALSE, dataout = TRUE, kmea = TRUE, sequential = TRUE,
numFreeCores = 2, maxIter = 1e3, tol = 1e-8,
reord = TRUE, fixedpars = 1,
NRmaxit = 100, NRtol = 1e-6, verbose = TRUE){
#
if(verbose){
pb=txtProgressBar(char="Cleaning data...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
check_inputs1(data,Y,iT,id_high,iM,Z,Zh)
nrow_data = nrow(data)
data = data[complete.cases(data[,Y]),c(Y,id_high,Z,Zh)]
approach = check_inputs2(data,Y,iT,id_high,iM,Z,Zh)
reord = as.numeric(reord)
fixed = fixedpars
if((fixed==2)&is.null(iM)) fixed=fixedpars=1
if(!is.null(id_high)){
id_high_levs = sort(unique(data[,id_high]))
data[,id_high] = as.numeric(factor(data[,id_high]))
data = data[order(data[,id_high]),]
id_high_name = id_high
id_high = data[,id_high]
vNj = table(id_high)
}
mY = as.matrix(data[,Y])
ivItemcat = apply(mY,2,function(x){length(unique(x))})
itemchar = any(apply(data[,Y],2,function(x){!is.numeric(x)}))
if(any(ivItemcat>2)&!itemchar){
mY = update_YmY(mY,Y,ivItemcat)
Y = mY$Y
mY = mY$mY
} else if(itemchar){
mY = update_YmY_nonnum(mY,Y,ivItemcat)
Y = mY$Y
mY = mY$mY
}
if(!is.null(Z)){
mZ = clean_cov(data,Z)
Z = mZ$covnam
mZ = mZ$cov
}
if(!is.null(Zh)){
mZh = clean_cov(data,Zh)
Zh = mZh$covnam
mZh = mZh$cov
}
if(approach=="direct"){
#
if(verbose){
if(fixed==0){
pb=txtProgressBar(char="Fitting LC model...",width=1)
} else{
pb=txtProgressBar(char="Fitting measurement model...",width=1)
}
setTxtProgressBar(pb,1)
close(pb)
}
#
if(any(ivItemcat>2)){
#
if(is.null(id_high)&is.null(Z)&is.null(Zh)){
out = LCA_fast_init_poly(mY,iT,ivItemcat,kmea,maxIter,tol,reord)
out = clean_output1(out,Y,iT,length(Y),extout,dataout,ivItemcat)
} else if(is.null(id_high)&!is.null(Z)&is.null(Zh)){
out = LCA_fast_init_wcov_poly(mY,mZ,iT,ivItemcat,kmea,maxIter,tol,fixed,reord,
NRtol,NRmaxit,verbose)
out = clean_output2(out,Y,iT,c("Intercept",Z),length(Y),length(Z)+1,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&is.null(Z)&is.null(Zh)){
init = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
out = MLTLCA_poly(mY,vNj,init$vOmega_start,init$mPi_start,init$mPhi_start,
ivItemcat,maxIter,tol,reord)
out = clean_output3(out,Y,iT,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&is.null(Zh)){
init1 = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
init2 = MLTLCA_poly(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
vOmega_start = init2$vOmega
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
#
mY = mY[complete.cases(mZ),]
id_high = id_high[complete.cases(mZ)]
vNj = table(id_high)
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_cov_poly(mY,mZ,vNj,vOmega_start,cGamma_start,mPhi_start,mStep1Var,
ivItemcat,maxIter,tol,fixedpars,1,NRtol,NRmaxit)
out = clean_output4(out,Y,iT,iM,c("Intercept",Z),mY,mZ,id_high,length(Y),P,id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&!is.null(Zh)){
init1 = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat,kmea)
init2 = MLTLCA_poly(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
P_high = ncol(mZh)
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
mDelta_start = matrix(0,iM-1,P_high)
mDelta_start[,1] = init2$vDelta
#
nomissing = complete.cases(mZ)&complete.cases(mZh)
mY = mY[nomissing,]
id_high = id_high[nomissing]
vNj = table(id_high)
mZ = mZ[nomissing,]
mZh = mZh[nomissing,]
mZh = mZh[!duplicated(id_high),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_covlowhigh_poly(mY,mZ,mZh,vNj,mDelta_start,cGamma_start,mPhi_start,
mStep1Var,ivItemcat,maxIter,tol,fixedpars,1,NRtol,NRmaxit)
out = clean_output5(out,Y,iT,iM,c("Intercept",Z),c("Intercept",Zh),mY,mZ,mZh,id_high,length(Y),P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat)
}
} else{
if(is.null(id_high)&is.null(Z)&is.null(Zh)){
out = LCA_fast_init(mY,iT,kmea,maxIter,tol,reord)
out = clean_output1(out,Y,iT,length(Y),extout,dataout,ivItemcat)
} else if(is.null(id_high)&!is.null(Z)&is.null(Zh)){
out = LCA_fast_init_wcov(mY,mZ,iT,kmea,maxIter,tol,fixed,reord,
NRtol,NRmaxit,verbose)
out = clean_output2(out,Y,iT,c("Intercept",Z),length(Y),length(Z)+1,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&is.null(Z)&is.null(Zh)){
init = meas_Init(mY,id_high,vNj,iM,iT,kmea)
out = MLTLCA(mY,vNj,init$vOmega_start,init$mPi_start,init$mPhi_start,
maxIter,tol,reord)
out = clean_output3(out,Y,iT,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&is.null(Zh)){
init1 = meas_Init(mY,id_high,vNj,iM,iT,kmea)
init2 = MLTLCA(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
maxIter,tol,reord)
P = ncol(mZ)
vOmega_start = init2$vOmega
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
#
mY = mY[complete.cases(mZ),]
id_high = id_high[complete.cases(mZ)]
vNj = table(id_high)
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_cov(mY,mZ,vNj,vOmega_start,cGamma_start,mPhi_start,mStep1Var,
maxIter,tol,fixedpars,NRtol,NRmaxit)
out = clean_output4(out,Y,iT,iM,c("Intercept",Z),mY,mZ,id_high,length(Y),P,id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(!is.null(id_high)&!is.null(Z)&!is.null(Zh)){
init1 = meas_Init(mY,id_high,vNj,iM,iT,kmea)
init2 = MLTLCA(mY,vNj,init1$vOmega_start,init1$mPi_start,init1$mPhi_start,
maxIter,tol,reord)
P = ncol(mZ)
P_high = ncol(mZh)
cGamma_start = array(c(rbind(init2$mGamma,matrix(0,(iT-1)*(P-1),iM))),c(iT-1,P,iM))
mPhi_start = init2$mPhi
mStep1Var = init2$Varmat
mDelta_start = matrix(0,iM-1,P_high)
mDelta_start[,1] = init2$vDelta
#
nomissing = complete.cases(mZ)&complete.cases(mZh)
mY = mY[nomissing,]
id_high = id_high[nomissing]
vNj = table(id_high)
mZ = mZ[nomissing,]
mZh = mZh[nomissing,]
mZh = mZh[!duplicated(id_high),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
out = MLTLCA_covlowhigh(mY,mZ,mZh,vNj,mDelta_start,cGamma_start,mPhi_start,
mStep1Var,maxIter,tol,fixedpars,NRtol,NRmaxit)
out = clean_output5(out,Y,iT,iM,c("Intercept",Z),c("Intercept",Zh),mY,mZ,mZh,id_high,length(Y),P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat)
}
}
} else if(approach=="model selection on low"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output1(out$outmuLCA,Y,out$iT_best,length(Y),extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output1(out$outmuLCA,Y,out$iT_best,length(Y),extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on low with high"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,out$iT_best,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iT_best,dimnames=list("iT=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,out$iT_best,iM,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on high"){
if(any(ivItemcat>2)){
out = sel_other_poly(mY,id_high,iT,iM,ivItemcat,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iM),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iM_best,dimnames=list("iM=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,iT,out$iM_best,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
} else{
out = sel_other(mY,id_high,iT,iM,approach,verbose)
out_mat = matrix(round(unlist(out[3:7]),2),length(iM),5,
dimnames=list(paste0("iT=",iT,",iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
out_mat[out_mat=="Inf"|out_mat=="-Inf"|is.na(out_mat)] = "-"
optimal = matrix(out$iM_best,dimnames=list("iM=",""))
list_sel = list(model_selection=out_mat,optimal=optimal)
out = clean_output3(out$outmuLCA,Y,iT,out$iM_best,mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
out$model_selection = list_sel
if(verbose)print(noquote(list_sel))
}
} else if(approach=="model selection on low and high"){
if(sequential){
if(any(ivItemcat>2)){
out = lukosel_fun_poly(mY,id_high,iT,iM,ivItemcat,verbose)
step1mat = matrix(as.character(round(unlist(out[4:8]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step1mat[step1mat=="Inf"|step1mat=="-Inf"|is.na(step1mat)] = "-"
step2mat = matrix(as.character(round(unlist(out[9:13]),2)),length(iM),5,
dimnames=list(paste0("iT*,","iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step2mat[step2mat=="Inf"|step2mat=="-Inf"|is.na(step2mat)] = "-"
step3mat = matrix(as.character(round(unlist(out[14:18]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM*"),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step3mat[step3mat=="Inf"|step3mat=="-Inf"|is.na(step3mat)] = "-"
optimal = matrix(c(out$iT_opt,out$iM_opt),dimnames=list(c("iT=","iM="),""))
list_luk = list(step1 = step1mat,step2 = step2mat,step3 = step3mat, optimal = optimal)
if(out$iM_opt>1){
out = clean_output3(out$outmuLCA_step3,Y,out$iT_opt,out$iM_opt,mY,id_high,length(Y),
id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(out$iM_opt==1){
out = clean_output1(out$outmuLCA_step3,Y,out$iT_opt,length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_luk
if(verbose)print(noquote(list_luk))
} else{
out = lukosel_fun(mY,id_high,iT,iM,verbose)
step1mat = matrix(as.character(round(unlist(out[4:8]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step1mat[step1mat=="Inf"|step1mat=="-Inf"|is.na(step1mat)] = "-"
step2mat = matrix(as.character(round(unlist(out[9:13]),2)),length(iM),5,
dimnames=list(paste0("iT*,","iM=",iM),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step2mat[step2mat=="Inf"|step2mat=="-Inf"|is.na(step2mat)] = "-"
step3mat = matrix(as.character(round(unlist(out[14:18]),2)),length(iT),5,
dimnames=list(paste0("iT=",iT,",iM*"),
c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")))
step3mat[step3mat=="Inf"|step3mat=="-Inf"|is.na(step3mat)] = "-"
optimal = matrix(c(out$iT_opt,out$iM_opt),dimnames=list(c("iT=","iM="),""))
list_luk = list(step1 = step1mat,step2 = step2mat,step3 = step3mat, optimal = optimal)
if(out$iM_opt>1){
out = clean_output3(out$outmuLCA_step3,Y,out$iT_opt,out$iM_opt,mY,id_high,length(Y),
id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(out$iM_opt==1){
out = clean_output1(out$outmuLCA_step3,Y,out$iT_opt,length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_luk
if(verbose)print(noquote(list_luk))
}
} else{
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement models simultaneously...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(any(ivItemcat>2)){
mYfoo = mY
id_highfoo = id_high
ivItemcatfoo = ivItemcat
select_mat = as.matrix(tidyr::expand_grid(iT,iM))
nmod = nrow(select_mat)
iV = parallel::detectCores()-numFreeCores
cluster = parallel::makeCluster(iV)
parallel::clusterExport(cluster,c("mYfoo","id_highfoo","ivItemcatfoo","select_mat","nmod"),
envir = environment())
parallel::clusterExport(cluster,
unclass(lsf.str(envir = asNamespace("multilevLCA"),
all = T)),
envir = as.environment(asNamespace("multilevLCA"))
)
parallel::clusterEvalQ(cluster, {
library(clustMixType)
library(multilevLCA)
library(mclust)
library(tictoc)
library(numDeriv)
library(klaR)
library(tidyverse)
library(MASS)})
simultaneous_out = parallel::parLapply(cluster,1:nmod,
function(x){
iFoo = select_mat[x,]
iT_curr = iFoo[1]
iM_curr = iFoo[2]
out_simult = simultsel_fun_poly(mYfoo,id_highfoo,iT_curr,iM_curr,ivItemcatfoo)
BIClow = out_simult$BIClow
BIChigh = out_simult$BIChigh
AIC = out_simult$AIC
ICLlow = out_simult$ICL_BIClow
ICLhigh = out_simult$ICL_BIChigh
nout = c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")
out = c(BIClow,BIChigh,AIC,ICLlow,ICLhigh)
names(out) = nout
return(list(out=out,fit=out_simult$modFIT))
}
)
parallel::stopCluster(cluster)
mod_sel = apply(t(sapply(simultaneous_out,function(x){x$out})),
2,
function(x){round(as.numeric(x),2)})
rownames(mod_sel) = paste0("iT=",select_mat[,1],",iM=",select_mat[,2])
r_optimal = which.min(mod_sel[,1])
optimal = matrix(c(select_mat[r_optimal,1],select_mat[r_optimal,2]),dimnames=list(c("iT=","iM="),""))
mod_sel[mod_sel==Inf|mod_sel==-Inf|is.na(mod_sel)] = "-"
list_simul = list(model_selection=mod_sel,optimal=optimal)
if(select_mat[r_optimal,2]>1){
out = clean_output3(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],select_mat[r_optimal,2],
mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(select_mat[r_optimal,2]==1){
out = clean_output1(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_simul
if(verbose)print(noquote(list_simul))
} else{
mYfoo = mY
id_highfoo = id_high
select_mat = as.matrix(tidyr::expand_grid(iT,iM))
nmod = nrow(select_mat)
iV = parallel::detectCores()-numFreeCores
cluster = parallel::makeCluster(iV)
parallel::clusterExport(cluster,c("mYfoo","id_highfoo","select_mat","nmod"), envir = environment())
parallel::clusterExport(cluster,
unclass(lsf.str(envir = asNamespace("multilevLCA"),
all = T)),
envir = as.environment(asNamespace("multilevLCA"))
)
parallel::clusterEvalQ(cluster, {
library(clustMixType)
library(multilevLCA)
library(mclust)
library(tictoc)
library(numDeriv)
library(klaR)
library(tidyverse)
library(MASS)})
simultaneous_out = parallel::parLapply(cluster,1:nmod,
function(x){
iFoo = select_mat[x,]
iT_curr = iFoo[1]
iM_curr = iFoo[2]
out_simult = simultsel_fun(mYfoo,id_highfoo,iT_curr,iM_curr)
BIClow = out_simult$BIClow
BIChigh = out_simult$BIChigh
AIC = out_simult$AIC
ICLlow = out_simult$ICL_BIClow
ICLhigh = out_simult$ICL_BIChigh
nout = c("BIClow","BIChigh","AIC","ICL_BIClow","ICL_BIChigh")
out = c(BIClow,BIChigh,AIC,ICLlow,ICLhigh)
names(out) = nout
return(list(out=out,fit=out_simult$modFIT))
}
)
parallel::stopCluster(cluster)
mod_sel = apply(t(sapply(simultaneous_out,function(x){x$out})),
2,
function(x){round(as.numeric(x),2)})
rownames(mod_sel) = paste0("iT=",select_mat[,1],",iM=",select_mat[,2])
r_optimal = which.min(mod_sel[,1])
optimal = matrix(c(select_mat[r_optimal,1],select_mat[r_optimal,2]),dimnames=list(c("iT=","iM="),""))
mod_sel[mod_sel==Inf|mod_sel==-Inf|is.na(mod_sel)] = "-"
list_simul = list(model_selection=mod_sel,optimal=optimal)
if(select_mat[r_optimal,2]>1){
out = clean_output3(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],select_mat[r_optimal,2],
mY,id_high,length(Y),id_high_levs,id_high_name,extout,dataout,ivItemcat)
} else if(select_mat[r_optimal,2]==1){
out = clean_output1(simultaneous_out[[r_optimal]]$fit,
Y,select_mat[r_optimal,1],length(Y),extout,dataout,ivItemcat)
}
out$model_selection = list_simul
if(verbose)print(noquote(list_simul))
}
}
}
if(nrow(data)<nrow_data){
if(fixed==1|fixed==2){
warning(paste("Missing values in columns for indicators,",
"sample size for estimation of measurement model:",
nrow(data)),call.=FALSE)
} else if(fixed==0&is.null(Z)){
warning(paste("Missing values in columns for indicators,",
"sample size for estimation of measurement model:",
nrow(data)),call.=FALSE)
} else if(fixed==0&!is.null(Z)){
if(nrow(na.omit(mZ))<nrow(data)){
warning(paste("Missing values in columns for indicators",
"and columns for covariates,",
"sample size for estimation:",
nrow(na.omit(mZ))),call.=FALSE)
} else{
warning(paste("Missing values in columns for indicators,",
"sample size for estimation:",
nrow(data)),call.=FALSE)
}
}
}
if(approach=="direct"&!is.null(Z)){
if(nrow(na.omit(mZ))<nrow(data)){
if(fixed==1|fixed==2){
warning(paste("Missing values in columns for covariates,",
"sample size for estimation of structural model:",
nrow(na.omit(mZ))),call.=FALSE)
} else if(fixed==0){
if(nrow(data)==nrow_data){
warning(paste("Missing values in columns for covariates,",
"sample size for estimation of structural model:",
nrow(na.omit(mZ))),call.=FALSE)
}
}
}
}
out$call = match.call()
class(out) = "multiLCA"
return(out)
}
#
LCA_fast_init = function(mY, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
out$mY_unique = mY_unique
out$freq = freq
return(out)
}
#
LCA_fast_init_wcov = function(mY, mZ, iT, kmea = T, maxIter = 1e3, tol = 1e-8, fixed = 0, reord = 1,
NRtol = 1e-6, NRmaxit = 100, verbose){
# mZ must include a column of ones!
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
P = ncol(mZ)
mBeta_init = matrix(0,P,iT-1)
mBeta_init[1,] = out$alphas
Step1Var = out$Varmat
mY = mY[complete.cases(mZ),]
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
outcov = LCAcov(mY,mZ,iT,out$mPhi,mBeta_init,Step1Var,fixed,maxIter,
tol,NRtol,NRmaxit)
return(list(out=out,outcov=outcov,mY=mY,mZ=mZ))
}
#
LCA_fast_init_whigh = function(mY, id_high, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
group_by_all = NULL
mY_df = data.frame(mY,id_high)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
mY_aggr = mY_aggr[order(mY_aggr[,iHf-1]),]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-c(iHf-1,iHf)]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast(mY_unique,freq,iT,mU,maxIter,tol,reord)
mU = out$mU[rep(1:nrow(mY_unique),freq),]
return(list(out=out,mU=mU))
}
#
meas_Init = function(mY, id_high, vNj, iM, iT, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# Fixed number of low- (iT) and high- (iM) level classes
iJ = length(vNj)
iN = dim(mY)[1]
iK = dim(mY)[2]
#
# Working out starting values at the low level first
#
out_LCA = LCA_fast_init_whigh(mY,id_high,iT,kmea,maxIter,tol,reord)
#
# Now turning to the higher level
#
lowlev_relclassprop = matrix(0,iJ,iT)
foo = 0
for(j in 1:iJ){
lowlev_relclassprop[j,] = colMeans(out_LCA$mU[foo+(1:vNj[j]),,drop=FALSE])
foo = foo +vNj[j]
}
high_out = kmeans(lowlev_relclassprop,centers = iM,iter.max=100,nstart=100)
vOmega_start = high_out$size/iJ
Wmodal_mat = vecTomatClass(high_out$cluster)
# Reordering in decreasing order
highreord = order(vOmega_start,decreasing = T)
Wmodal_mat_reord = Wmodal_mat[,highreord]
Wmodal = apply(Wmodal_mat_reord,1,which.max)
vOmega_start = vOmega_start[highreord]
index_indiv_high = rep(Wmodal,times=vNj)
#
mPhi_start = out_LCA$out$mPhi
mPi_fast = table(index_indiv_high,rep(out_LCA$out$vModalAssnm+1,times=out_LCA$out$freq))
mPi_start = t(mPi_fast/rowSums(mPi_fast))
if(nrow(mPi_start)!=iT){
warning(paste0("Initialization suggests less than iT=",iT," lower-level classes."),call.=FALSE)
mPi_start_replace = matrix(NA,nrow=iT,ncol=iM)
mPi_start_replace[as.numeric(rownames(mPi_start)),] = mPi_start
mPi_start_replace[is.na(mPi_start_replace)] = 1e-4
mPi_start_replace = apply(mPi_start_replace,2,function(x){x/sum(x)})
mPi_start = mPi_start_replace
}
#
return(list(vOmega_start=vOmega_start, mPi_start = mPi_start, mPhi_start = mPhi_start))
}
#
LCA_fast_init_poly = function(mY, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
out$mY_unique = mY_unique
out$freq = freq
return(out)
}
#
LCA_fast_init_wcov_poly = function(mY, mZ, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, fixed = 0, reord = 1,
NRtol = 1e-6, NRmaxit = 100, verbose){
# mZ must include a column of ones!
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-iHf]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
P = ncol(mZ)
mBeta_init = matrix(0,P,iT-1)
mBeta_init[1,] = out$alphas
Step1Var = out$Varmat
mY = mY[complete.cases(mZ),]
mZ = mZ[complete.cases(mZ),]
#
if(verbose){
if(fixed==1|fixed==2){
pb=txtProgressBar(char="Fitting structural model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
}
#
outcov = LCAcov_poly(mY,mZ,iT,out$mPhi,mBeta_init,Step1Var,ivItemcat,fixed,maxIter,
tol,NRtol,NRmaxit)
return(list(out=out,outcov=outcov,mY=mY,mZ=mZ))
}
#
LCA_fast_init_whigh_poly = function(mY, id_high, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# ivItemcat is the vector of number of categories for each item
group_by_all = NULL
mY_df = data.frame(mY,id_high)
mY_aggr = as.matrix(mY_df%>%group_by_all%>%count)
iHf = dim(mY_aggr)[2]
mY_aggr = mY_aggr[order(mY_aggr[,iHf-1]),]
freq = mY_aggr[,iHf]
mY_unique = mY_aggr[,-c(iHf-1,iHf)]
if(kmea==FALSE){
clusfoo = klaR::kmodes(mY_unique,modes=iT,iter.max=100)$cluster
} else{
prscores = prcomp(mY_unique)
num_pc = max(round(ncol(mY_unique)/2),(sum(cumsum(prscores$sdev^2/sum(prscores$sdev^2))<0.85)+1))
prscores = prscores$x[,1:num_pc]
spectclust = kmeans(prscores,centers=iT,iter.max=100,nstart=100)
clusfoo = spectclust$cluster
}
mU = vecTomatClass(clusfoo)
out = LCA_fast_poly(mY_unique,freq,iT,mU,ivItemcat,maxIter,tol,reord)
mU = out$mU[rep(1:nrow(mY_unique),freq),]
return(list(out=out,mU=mU))
}
#
meas_Init_poly = function(mY, id_high, vNj, iM, iT, ivItemcat, kmea = T, maxIter = 1e3, tol = 1e-8, reord = 1){
# Fixed number of low- (iT) and high- (iM) level classes
# ivItemcat is the vector of number of categories for each item
iJ = length(vNj)
iN = dim(mY)[1]
iK = dim(mY)[2]
#
# Working out starting values at the low level first
#
out_LCA = LCA_fast_init_whigh_poly(mY,id_high,iT,ivItemcat,kmea,maxIter,tol,reord)
#
# now turning to higher level
#
lowlev_relclassprop = matrix(0,iJ,iT)
foo = 0
for(j in 1:iJ){
lowlev_relclassprop[j,] = colMeans(out_LCA$mU[foo+(1:vNj[j]),,drop=FALSE])
foo = foo +vNj[j]
}
high_out = kmeans(lowlev_relclassprop,centers = iM,iter.max=100,nstart=100)
vOmega_start = high_out$size/iJ
Wmodal_mat = vecTomatClass(high_out$cluster)
# Reordering in decreasing order
highreord = order(vOmega_start,decreasing = T)
Wmodal_mat_reord = Wmodal_mat[,highreord]
Wmodal = apply(Wmodal_mat_reord,1,which.max)
vOmega_start = vOmega_start[highreord]
index_indiv_high = rep(Wmodal,times=vNj)
#
mPhi_start = out_LCA$out$mPhi
mPi_fast = table(index_indiv_high,rep(out_LCA$out$vModalAssnm+1,times=out_LCA$out$freq))
mPi_start = t(mPi_fast/rowSums(mPi_fast))
if(nrow(mPi_start)!=iT){
warning(paste0("Initialization suggests less than iT=",iT," lower-level classes."),call.=FALSE)
mPi_start_replace = matrix(NA,nrow=iT,ncol=iM)
mPi_start_replace[as.numeric(rownames(mPi_start)),] = mPi_start
mPi_start_replace[is.na(mPi_start_replace)] = 1e-4
mPi_start_replace = apply(mPi_start_replace,2,function(x){x/sum(x)})
mPi_start = mPi_start_replace
}
#
return(list(vOmega_start=vOmega_start, mPi_start = mPi_start, mPhi_start = mPhi_start))
}
#
#
#
simultsel_fun = function(mY,id_high,iT,iM){
LCAout = NULL
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
if(iM ==1 & iT ==1){
ll = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow = -2*ll + K*log(sum(vNj))
BIChigh = -2*ll + K*log(iN)
AIC = -2*ll + 2*K
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM > 1 & iT == 1){
ll = -Inf
BIClow = Inf
BIChigh = Inf
AIC = Inf
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM==1 & iT > 1){
LCAout = LCA_fast_init(mY,iT)
ll = tail(LCAout$LLKSeries,1)
npar = iT*iH + iT - 1
BIClow = LCAout$BIC
BIChigh = -2*ll + npar*log(iN)
AIC = -2*ll + 2*npar
ICL_BIClow = BIClow + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh = Inf
outmuLCA = LCAout
}else{
# Actual multilevel LC fit
start = meas_Init(mY,id_high,vNj,iM,iT)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA = MLTLCA(mY, vNj, vOmegast, mPi, start$mPhi_start)
ll = tail(outmuLCA$LLKSeries,1)
BIClow = outmuLCA$BIClow
BIChigh = outmuLCA$BIChigh
AIC = outmuLCA$AIC
ICL_BIClow = outmuLCA$ICL_BIClow
ICL_BIChigh = outmuLCA$ICL_BIChigh
}
return(list(modFIT = outmuLCA,ll=ll, BIClow = BIClow, BIChigh = BIChigh, AIC = AIC,
ICL_BIClow = ICL_BIClow, ICL_BIChigh = ICL_BIChigh))
}
#
lukosel_fun = function(mY,id_high,iT_range,iM_range,verbose){
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
# Step 1 - lower level
LCAout = list()
BIClow_step1 = rep(NA,num_iT)
BIChigh_step1 = rep(NA,num_iT)
AIC_step1 = rep(NA,num_iT)
ICL_BIClow_step1 = rep(NA,num_iT)
ICL_BIChigh_step1 = rep(NA,num_iT)
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting single-level measurement model (step 1)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow_step1[1] = -2*llfoo + K*log(sum(vNj))
BIChigh_step1[1] = -2*llfoo + K*log(iN)
AIC_step1[1] = -2*llfoo + 2*K
ICL_BIClow_step1[1] = Inf
ICL_BIChigh_step1[1] = Inf
} else{
LCAout = LCA_fast_init(mY,i)
ll = tail(LCAout$LLKSeries,1)
npar = i*iH + i - 1
BIClow_step1[1+i-iT_min] = LCAout$BIC
BIChigh_step1[1+i-iT_min] = -2*ll + npar*log(iN)
AIC_step1[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow_step1[1+i-iT_min] = BIClow_step1[i] + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh_step1[1+i-iT_min] = Inf
}
}
iT_currbest = which.min(BIClow_step1)+iT_min-1
if(iT_currbest==1) stop("iT=1 optimal.",call.=FALSE)
# Step 2 - higher level
outmuLCA2 = list()
BIClow_step2 = rep(NA,num_iM)
BIChigh_step2 = rep(NA,num_iM)
AIC_step2 = rep(NA,num_iM)
ICL_BIClow_step2 = rep(NA,num_iM)
ICL_BIChigh_step2 = rep(NA,num_iM)
if(iT_currbest > 1){
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 2)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
BIClow_step2[1] = BIClow_step1[1+iT_currbest-iT_min]
BIChigh_step2[1] = BIChigh_step1[1+iT_currbest-iT_min]
AIC_step2[1] = AIC_step1[1+iT_currbest-iT_min]
ICL_BIClow_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
ICL_BIChigh_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
} else{
start = meas_Init(mY,id_high,vNj,i,iT_currbest)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA2[[1+i-iM_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIClow
BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIChigh
AIC_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$AIC
ICL_BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_currbest = which.min(BIChigh_step2)+iM_min-1
} else{
iM_currbest = 1
}
# Step 3 - revisiting lower level
outmuLCA3 = list()
BIClow_step3 = rep(NA,num_iT)
BIChigh_step3 = rep(NA,num_iT)
AIC_step3 = rep(NA,num_iT)
ICL_BIClow_step3 = rep(NA,num_iT)
ICL_BIChigh_step3 = rep(NA,num_iT)
if(iM_currbest > 1){
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 3)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init(mY,id_high,vNj,iM_currbest,i)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA3[[1+i-iT_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIClow
BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIChigh
AIC_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$AIC
ICL_BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_currbest = which(BIClow_step3==min(BIClow_step3,na.rm=T))+iT_min-1
outmuLCA_step3 = outmuLCA3[[1+iT_currbest-iT_min]]
} else{
outmuLCA_step3 = LCA_fast_init(mY,iT_currbest)
}
return(list(outmuLCA_step3=outmuLCA_step3,iT_opt=iT_currbest, iM_opt=iM_currbest,
BIClow_step1 = BIClow_step1, BIChigh_step1 = BIChigh_step1, AIC_step1 = AIC_step1,
ICL_BIClow_step1 = ICL_BIClow_step1, ICL_BIChigh_step1 = ICL_BIChigh_step1,
BIClow_step2 = BIClow_step2, BIChigh_step2 = BIChigh_step2, AIC_step2 = AIC_step2,
ICL_BIClow_step2 = ICL_BIClow_step2, ICL_BIChigh_step2 = ICL_BIChigh_step2,
BIClow_step3 = BIClow_step3, BIChigh_step3 = BIChigh_step3, AIC_step3 = AIC_step3,
ICL_BIClow_step3 = ICL_BIClow_step3, ICL_BIChigh_step3 = ICL_BIChigh_step3))
}
#
sel_other = function(mY,id_high,iT_range,iM_range,approach,verbose){
iH = ncol(mY)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
if(!is.null(iM_range)){
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
if(num_iM==1){
if(iM_range==1){
iM_range=NULL
}
}
}
#
if(approach=="model selection on low"|approach=="model selection on low with high"){
BIClow = rep(NA,num_iT)
BIChigh = rep(NA,num_iT)
AIC = rep(NA,num_iT)
ICL_BIClow = rep(NA,num_iT)
ICL_BIChigh = rep(NA,num_iT)
if(is.null(iM_range)){
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = NA
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow[1] = -2*llfoo + K*log(sum(vNj))
if(is.null(iM_range)) BIClow[1] = -2*llfoo + K*log(nrow(mY))
if(!is.null(iM_range)) BIChigh[1] = -2*llfoo + K*log(iN)
AIC[1] = -2*llfoo + 2*K
ICL_BIClow[1] = Inf
ICL_BIChigh[1] = Inf
} else{
outmuLCA[[1+i-iT_min]] = LCA_fast_init(mY,i)
ll = tail(outmuLCA[[1+i-iT_min]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIC
BIChigh[1+i-iT_min] = -2*ll + npar*log(iN)
AIC[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow[1+i-iT_min] = BIClow[i] + 2*sum(outmuLCA[[1+i-iT_min]]$freq*apply(-outmuLCA[[1+i-iT_min]]$mU*log(outmuLCA[[1+i-iT_min]]$mU),1,sum))
ICL_BIChigh[1+i-iT_min] = Inf
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
} else{
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init(mY,id_high,vNj,iM_range,i)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iT_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIClow
BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIChigh
AIC[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$AIC
ICL_BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
}
if(iT_best==1) stop("iT=1 optimal.",call.=FALSE)
return(list(outmuLCA=outmuLCA,iT_best=iT_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
} else if(approach=="model selection on high"){
BIClow = rep(NA,num_iM)
BIChigh = rep(NA,num_iM)
AIC = rep(NA,num_iM)
ICL_BIClow = rep(NA,num_iM)
ICL_BIChigh = rep(NA,num_iM)
if(iT_range > 1){
outmuLCA = list()
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = LCA_fast_init(mY,iT_range)
ll = tail(outmuLCA[[1]]$LLKSeries,1)
npar = iT_range*iH + iT_range - 1
BIClow[1] = outmuLCA[[1]]$BIC
BIChigh[1] = -2*ll + npar*log(iN)
AIC[1] = -2*ll + 2*npar
ICL_BIClow[1] = BIClow[1] + 2*sum(outmuLCA[[1]]$freq*apply(-outmuLCA[[1]]$mU*log(outmuLCA[[1]]$mU),1,sum))
ICL_BIChigh[1] = Inf
} else{
start = meas_Init(mY,id_high,vNj,i,iT_range)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iM_min]] = MLTLCA(mY,vNj,vOmegast,mPi,start$mPhi_start)
BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIClow
BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIChigh
AIC[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$AIC
ICL_BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_best = which.min(BIChigh)+iM_min-1
outmuLCA = outmuLCA[[1+iM_best-iM_min]]
} else{
stop("iM=1 optimal.",call.=FALSE)
}
return(list(outmuLCA=outmuLCA,iM_best=iM_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
}
}
#
simultsel_fun_poly = function(mY,id_high,iT,iM,ivItemcat){
LCAout = NULL
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
if(iM ==1 & iT ==1){
ll = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow = -2*ll + K*log(sum(vNj))
BIChigh = -2*ll + K*log(iN)
AIC = -2*ll + 2*K
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM > 1 & iT == 1){
ll = -Inf
BIClow = Inf
BIChigh = Inf
AIC = Inf
ICL_BIClow = Inf
ICL_BIChigh = Inf
outmuLCA = list()
}else if(iM==1 & iT > 1){
LCAout = LCA_fast_init_poly(mY,iT,ivItemcat)
ll = tail(LCAout$LLKSeries,1)
npar = iT*iH + iT - 1
BIClow = LCAout$BIC
BIChigh = -2*ll + npar*log(iN)
AIC = -2*ll + 2*npar
ICL_BIClow = BIClow + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh = Inf
outmuLCA = LCAout
}else{
# Actual multilevel LC fit
start = meas_Init_poly(mY,id_high,vNj,iM,iT,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA = MLTLCA_poly(mY,vNj,vOmegast,mPi,
start$mPhi_start,ivItemcat)
ll = tail(outmuLCA$LLKSeries,1)
BIClow = outmuLCA$BIClow
BIChigh = outmuLCA$BIChigh
AIC = outmuLCA$AIC
ICL_BIClow = outmuLCA$ICL_BIClow
ICL_BIChigh = outmuLCA$ICL_BIChigh
}
return(list(modFIT = outmuLCA,ll=ll, BIClow = BIClow, BIChigh = BIChigh, AIC = AIC,
ICL_BIClow = ICL_BIClow, ICL_BIChigh = ICL_BIChigh))
}
#
lukosel_fun_poly = function(mY,id_high,iT_range,iM_range,ivItemcat,verbose){
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
# Step 1 - lower level
LCAout = list()
BIClow_step1 = rep(NA,num_iT)
BIChigh_step1 = rep(NA,num_iT)
AIC_step1 = rep(NA,num_iT)
ICL_BIClow_step1 = rep(NA,num_iT)
ICL_BIChigh_step1 = rep(NA,num_iT)
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting single-level measurement model (step 1)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow_step1[1] = -2*llfoo + K*log(sum(vNj))
BIChigh_step1[1] = -2*llfoo + K*log(iN)
AIC_step1[1] = -2*llfoo + 2*K
ICL_BIClow_step1[1] = Inf
ICL_BIChigh_step1[1] = Inf
} else{
LCAout = LCA_fast_init_poly(mY,i,ivItemcat)
ll = tail(LCAout$LLKSeries,1)
npar = i*iH + i - 1
BIClow_step1[1+i-iT_min] = LCAout$BIC
BIChigh_step1[1+i-iT_min] = -2*ll + npar*log(iN)
AIC_step1[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow_step1[1+i-iT_min] = BIClow_step1[i] + 2*sum(LCAout$freq*apply(-LCAout$mU*log(LCAout$mU),1,sum))
ICL_BIChigh_step1[1+i-iT_min] = Inf
}
}
iT_currbest = which.min(BIClow_step1)+iT_min-1
if(iT_currbest==1) stop("iT=1 optimal.",call.=FALSE)
# Step 2 - higher level
outmuLCA2 = list()
BIClow_step2 = rep(NA,num_iM)
BIChigh_step2 = rep(NA,num_iM)
AIC_step2 = rep(NA,num_iM)
ICL_BIClow_step2 = rep(NA,num_iM)
ICL_BIChigh_step2 = rep(NA,num_iM)
if(iT_currbest > 1){
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 2)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
BIClow_step2[1] = BIClow_step1[1+iT_currbest-iT_min]
BIChigh_step2[1] = BIChigh_step1[1+iT_currbest-iT_min]
AIC_step2[1] = AIC_step1[1+iT_currbest-iT_min]
ICL_BIClow_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
ICL_BIChigh_step2[1] = ICL_BIClow_step1[1+iT_currbest-iT_min]
} else{
start = meas_Init_poly(mY,id_high,vNj,i,iT_currbest,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA2[[1+i-iM_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIClow
BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$BIChigh
AIC_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$AIC
ICL_BIClow_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh_step2[1+i-iM_min] = outmuLCA2[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_currbest = which.min(BIChigh_step2)+iM_min-1
} else{
iM_currbest = 1
}
# Step 3 - revisiting lower level
outmuLCA3 = list()
BIClow_step3 = rep(NA,num_iT)
BIChigh_step3 = rep(NA,num_iT)
AIC_step3 = rep(NA,num_iT)
ICL_BIClow_step3 = rep(NA,num_iT)
ICL_BIChigh_step3 = rep(NA,num_iT)
if(iM_currbest > 1){
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting multilevel measurement model (step 3)...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>2){
start = meas_Init_poly(mY,id_high,vNj,iM_currbest,i,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA3[[1+i-iT_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIClow
BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$BIChigh
AIC_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$AIC
ICL_BIClow_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh_step3[1+i-iT_min] = outmuLCA3[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_currbest = which(BIClow_step3==min(BIClow_step3,na.rm=T))+iT_min-1
outmuLCA_step3 = outmuLCA3[[1+iT_currbest-iT_min]]
} else{
outmuLCA_step3 = LCAout
}
return(list(outmuLCA_step3=outmuLCA_step3,iT_opt=iT_currbest, iM_opt=iM_currbest,
BIClow_step1 = BIClow_step1, BIChigh_step1 = BIChigh_step1, AIC_step1 = AIC_step1,
ICL_BIClow_step1 = ICL_BIClow_step1, ICL_BIChigh_step1 = ICL_BIChigh_step1,
BIClow_step2 = BIClow_step2, BIChigh_step2 = BIChigh_step2, AIC_step2 = AIC_step2,
ICL_BIClow_step2 = ICL_BIClow_step2, ICL_BIChigh_step2 = ICL_BIChigh_step2,
BIClow_step3 = BIClow_step3, BIChigh_step3 = BIChigh_step3, AIC_step3 = AIC_step3,
ICL_BIClow_step3 = ICL_BIClow_step3, ICL_BIChigh_step3 = ICL_BIChigh_step3))
}
#
sel_other_poly = function(mY,id_high,iT_range,iM_range,ivItemcat,approach,verbose){
iH = sum(ivItemcat-1)
K = iH
vNj = table(id_high-1)
iN = length(vNj)
iT_max = max(iT_range)
iT_min = min(iT_range)
num_iT = iT_max-iT_min+1
if(!is.null(iM_range)){
iM_max = max(iM_range)
iM_min = min(iM_range)
num_iM = iM_max-iM_min+1
if(num_iM==1){
if(iM_range==1){
iM_range=NULL
}
}
}
#
if(approach=="model selection on low"|approach=="model selection on low with high"){
BIClow = rep(NA,num_iT)
BIChigh = rep(NA,num_iT)
AIC = rep(NA,num_iT)
ICL_BIClow = rep(NA,num_iT)
ICL_BIChigh = rep(NA,num_iT)
if(is.null(iM_range)){
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = NA
llfoo = sum(apply(mY,2,function(x){dbinom(x,1,mean(x),log=T)}))
BIClow[1] = -2*llfoo + K*log(sum(vNj))
if(is.null(iM_range)) BIClow[1] = -2*llfoo + K*log(nrow(mY))
if(!is.null(iM_range)) BIChigh[1] = -2*llfoo + K*log(iN)
AIC[1] = -2*llfoo + 2*K
ICL_BIClow[1] = Inf
ICL_BIChigh[1] = Inf
} else{
outmuLCA[[1+i-iT_min]] = LCA_fast_init_poly(mY,i,ivItemcat)
ll = tail(outmuLCA[[1+i-iT_min]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIC
BIChigh[1+i-iT_min] = -2*ll + npar*log(iN)
AIC[1+i-iT_min] = -2*ll + 2*npar
ICL_BIClow[1+i-iT_min] = BIClow[i] + 2*sum(outmuLCA[[1+i-iT_min]]$freq*apply(-outmuLCA[[1+i-iT_min]]$mU*log(outmuLCA[[1+i-iT_min]]$mU),1,sum))
ICL_BIChigh[1+i-iT_min] = Inf
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
} else{
outmuLCA = list()
for(i in iT_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i>1){
start = meas_Init_poly(mY,id_high,vNj,iM_range,i,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iT_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIClow
BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$BIChigh
AIC[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$AIC
ICL_BIClow[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIClow
ICL_BIChigh[1+i-iT_min] = outmuLCA[[1+i-iT_min]]$ICL_BIChigh
}
}
iT_best = which(BIClow==min(BIClow,na.rm=T))+iT_min-1
outmuLCA = outmuLCA[[1+iT_best-iT_min]]
}
if(iT_best==1) stop("iT=1 optimal.",call.=FALSE)
return(list(outmuLCA=outmuLCA,iT_best=iT_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
} else if(approach=="model selection on high"){
BIClow = rep(NA,num_iM)
BIChigh = rep(NA,num_iM)
AIC = rep(NA,num_iM)
ICL_BIClow = rep(NA,num_iM)
ICL_BIChigh = rep(NA,num_iM)
if(iT_range > 1){
outmuLCA = list()
for(i in iM_range){
#
if(verbose){
pb=txtProgressBar(char="Fitting measurement model...",width=1)
setTxtProgressBar(pb,1)
close(pb)
}
#
if(i==1){
outmuLCA[[1]] = LCA_fast_init_poly(mY,iT_range,ivItemcat)
ll = tail(outmuLCA[[1]]$LLKSeries,1)
npar = i*iH + i - 1
BIClow[1] = outmuLCA[[1]]$BIC
BIChigh[1] = -2*ll + npar*log(iN)
AIC[1] = -2*ll + 2*npar
ICL_BIClow[1] = BIClow[i] + 2*sum(outmuLCA[[1]]$freq*apply(-outmuLCA[[1]]$mU*log(outmuLCA[[1]]$mU),1,sum))
ICL_BIChigh[1] = Inf
} else{
start = meas_Init_poly(mY,id_high,vNj,i,iT_range,ivItemcat)
vOmegast = start$vOmega_start
mPi = start$mPi_start
outmuLCA[[1+i-iM_min]] = MLTLCA_poly(mY,vNj,vOmegast,mPi,start$mPhi_start,ivItemcat)
BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIClow
BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$BIChigh
AIC[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$AIC
ICL_BIClow[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIClow
ICL_BIChigh[1+i-iM_min] = outmuLCA[[1+i-iM_min]]$ICL_BIChigh
}
}
iM_best = which.min(BIChigh)+iM_min-1
outmuLCA = outmuLCA[[1+iM_best-iM_min]]
} else{
stop("iM=1 optimal.",call.=FALSE)
}
return(list(outmuLCA=outmuLCA,iM_best=iM_best,
BIClow=BIClow,BIChigh=BIChigh,AIC=AIC,
ICL_BIClow=ICL_BIClow,ICL_BIChigh=ICL_BIChigh))
}
}
#
#
#
clean_output1 = function(output,Y,iT,iH,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add vector of class proportions
out$vPi = output$pg
rownames(out$vPi) = paste0("P(C",1:iT,")")
colnames(out$vPi) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add vector of class proportions
out$vPi = output$pg
rownames(out$vPi) = paste0("P(C",1:iT,")")
colnames(out$vPi) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add matrix of posterior class membership probabilities
out$mU = output$mU[rep(1:dim(output$mU)[1],output$freq),]
colnames(out$mU) = colnames(out$mPhi)
if(dataout){
out$mU = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$mU)
rownames(out$mU) = NULL
}
# Add matrix of modal class assignment
out$mU_modal = output$mModalAssnm[rep(1:dim(output$mModalAssnm)[1],output$freq),]
colnames(out$mU_modal) = colnames(out$mPhi)
if(dataout){
out$mU_modal = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$mU_modal)
rownames(out$mU_modal) = NULL
}
# Add vector of modal class assignment
out$vU_modal = output$vModalAssnm[rep(1:dim(output$vModalAssnm)[1],output$freq),,drop=FALSE]+1
colnames(out$vU_modal) = "C"
if(dataout){
out$vU_modal = cbind(output$mY_unique[rep(1:dim(output$mY_unique)[1],output$freq),],out$vU_modal)
rownames(out$vU_modal) = NULL
}
# Add matrix of classification errors
out$mClassErr = output$mClassErr
rownames(out$mClassErr) = paste0("C",1:iT,"_true")
colnames(out$mClassErr) = paste0("C",1:iT,"_pred")
# Add matrix of proportion of classification errors
out$mClassErrProb = output$mClassErrProb
rownames(out$mClassErrProb) = rownames(out$mClassErr)
colnames(out$mClassErrProb) = colnames(out$mClassErr)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$vGamma = output$alphas
rownames(out$vGamma) = paste0("gamma(C",2:iT,")")
colnames(out$vGamma) = ""
# Add betas
out$mBeta = output$gamma
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vGamma),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of standard errors
out$SEs = output$SEs
rownames(out$SEs) = rownames(out$parvec)
colnames(out$SEs) = ""
# Add variance-covariance matrix
out$Varmat = output$Varmat
rownames(out$Varmat) = colnames(out$Varmat) = rownames(out$parvec)
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add specification
out$spec = as.matrix("Single-level LC model")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output2 = function(output,Y,iT,Z,iH,P,extout,dataout,ivItemcat){
mY = output$mY
mZ = output$mZ
output = output$outcov
if(!extout){
# Create empty list for output
out = list()
# Add vector of sample means of class proportions
out$vPi_avg = as.matrix(apply(output$mPg,2,mean))
rownames(out$vPi_avg) = paste0("P(C",1:iT,")")
colnames(out$vPi_avg) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = output$beta
rownames(out$cGamma) = paste0("gamma(",Z,"|C)")
colnames(out$cGamma) = paste0("C",2:iT)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = matrix(output$SEs_cor[1:((iT-1)*P),],P,iT-1)
rownames(out$SEs_cor_gamma) = rownames(out$cGamma)
colnames(out$SEs_cor_gamma) = colnames(out$cGamma)
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model with covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of class proportions
out$mPi = output$mPg
colnames(out$mPi) = paste0("P(C",1:iT,")")
# Add vector of sample means of class proportions
out$vPi_avg = as.matrix(apply(output$mPg,2,mean))
rownames(out$vPi_avg) = colnames(out$mPi)
colnames(out$vPi_avg) = ""
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add matrix of posterior class membership probabilities
out$mU = output$mU
colnames(out$mU) = colnames(out$mPhi)
if(dataout){
data = cbind(mY,mZ)
rownames(data) = NULL
colnames(data) = c(Y,Z)
data = data[,-which(colnames(data)=="Intercept"),drop=FALSE]
out$mU = cbind(data,out$mU)
rownames(out$mU) = NULL
}
# Add matrix of classification errors
out$mClassErr = output$mClassErr
rownames(out$mClassErr) = paste0("C",1:iT,"_true")
colnames(out$mClassErr) = paste0("C",1:iT,"_pred")
# Add matrix of proportion of classification errors
out$mClassErrProb = output$mClassErrProb
rownames(out$mClassErrProb) = paste0("C",1:iT,"_true")
colnames(out$mClassErrProb) = paste0("C",1:iT,"_pred")
# Add average proportion of classification errors
out$AvgClassErrProb = output$dClassErr_tot
# Add entropy R-sqr
out$R2entr = output$R2entr
# Add Bayesian information criterion
out$BIC = output$BIC
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = output$beta
rownames(out$cGamma) = paste0("gamma(",Z,"|C)")
colnames(out$cGamma) = paste0("C",2:iT)
# Add betas
out$mBeta = output$gamma
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(paste0(rep(substr(rownames(out$cGamma),1,nchar(rownames(out$cGamma))-2),iT-1),rep(colnames(out$cGamma),rep(P,iT-1)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:((iT-1)*P),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:((iT-1)*P),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for gamma
out$SEs_cor_gamma = matrix(output$SEs_cor[1:((iT-1)*P),],P,iT-1)
rownames(out$SEs_cor_gamma) = rownames(out$cGamma)
colnames(out$SEs_cor_gamma) = colnames(out$cGamma)
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mQ) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat_unc[1:((iT-1)*P),1:((iT-1)*P)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$Varmat_cor
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$mV2)
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = paste0(rep(substr(rownames(out$cGamma),1,nchar(rownames(out$cGamma))-2),iT-1),rep(colnames(out$cGamma),rep(P,iT-1)),")")
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Single-level LC model with covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output3 = function(output,Y,iT,iM,mY,id_high,iH,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = output$mPi
rownames(out$mPi) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = output$mPi
rownames(out$mPi) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi)))
if(dataout){
data = cbind(mY,id_high)
rownames(data) = NULL
colnames(data) = c(Y,id_high_name)
for(i in 1:length(id_high_levs)){
data[data[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX,dim(output$cLogPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX,dim(output$cPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi)
if(dataout){
out$mSumPX = cbind(data,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi)
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi)
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi)
if(dataout){
out$mPW_N = cbind(data,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$vAlpha = output$vDelta
rownames(out$vAlpha) = paste0("alpha(G",2:iM,")")
colnames(out$vAlpha) = ""
# Add gammas
out$mGamma = output$mGamma
rownames(out$mGamma) = paste0("gamma(C",2:iT,"|G)")
colnames(out$mGamma) = colnames(out$mPi)
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vAlpha),paste0(rep(substr(rownames(out$mGamma),1,nchar(rownames(out$mGamma))-2),iM),rep(colnames(out$mGamma),rep(iT-1,iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of standard errors
out$SEs = output$SEs
rownames(out$SEs) = rownames(out$parvec)
colnames(out$SEs) = ""
# Add variance-covariance matrix
out$Varmat = output$Varmat
rownames(out$Varmat) = colnames(out$Varmat) = rownames(out$parvec)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = colnames(out$Infomat) = rownames(out$parvec)
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add specification
out$spec = as.matrix("Multilevel LC model")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output4 = function(output,Y,iT,iM,Z,mY,mZ,id_high,iH,P,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add vector of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[iM:(iM-1+(P*(iT-1)*iM))],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model with lower-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$vOmega = output$vOmega
rownames(out$vOmega) = paste0("P(G",1:iM,")")
colnames(out$vOmega) = ""
# Add matrix of conditional low-level class proportions
out$mPi = matrix(output$cPi,nrow(output$cPi),iT*iM)
colnames(out$mPi) = paste0(rep(paste0("P(C",1:iT,"|"), iM),rep(paste0("G",1:iM,")"),rep(iT,iM)))
# Add vector of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi_avg)))
if(dataout){
data = cbind(mY,id_high,mZ)
rownames(data) = NULL
colnames(data) = c(Y,id_high_name,Z)
for(i in 1:length(id_high_levs)){
data[data[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
data = data[,-which(colnames(data)=="Intercept"),drop=FALSE]
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX, dim(output$cLogPMX), dimnames=list(NULL,paste0("C",1:iT,",G"),colnames(out$mPi_avg)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX, dim(output$cPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi_avg)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C",1:iT,"|G"),colnames(out$mPi_avg)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi_avg)
if(dataout){
out$mSumPX = cbind(data,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi_avg)
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi_avg)
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi_avg)
if(dataout){
out$mPW_N = cbind(data,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$vAlpha = output$vDelta
rownames(out$vAlpha) = paste0("alpha(G",2:iM,")")
colnames(out$vAlpha) = ""
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),paste0("G",1:iM)))
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(rownames(out$vAlpha),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:(iM-1+P*(iT-1)*iM),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:(iM-1+P*(iT-1)*iM),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[iM:(iM-1+(P*(iT-1)*iM))],c(iT-1,P,iM)),3,t),dim(out$cGamma), dimnames=dimnames(out$cGamma))
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mQ) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat[1:(iM-1+P*(iT-1)*iM),1:(iM-1+P*(iT-1)*iM)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$mVar_corr
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$mV2)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = rownames(out$parvec)
# Add cube of Fisher information matrix for gamma
out$cGamma_Info = array(output$cGamma_Info,dim(output$cGamma_Info),dimnames=list(paste0("gamma(",Z,"|C,G)"),paste0("gamma(",Z,"|C,G)"),paste0("C", 2:iT,rep(paste0(",G", 1:iM),rep(iT-1,iM)))))
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = c(rownames(out$vAlpha),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"))
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add subset of matrix of model parameter contributions to log-likelihood score for gamma
out$mGamma_Score = output$mGamma_Score
colnames(out$mGamma_Score) = paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")")
# Add specification
out$spec = as.matrix("Multilevel LC model with lower-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
clean_output5 = function(output,Y,iT,iM,Z,Zh,mY,mZ,mZh,id_high,iH,P,P_high,id_high_levs,id_high_name,extout,dataout,ivItemcat){
if(!extout){
# Create empty list for output
out = list()
# Add matrix of sample means of high-level class proportions
out$vOmega_avg = as.matrix(apply(output$mOmega,2,mean))
rownames(out$vOmega_avg) = paste0("P(G",1:iM,")")
colnames(out$vOmega_avg) = ""
# Add matrix of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$mAlpha = t(output$mDelta)
rownames(out$mAlpha) = paste0("alpha(",Zh,"|G)")
colnames(out$mAlpha) = paste0("G",2:iM)
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add corrected standard errors for alpha
out$SEs_cor_alpha = t(matrix(output$SEs_cor[1:(P_high*(iM-1))],iM-1,P_high))
rownames(out$SEs_cor_alpha) = rownames(out$mAlpha)
colnames(out$SEs_cor_alpha) = colnames(out$mAlpha)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[(1+P_high*(iM-1)):(P_high*(iM-1)+P*(iT-1)*iM)],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add number of iterations
out$iter = output$iter
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add specification
out$spec = as.matrix("Multilevel LC model with lower- and higher-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
} else{
Ylab_beta = unlist(lapply(ivItemcat,function(x){if(x==2){TRUE}else{c(FALSE,rep(TRUE,x-1))}}))
# Create empty list for output
out = list()
# Add matrix of high-level class proportions
out$mOmega = output$mOmega
colnames(out$mOmega) = paste0("P(G",1:iM,")")
# Add matrix of sample means of high-level class proportions
out$vOmega_avg = as.matrix(apply(output$mOmega,2,mean))
rownames(out$vOmega_avg) = colnames(out$mOmega)
colnames(out$vOmega_avg) = ""
# Add matrix of conditional low-level class proportions
out$mPi = matrix(output$cPi,nrow(output$cPi),iT*iM)
colnames(out$mPi) = paste0(rep(paste0("P(C",1:iT,"|"),iM),rep(paste0("G",1:iM,")"),rep(iT,iM)))
# Add matrix of sample means of conditional low-level class proportions
out$mPi_avg = apply(output$cPi,3,function(x){apply(x,2,mean)})
rownames(out$mPi_avg) = paste0("P(C",1:iT,"|G)")
colnames(out$mPi_avg) = paste0("G",1:iM)
# Add matrix of conditional response probabilities
out$mPhi = output$mPhi
rownames(out$mPhi) = paste0("P(",Y,"|C)")
colnames(out$mPhi) = paste0("C",1:iT)
# Add cube of joint posterior class membership probabilities
out$cPMX = array(output$cPMX,dim(output$cPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),paste0("G", 1:iM)))
if(dataout){
data_low = cbind(mY,id_high,mZ)
rownames(data_low) = NULL
colnames(data_low) = c(Y,id_high_name,Z)
for(i in 1:length(id_high_levs)){
data_low[data_low[,id_high_name]==i,id_high_name] = id_high_levs[i]
}
data_low = data_low[,-which(colnames(data_low)=="Intercept"),drop=FALSE]
out$cPMX = array(unlist(lapply(1:dim(out$cPMX)[3],function(x){cbind(data_low,out$cPMX[,,x])})),
c(dim(out$cPMX)[1],ncol(data_low)+ncol(out$cPMX),dim(out$cPMX)[3]),
dimnames=list(dimnames(out$cPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPMX)[[2]]),dimnames(out$cPMX)[[3]]))
}
# Add cube of log of joint posterior class membership probabilities
out$cLogPMX = array(output$cLogPMX,dim(output$cLogPMX),dimnames=list(NULL,paste0("C",1:iT,",G"),paste0("G", 1:iM)))
if(dataout){
out$cLogPMX = array(unlist(lapply(1:dim(out$cLogPMX)[3],function(x){cbind(data_low,out$cLogPMX[,,x])})),
c(dim(out$cLogPMX)[1],ncol(data_low)+ncol(out$cLogPMX),dim(out$cLogPMX)[3]),
dimnames=list(dimnames(out$cLogPMX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPMX)[[2]]),dimnames(out$cLogPMX)[[3]]))
}
# Add cube of conditional posterior low-level class membership probabilities
out$cPX = array(output$cPX,dim(output$cPX),dimnames=list(NULL,paste0("C", 1:iT,"|G"),paste0("G", 1:iM)))
if(dataout){
out$cPX = array(unlist(lapply(1:dim(out$cPX)[3],function(x){cbind(data_low,out$cPX[,,x])})),
c(dim(out$cPX)[1],ncol(data_low)+ncol(out$cPX),dim(out$cPX)[3]),
dimnames=list(dimnames(out$cPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cPX)[[2]]),dimnames(out$cPX)[[3]]))
}
# Add cube of log of conditional posterior low-level class membership probabilities
out$cLogPX = array(output$cLogPX,dim(output$cLogPX),dimnames=list(NULL,paste0("C", 1:iT,"|G"),paste0("G", 1:iM)))
if(dataout){
out$cLogPX = array(unlist(lapply(1:dim(out$cLogPX)[3],function(x){cbind(data_low,out$cLogPX[,,x])})),
c(dim(out$cLogPX)[1],ncol(data_low)+ncol(out$cLogPX),dim(out$cLogPX)[3]),
dimnames=list(dimnames(out$cLogPX)[[1]],c(Y,id_high_name,Z[-1],dimnames(out$cLogPX)[[2]]),dimnames(out$cLogPX)[[3]]))
}
# Add matrix of posterior high-level class membership probabilities for low-level units after marginalizing over low-level classes
out$mSumPX = output$mSumPX
colnames(out$mSumPX) = colnames(out$mPi_avg)
if(dataout){
out$mSumPX = cbind(data_low,out$mSumPX)
rownames(out$mSumPX) = NULL
}
# Add matrix of posterior high-level class membership probabilities for high-level units
out$mPW = output$mPW
rownames(out$mPW) = id_high_levs
colnames(out$mPW) = colnames(out$mPi_avg)
if(dataout){
data_high = mZh
rownames(data_high) = id_high_levs
colnames(data_high) = Zh
data_high = data_high[,-which(colnames(data_high)=="Intercept"),drop=FALSE]
out$mPW = cbind(data_high,out$mPW)
}
# Add matrix of log of posterior high-level class membership probabilities for high-level units
out$mlogPW = output$mlogPW
rownames(out$mlogPW) = id_high_levs
colnames(out$mlogPW) = colnames(out$mPi_avg)
if(dataout){
out$mlogPW = cbind(data_high,out$mlogPW)
}
# Add matrix of posterior high-level class membership probabilities for low-level units
out$mPW_N = output$mPW_N
colnames(out$mPW_N) = colnames(out$mPi_avg)
if(dataout){
out$mPW_N = cbind(data_low,out$mPW_N)
rownames(out$mPW_N) = NULL
}
# Add matrix of posterior low-level class membership probabilities for low-level units after marginalizing over high-level classes
out$mPMsumX = output$mPMsumX
colnames(out$mPMsumX) = colnames(out$mPhi)
if(dataout){
out$mPMsumX = cbind(data_low,out$mPMsumX)
rownames(out$mPMsumX) = NULL
}
# Add low-level entropy R-sqr
out$R2entr_low = output$R2entr_low
# Add high-level entropy R-sqr
out$R2entr_high = output$R2entr_high
# Add low-level Bayesian information criterion
out$BIClow = output$BIClow
# Add high-level Bayesian information criterion
out$BIChigh = output$BIChigh
# Add low-level integrated completed likelihood Bayesian information criterion
out$ICL_BIClow = output$ICL_BIClow
# Add high-level integrated completed likelihood Bayesian information criterion
out$ICL_BIChigh = output$ICL_BIChigh
# Add Akaike information criterion
out$AIC = output$AIC
# Add alphas
out$mAlpha = t(output$mDelta)
rownames(out$mAlpha) = paste0("alpha(",Zh,"|G)")
colnames(out$mAlpha) = paste0("G",2:iM)
# Add gammas
out$cGamma = array(apply(output$cGamma,3,t),c(P,iT-1,iM),dimnames=list(paste0("gamma(",Z,"|C)"),paste0("C",2:iT,",G"),colnames(out$mPi_avg)))
# Add betas
out$mBeta = output$mBeta
rownames(out$mBeta) = paste0("beta(",Y[Ylab_beta],"|C)")
colnames(out$mBeta) = colnames(out$mPhi)
# Add vector of model parameters
out$parvec = output$parvec
rownames(out$parvec) = c(paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")"),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"),paste0(rep(substr(rownames(out$mBeta),1,nchar(rownames(out$mBeta))-2),iT),rep(colnames(out$mBeta),rep(nrow(out$mBeta),iT)),")"))
colnames(out$parvec) = ""
# Add vector of uncorrected standard errors
out$SEs_unc = as.matrix(output$SEs_unc[1:(P_high*(iM-1)+P*(iT-1)*iM),])
rownames(out$SEs_unc) = rownames(out$mV2)
colnames(out$SEs_unc) = ""
# Add vector of corrected standard errors
out$SEs_cor = as.matrix(output$SEs_cor[1:(P_high*(iM-1)+P*(iT-1)*iM),])
rownames(out$SEs_cor) = rownames(out$mV2)
colnames(out$SEs_cor) = ""
# Add corrected standard errors for alpha
out$SEs_cor_alpha = t(matrix(output$SEs_cor[1:(P_high*(iM-1))],iM-1,P_high))
rownames(out$SEs_cor_alpha) = rownames(out$mAlpha)
colnames(out$SEs_cor_alpha) = colnames(out$mAlpha)
# Add corrected standard errors for gamma
out$SEs_cor_gamma = array(apply(array(output$SEs_cor[(1+P_high*(iM-1)):(P_high*(iM-1)+P*(iT-1)*iM)],c(iT-1,P,iM)),3,t),dim(out$cGamma),dimnames=dimnames(out$cGamma))
# Add mQ
out$mQ = output$mQ
rownames(out$mQ) = colnames(out$mV2) = rownames(out$mV2)
# Add uncorrected variance-covariance matrix
out$Varmat_unc = output$Varmat[1:(P_high*(iM-1)+P*(iT-1)*iM),1:(P_high*(iM-1)+P*(iT-1)*iM)]
rownames(out$Varmat_unc) = colnames(out$Varmat_unc) = rownames(out$mV2)
# Add corrected variance-covariance matrix
out$Varmat_cor = output$mVar_corr
rownames(out$Varmat_cor) = colnames(out$Varmat_cor) = rownames(out$Varmat_unc)
# Add Fisher information matrix
out$Infomat = output$Infomat
rownames(out$Infomat) = colnames(out$Infomat) = rownames(out$parvec)
# Add cube of Fisher information matrix for alpha
out$cAlpha_Info = array(output$cDelta_Info, dim(output$cDelta_Info), dimnames=list(paste0("alpha(",Zh,"|G)"),paste0("alpha(",Zh,"|G)"),colnames(out$mAlpha)))
# Add cube of Fisher information matrix for gamma
out$cGamma_Info = array(output$cGamma_Info, dim(output$cGamma_Info), dimnames=list(paste0("gamma(",Z,"|C,G)"),paste0("gamma(",Z,"|C,G)"),paste0("C",2:iT,rep(paste0(",G",1:iM),rep(iT-1,iM)))))
# Add inverse of the information matrix from the second step
out$mV2 = output$mV2
rownames(out$mV2) = colnames(out$mV2) = c(paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")"),paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")"))
# Add number of iterations
out$iter = output$iter
# Add epsilon
out$eps = output$eps
# Add log-likelihood series
out$LLKSeries = output$LLKSeries
# Add current log-likelihood for high-level units
out$vLLK = output$vLLK
rownames(out$vLLK) = id_high_levs
colnames(out$vLLK) = ""
# Add matrix of model parameter contributions to log-likelihood score
out$mScore = output$mScore
colnames(out$mScore) = rownames(out$parvec)
# Add subset of matrix of model parameter contributions to log-likelihood score for alpha
out$mAlpha_Score = output$mDelta_Score
rownames(out$mAlpha_Score) = id_high_levs
colnames(out$mAlpha_Score) = paste0(rep(substr(rownames(out$mAlpha),1,nchar(rownames(out$mAlpha))-2),rep(iM-1,P_high)),rep(colnames(out$mAlpha),P_high),")")
# Add subset of matrix of model parameter contributions to log-likelihood score for gamma
out$mGamma_Score = output$mGamma_Score
colnames(out$mGamma_Score) = paste0(apply(out$cGamma,3,function(x){paste0(rep(substr(rownames(x),1,nchar(rownames(x))-2),rep(iT-1,P)),rep(substr(colnames(x),1,nchar(colnames(x))-2),P),",")}),rep(unlist(dimnames(out$cGamma)[3]),rep(P*(iT-1),iM)),")")
# Add specification
out$spec = as.matrix("Multilevel LC model with lower- and higher-level covariates")
rownames(out$spec) = colnames(out$spec) = ""
}
return(out)
}
#
clean_cov = function(data,covnam){
cov = data[,covnam,drop=FALSE]
if(any(!sapply(cov,is.numeric))){
covnam_upd = c()
for(i in covnam){
if(!sapply(cov,is.numeric)[which(colnames(cov)==i)]){
covnam_lev = sort(unique(cov[,i]))
cov_missing = is.na(as.numeric(factor(cov[,i])))
cov_upd = vecTomatClass(as.numeric(factor(cov[,i])))
cov_upd[cov_missing,] = NA
colnames(cov_upd) = paste0(i,".",covnam_lev)
cov_upd = cov_upd[,-1,drop=FALSE]
covnam_upd = c(covnam_upd,colnames(cov_upd))
} else{
cov_upd = cov[,i,drop=FALSE]
covnam_upd = c(covnam_upd,i)
}
cov = cov[,-which(colnames(cov)==i),drop=FALSE]
cov = cbind(cov,cov_upd)
}
covnam = covnam_upd
colnames(cov) = covnam
}
cov = as.matrix(cbind(1,cov))
return(list(cov=cov,covnam=covnam))
}
#
update_YmY = function(mY,Y,ivItemcat){
Y_upd = c()
for(i in Y){
if(ivItemcat[i]>2){
mY_upd = vecTomatClass(mY[,i]+1)
colnames(mY_upd) = paste0(i,".",1:ncol(mY_upd)-1)
Y_upd = c(Y_upd,colnames(mY_upd))
} else{
mY_upd = mY[,i,drop=FALSE]
Y_upd = c(Y_upd,i)
}
mY = mY[,-which(colnames(mY)==i)]
mY = cbind(mY,mY_upd)
}
Y = Y_upd
colnames(mY) = Y
return(list(mY=mY,Y=Y))
}
#
update_YmY_nonnum = function(mY,Y,ivItemcat){
itemnam = lapply(data.frame(mY),function(x){levels(factor(x))})
mY = apply(mY,2,function(x){as.numeric(factor(x))-1})
Y_upd = c()
for(i in Y){
if(ivItemcat[i]>2){
mY_upd = vecTomatClass(mY[,i]+1)
colnames(mY_upd) = paste0(i,".",itemnam[[i]])
Y_upd = c(Y_upd,colnames(mY_upd))
} else{
mY_upd = mY[,i,drop=FALSE]
Y_upd = c(Y_upd,paste0(i,".",itemnam[[i]][2]))
}
mY = mY[,-which(colnames(mY)==i)]
mY = cbind(mY,mY_upd)
}
Y = Y_upd
colnames(mY) = Y
return(list(mY=mY,Y=Y))
}
#
check_inputs1 = function(data,Y,iT,id_high,iM,Z,Zh){
# <data> is matrix or dataframe
if(!(is.matrix(data)|is.data.frame(data))){
stop("data must be matrix or dataframe.",call.=FALSE)
}
# <Y>, <id_high>, <Z> and <Zh> are column names in <data>
if(any(!c(Y,id_high,Z,Zh)%in%colnames(data))){
if(any(!Y%in%colnames(data))){
stop("Not all Y in data.",call.=FALSE)
}
if(!is.null(id_high)&any(!id_high%in%colnames(data))){
stop("id_high not in data.",call.=FALSE)
}
if(!is.null(Z)&any(!Z%in%colnames(data))){
stop("Not all Z in data.",call.=FALSE)
}
if(!is.null(Zh)&any(!Zh%in%colnames(data))){
stop("Not all Zh in data.",call.=FALSE)
}
}
}
#
check_inputs2 = function(data,Y,iT,id_high,iM,Z,Zh){
# Single <id_high>
if(!is.null(id_high)){
if(length(id_high)!=1){
stop("Invalid id_high.",call.=FALSE)
}
}
# <iT> positive integer or sequential positive integers
if(!length(iT)>=1){
stop("Invalid iT.",call.=FALSE)
} else{
if(length(iT)==1){
if(!is.numeric(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(iT!=round(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(iT==1){
stop("Invalid iT: not a mixture model.",call.=FALSE)
} else if(iT<1){
stop("Invalid iT.",call.=FALSE)
}
} else{
if(!is.vector(iT)|is.list(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(!is.numeric(iT)){
stop("Invalid iT.",call.=FALSE)
} else if(any(!iT==round(iT))){
stop("Invalid iT.",call.=FALSE)
} else if(min(iT)<1){
stop("Invalid iT.",call.=FALSE)
} else if(length(iT)!=length(min(iT):max(iT))){
stop("Invalid iT.",call.=FALSE)
} else if(any(!iT==min(iT):max(iT))){
stop("Invalid iT.",call.=FALSE)
}
}
}
# <iM> positive integer or sequential positive integers
if(!is.null(iM)&!length(iM)>=1){
stop("Invalid iM.",call.=FALSE)
} else if(!is.null(iM)){
if(length(iM)==1){
if(!is.numeric(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(iM!=round(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(iM<1){
stop("Invalid iM.",call.=FALSE)
}
} else{
if(!is.vector(iM)|is.list(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(!is.numeric(iM)){
stop("Invalid iM.",call.=FALSE)
} else if(any(!iM==round(iM))){
stop("Invalid iM.",call.=FALSE)
} else if(min(iM)<1){
stop("Invalid iM.",call.=FALSE)
} else if(length(iM)!=length(min(iM):max(iM))){
stop("Invalid iM.",call.=FALSE)
} else if(any(!iM==min(iM):max(iM))){
stop("Invalid iM.",call.=FALSE)
}
}
}
# <Y> consecutive integers from 0 or be non-numeric character
if(any(sapply(as.data.frame(data[,Y]),function(x){(any(sort(unique(x))!=0:(length(unique(x))-1)))&!is.character(x)}))){
stop("Items must contain consecutive integers from 0 or be non-numeric character.",call.=FALSE)
}
# No missings in <id_high>
if(any(is.na(data[,id_high]))){
stop("id_high cannot contain missing values.",call.=FALSE)
}
# <id_high> not labelled "Intercept"
if(!is.null(id_high)){
if(id_high=="Intercept"){
stop("id_high may not be labelled 'Intercept'.",call.=FALSE)
}
}
# No duplicate <Y>, <Z> or >Zh>
if(any(duplicated(as.list(as.data.frame(data[,Y]))))){
stop("Duplicate items.",call.=FALSE)
}
if(!is.null(Z)){
if(length(Z)>2){
if(any(duplicated(as.list(data[,Z])))){
stop("Duplicate lower-level covariates.",call.=FALSE)
}
}
}
if(!is.null(Zh)){
if(length(Zh)>2){
if(any(duplicated(as.list(data[,Zh])))){
stop("Duplicate higher-level covariates.",call.=FALSE)
}
}
}
# No constant <Y>, <Z> or <Zh>
if(any(apply(data[,Y],2,function(x){length(unique(x))==1}))){
stop("Constant in indicators.",call.=FALSE)
}
if(!is.null(Z)){
if(any(apply(data[,Z,drop=FALSE],2,function(x){length(unique(x))==1}))){
stop("Constant in lower-level covariates.",call.=FALSE)
}
}
if(!is.null(Zh)){
if(any(apply(data[,Zh,drop=FALSE],2,function(x){length(unique(x))==1}))){
stop("Constant in higher-level covariates.",call.=FALSE)
}
}
# Identify specification
if(is.null(id_high)&is.null(iM)){
specification = "single-level"
} else if(!is.null(id_high)&!is.null(iM)){
specification = "multilevel"
} else{
if(is.null(id_high)){
stop("iM specified, id_high not.",call.=FALSE)
}
if(is.null(iM)){
stop("id_high specified, iM not.",call.=FALSE)
}
}
if(is.null(Z)){
specification = paste(specification,"measurement")
if(!is.null(Zh)){
stop("Zh specified, Z not.",call.=FALSE)
}
} else if(!is.null(Z)){
specification = paste(specification,"structural")
if(specification=="single-level structural"&!is.null(Zh)){
stop("Zh specified, iM and id_high not.",call.=FALSE)
}
}
# Identify approach
if(length(iT)==1&is.null(iM)){
approach = "direct"
} else if(length(iT)==1&length(iM)==1){
approach = "direct"
} else if(length(iT)>1&is.null(iM)){
approach = "model selection on low"
} else if(length(iT)>1&length(iM)==1){
if(iM==1){
approach = "model selection on low"
} else{
approach = "model selection on low with high"
}
} else if(length(iT)==1&length(iM)>1){
approach = "model selection on high"
} else if(length(iT)>1&length(iM)>1){
approach = "model selection on low and high"
}
return(approach)
}
#
#
#
print.multiLCA = function(x,...){
out = NULL
if(x$spec == "Single-level LC model"){
######################################
### Single-level measurement model ###
######################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nCLASS PROPORTIONS:\n")
print(round(x$vPi, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$AvgClassErrProb, 4)), as.character(round(x$R2entr, 4)), as.character(round(x$BIC, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("ClassErr", "EntR-sqr", "BIC", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
} else if(x$spec == "Single-level LC model with covariates"){
#####################################
### Single-level structural model ###
#####################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n")
print(round(x$vPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$AvgClassErrProb, 4)), as.character(round(x$R2entr, 4)), as.character(round(x$BIC, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("ClassErr", "EntR-sqr", "BIC", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR CLASS MEMBERSHIP:\n\n")
gammas = format(round(as.matrix(x$cGamma), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma)/as.matrix(x$SEs_cor_gamma), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma)/as.matrix(x$SEs_cor_gamma)))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (i in 1:ncol(gammas)){
C = paste0("C", 1 + i)
logit_params = noquote(cbind(gammas[,i], SE[,i], Zscore[,i], matrix(paste0(pval[,i], psignf[,i]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma)), 1, nchar(rownames(as.matrix(x$cGamma))) - 1), 1 + i, ")")
cat("\nMODEL FOR", C, "(BASE C1)\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
} else if(x$spec == "Multilevel LC model"){
####################################
### Multilevel measurement model ###
####################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS:\n")
print(round(x$vOmega, 4))
cat("\nCLASS PROPORTIONS:\n")
print(round(x$mPi, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
} else if(x$spec == "Multilevel LC model with lower-level covariates"){
###############################################################
### Multilevel structural model with lower-level covariates ###
###############################################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS:\n")
print(round(x$vOmega, 4))
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n\n")
print(round(x$mPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR LOWER-LEVEL CLASS MEMBERSHIP:\n\n")
for (i in 1:dim(x$cGamma)[3]){
gammas = format(round(as.matrix(x$cGamma[,,i]), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i])))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (j in 1:ncol(gammas)){
G = paste0("G", i)
C = paste0("C", 1 + j)
logit_params = noquote(cbind(gammas[,j], SE[,j], Zscore[,j], matrix(paste0(pval[,j], psignf[,j]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma[,,i])), 1, nchar(rownames(as.matrix(x$cGamma[,,i]))) - 1), j + 1, ",G", i, ")")
cat("\nMODEL FOR", C, "(BASE C1) GIVEN", G, "\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
}
} else if(x$spec == "Multilevel LC model with lower- and higher-level covariates"){
#########################################################################
### Multilevel structural model with low- and higher-level covariates ###
#########################################################################
cat("\nCALL:\n")
print(x$call)
cat("\nSPECIFICATION:\n")
print(noquote(x$spec))
iter = t(matrix(c(x$iter, x$LLKSeries[1], tail(x$LLKSeries, 1))))
colnames(iter) = c("EMiter", "LLfirst", "LLlast")
rownames(iter) = ""
cat("\nESTIMATION DETAILS:\n\n")
print(iter)
cat("\n---------------------------\n")
cat("\nGROUP PROPORTIONS (SAMPLE MEAN):\n")
print(round(x$vOmega_avg, 4))
cat("\nCLASS PROPORTIONS (SAMPLE MEAN):\n\n")
print(round(x$mPi_avg, 4))
cat("\nRESPONSE PROBABILITIES:\n\n")
print(round(x$mPhi, 4))
cat("\n---------------------------\n")
cat("\nMODEL AND CLASSIFICATION STATISTICS:\n")
stat = as.matrix(c(as.character(round(x$R2entr_low, 4)), as.character(round(x$R2entr_high, 4)), as.character(round(x$BIClow, 4)), as.character(round(x$BIChigh, 4)), as.character(round(x$ICL_BIClow, 4)), as.character(round(x$ICL_BIChigh, 4)), as.character(round(x$AIC, 4))))
rownames(stat) = c("R2entrlow", "R2entrhigh", "BIClow", "BIChigh", "ICLBIClow", "ICLBIChigh", "AIC")
colnames(stat) = ""
print(noquote(stat), right = TRUE)
cat("\n")
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR HIGHER-LEVEL CLASS MEMBERSHIP:\n\n")
alphas = format(round(as.matrix(x$mAlpha), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_alpha), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$mAlpha)/as.matrix(x$SEs_cor_alpha), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$mAlpha)/as.matrix(x$SEs_cor_alpha)))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(alphas), ncol(alphas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (i in 1:ncol(alphas)){
G = paste0("G", 1 + i)
logit_params = noquote(cbind(alphas[,i], SE[,i], Zscore[,i], matrix(paste0(pval[,i], psignf[,i]))))
colnames(logit_params) = c("Alpha", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$mAlpha)), 1, nchar(rownames(as.matrix(x$mAlpha))) - 1), 1 + i, ")")
cat("\nMODEL FOR", G, "(BASE G1)\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
cat("\n---------------------------\n")
cat("\nLOGISTIC MODEL FOR LOWER-LEVEL CLASS MEMBERSHIP:\n\n")
for (i in 1:dim(x$cGamma)[3]){
gammas = format(round(as.matrix(x$cGamma[,,i]), 4), nsmall = 4, scientific = FALSE)
SE = format(round(as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
Zscore = format(round(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i]), 4), nsmall = 4, scientific = FALSE)
pval = format(round(2*(1 - pnorm(abs(as.matrix(x$cGamma[,,i])/as.matrix(x$SEs_cor_gamma[,,i])))), 4), nsmall = 4, scientific = FALSE)
psignf = matrix(" ", nrow(gammas), ncol(gammas))
psignf[pval < 0.1 & pval >= 0.05] = "* "
psignf[pval < 0.05 & pval >= 0.01] = "** "
psignf[pval < 0.01] = "***"
for (j in 1:ncol(gammas)){
G = paste0("G", i)
C = paste0("C", 1 + j)
logit_params = noquote(cbind(gammas[,j], SE[,j], Zscore[,j], matrix(paste0(pval[,j], psignf[,j]))))
colnames(logit_params) = c("Gamma", "S.E.", "Z-score", "p-value")
rownames(logit_params) = paste0(substr(rownames(as.matrix(x$cGamma[,,i])), 1, nchar(rownames(as.matrix(x$cGamma[,,i]))) - 1), j + 1, ",G", i, ")")
cat("\nMODEL FOR", C, "(BASE C1) GIVEN", G, "\n\n")
print(logit_params, right = TRUE)
cat("\n")
cat("\n", "*** p < 0.01, ** p < 0.05, * p < 0.1")
cat("\n")
}
}
}
}
#
plot.multiLCA = function(x, horiz = FALSE, clab = NULL, ...){
out = NULL
iT = ncol(x$mPhi)
items = 1:nrow(x$mPhi)
itemnames = substr(rownames(x$mPhi), 3, nchar(rownames(x$mPhi)) - 3)
if(horiz == TRUE){
las = 1
} else{
las = 3
}
if(is.null(clab)){
legend = paste0("Class ", 1:iT)
} else{
if(!is.vector(clab) | is.list(clab)){
stop("Invalid clab.",call.=FALSE)
} else if(length(clab) != iT | !is.character(clab)){
stop("Invalid clab.",call.=FALSE)
}
legend = clab
}
oldpar = par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar = c(5, 5, 3, 8), xpd = TRUE)
plot(items, x$mPhi[,1], type = "b", xaxt = "n",yaxt = "n", pch = 1,
xlab = "", ylab = "Response probability", frame.plot = FALSE, ylim = c(0,1), lty = 1, ...)
axis(1, at = items, labels = itemnames, las = las)
axis(2, las = 1)
for(h in 2:iT){
lines(items, x$mPhi[,h], pch = h, type = "b", lty = h)
}
legend("topright", legend = legend, inset = c(-0.4,0),
lty = (1:iT), pch = 1:iT, cex = 0.8)
}
Try the multilevLCA package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.