Nothing
## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## ------------------------------------------------------------------------
alpha = rep(0.25, 4) ## equal
p1 = 0.85
p2 = 0.1
p3 = 0.2
P = matrix(c(
p1,p2,p2,p2,
p1,p3,p3,p3,
p2,p1,p2,p2,
p3,p1,p3,p3,
p2,p2,p1,p2,
p3,p3,p1,p3,
p2,p2,p2,p1,
p3,p3,p3,p1), nrow=8, ncol=4, byrow = TRUE)
## ------------------------------------------------------------------------
simLCA <- function(alpha, ProbMat, sampleSize, nsim) {
## Reformat the probability Matrix to a form required by poLCA.simdata
probs = vector("list", nrow(ProbMat))
for (i in 1:nrow(ProbMat)) {
probs[[i]] = cbind(ProbMat[i,], 1 - ProbMat[i,])
}
X = poLCA.simdata(N = sampleSize*nsim, probs = probs, P = alpha)$dat
split(X, rep(1:nsim, each=sampleSize))
}
## ------------------------------------------------------------------------
library(poLCA)
set.seed(37421)
X = simLCA(alpha, P, sampleSize=50, nsim=1)[[1]]
head(X)
## ------------------------------------------------------------------------
f = cbind(Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8) ~ 1
poLCA(f, X, nclass = 3, verbose=FALSE)
## ------------------------------------------------------------------------
library(sBIC)
lcas = LCAs(maxNumClasses=6, numVariables=8, numStatesForVariables=2)
results = sBIC(X, lcas)
## ------------------------------------------------------------------------
phivec = (6:10)/2
## ---- fig.cap="Influence of penalty parameter phi on sBIC for LCA"-------
plot(1:6, results$sBIC - max(results$sBIC), type="n", xlab="number of latent classes",
ylab=expression(bar(SBIC)[phi]))
for (i in seq_along(phivec)) {
setPhi(lcas, phivec[i])
results = sBIC(NULL, lcas)
lines(1:6, results$sBIC - max(results$sBIC), pch=i, lty=i, type="b")
}
legend("topleft", pch=1:6, lty=1:6, legend=paste("phi=", phivec))
## ------------------------------------------------------------------------
fitBIC <- function(X, maxNumClasses, phis) {
lcas = LCAs(maxNumClasses, numVariables=ncol(X), numStatesForVariables=2)
results = sBIC(X, lcas)
nclass.bic <- which.max(results$BIC)
nclass.sbic = numeric(length(phis))
for (k in seq_along(phis)) {
lcas$setPhi(phis[k])
results = sBIC(NULL, lcas)
nclass.sbic[k] = which.max(results$sBIC)
}
c(nclass.bic, nclass.sbic)
}
## ------------------------------------------------------------------------
createTable = function(Xlist, ...) {
## Set up cluster
cl <- makeCluster(detectCores() - 1)
clusterEvalQ(cl, library(sBIC))
clusterEvalQ(cl, library(poLCA))
clusterSetRNGStream(cl)
## Fit LCA models to simulated data and tabulate
bicResults = parSapply(cl, Xlist, fitBIC, ...)
bicTab = apply(bicResults, 1, tabulate, nbin=maxNumClasses)
dimnames(bicTab) = list(1:maxNumClasses, c("BIC", paste0("sBIC", phis)))
## End cluster
stopCluster(cl)
bicTab
}
## ---- eval=FALSE---------------------------------------------------------
# library(parallel)
# set.seed(1234)
# Xlist = simLCA(alpha, Pr, sampleSize=50, nsim=100)
# bictab = createTable(Xlist, maxNumClasses=6, phis=phivec)
# knitr::kable(bictab, row.names=TRUE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.