Nothing
glca_encode <- function(
call, terms, mf,
nclass, ncluster,
measure.inv, coeff.inv,
na.rm, verbose
)
{
# Import data
Y <- stats::model.response(mf)
Yattr <- attributes(Y)
if (is.null(dim(Y)))
dim(Y) <- c(length(Y), 1L)
if (!inherits(Y, "items"))
stop("Manifest items should be indicated by item function.\n")
grp <- stats::model.extract(mf, "group")
if (is.null(grp)) grp <- factor(rep("ALL", nrow(Y)))
else grp <- droplevels(factor(grp))
dataN <- Yattr$dataN
modelN <- nrow(mf)
isna <- which(rowSums(Y == 0L) > 0L)
tmis <- which(rowSums(Y != 0L) == 0L)
if (na.rm & length(isna) > 0L) {
Y <- Y[-isna, , drop = FALSE]
mf <- mf[-isna, , drop = FALSE]
grp <- droplevels(grp[-isna])
} else if (length(tmis) > 0L) {
Y <- Y[-tmis, , drop = FALSE]
mf <- mf[-tmis, , drop = FALSE]
grp <- droplevels(grp[-tmis])
}
lbs <- labels(terms)
if (length(lbs) > 0L) cvs <- stats::get_all_vars(stats::reformulate(lbs), mf)
else cvs <- NULL
# Model / variable name specification
N <- nrow(Y)
G <- nlevels(grp)
M <- ncol(Y)
R <- sapply(Yattr$y.level, length)
if (nclass < 2L)
stop("Number of latent classes should be greater than 1.")
if (!is.null(ncluster)) {
W <- ncluster
if (W < 2L) {
if (verbose)
cat("Number of latent clusters should be greater than 1.\n",
"\"ncluster = ", ncluster, "\" has been ignored.\n\n")
W <- 0L
} else if (G <= ncluster) {
if (verbose)
cat("Number of latent clusters should be less than number of groups.\n",
"\"ncluster = ", ncluster, "\" has been ignored.\n\n")
W <- 0L
} else {
W <- floor(ncluster)
if (measure.inv == FALSE && verbose)
cat("For nonparametric multilevel latent class analysis, always \"Meas.inv = TRUE\"")
measure.inv <- TRUE
}
} else W <- 0L
C <- floor(nclass)
y.names <- Yattr$y.names
r.names <- Yattr$y.level
resp.name <- data.frame(
matrix("", M, max(R)), row.names = y.names,
stringsAsFactors = FALSE
)
names(resp.name) <- paste0("Y = ", seq_len(max(R)))
for (m in seq_len(M))
resp.name[m, seq_len(R[m])] <- c(r.names[[m]])
if (!is.null(cvs)) {
cvs <- stats::get_all_vars(stats::reformulate(lbs), mf)
Zind <- apply(cvs, 2L, function(x)
all(by(x, grp, function(y) length(unique(y)) == 1L)))
if (sum(!Zind) == 0) X <- stats::model.matrix(~ 1, mf)
else X <- stats::model.matrix(stats::reformulate(lbs[!Zind]), mf)
if (any(Zind))
Z <- stats::model.matrix(stats::reformulate(lbs[ Zind]), mf)[, -1L, drop = FALSE]
else if (W > 1) Z <- mf[FALSE]
else Z <- NULL
} else {
X <- stats::model.matrix(terms, mf)
Z <- if (W > 1) Z <- mf[FALSE] else NULL
}
P <- ncol(X)
Q <- if (!is.null(Z)) ncol(Z) else 0L
# Covariate names
X.names <- if (!is.null(cvs)) lbs[!Zind] else NULL
Z.names <- if (!is.null(cvs)) lbs[ Zind] else NULL
x.names <- colnames(X)
z.names <- colnames(Z)
g.names <- levels(grp)
# Grouping data
grp <- as.numeric(grp)
Ng <- sapply(seq_len(G), function(g) sum(grp == g))
y <- lapply(seq_len(G), function(g) as.matrix(Y[grp == g, , drop = FALSE]))
x <- lapply(seq_len(G), function(g) as.matrix(X[grp == g, , drop = FALSE]))
if (W > 1L)
z <- lapply(seq_len(G), function(g) as.matrix(Z[grp == g, , drop = FALSE]))
else z <- NULL
# Variable print
if (verbose) {
cat("Manifest items :\n", y.names, "\n")
if (!is.null(call$group)) cat("Grouping variable :", call$group, "\n")
if (Q > 0L) cat("Covariates (Level 2) : \n", Z.names, "\n")
if (P > 1L) cat("Covariates (Level 1) : \n", X.names, "\n")
cat("\nDeleted observation(s) : \n")
if (na.rm)
cat(dataN - modelN + length(isna),
"observation(s) for missing at least 1 variable \n\n")
else {
cat(length(tmis), "observation(s) for missing all manifest items\n")
cat(dataN - modelN, "observation(s) for missing at least 1 covariates\n\n")
}
}
# Satuated model
grpx <- cbind(grp, X)
uniqH <- unique(grpx)
hind <- match(data.frame(t(grpx)), data.frame(t(uniqH)))
H <- nrow(uniqH)
fulldf <- prod(R) * H
if (prod(R) < 1e+6 & na.rm != TRUE & length(isna) != 0L) {
loglikh <- numeric(H)
for (h in seq_len(H)) {
Yh <- Y[hind == h, , drop = FALSE]
loglikh[h] <- ObsLik(as.matrix(Yh), nrow(Yh), M, R, 1000, 1e-8)
}
loglik0 <- sum(loglikh)
} else {
Y0 <- Y[rowSums(Y == 0L) == 0L, , drop = FALSE]
h0 <- hind[rowSums(Y == 0L) == 0L]
loglikh <- numeric(H)
for (h in seq_len(H)) {
Yh <- Y0[h0 == h, , drop = FALSE]
Y.sorted <- Yh[do.call(order, data.frame(Yh)[M:1]), , drop = FALSE]
pattern <- as.matrix(unique(Y.sorted))
obsvd <- ObsCell2(as.matrix(Y.sorted), pattern,
nrow(Y.sorted), nrow(pattern))
loglikh[h] <- sum(obsvd * log(obsvd / sum(obsvd)))
}
loglik0 <- sum(loglikh)
}
if (W == 0L) {
if (G == 1L)
type <- "Latent class analysis"
else
type <- "Multiple-group latent class analysis"
if (P == 1L) {
if (measure.inv)
npar <- G * (C - 1L) + C * sum(R - 1L)
else
npar <- G * (C - 1L) + G * C * sum(R - 1L)
} else {
if (coeff.inv)
npar <- G * (C - 1L) + (C - 1L) * (P - 1L)
else
npar <- G * (C - 1L) * P
if (measure.inv)
npar <- npar + C * sum(R - 1L)
else
npar <- npar + G * C * sum(R - 1L)
}
} else {
type <- "Nonparametric multilevel latent class analysis"
if (P == 1L && Q == 0)
npar <- W - 1L + W * (C - 1L) + C * sum(R - 1L)
else {
if (coeff.inv)
npar <- W - 1L + W * (C - 1L) + (P - 1L + Q) * (C - 1L) +
C * sum(R - 1L)
else
npar <- W - 1L + (W * P + Q) * (C - 1L) + C * sum(R - 1L)
}
}
return(
list(datalist = list(y = y, x = x, z = z, group = grp,
loglik0 = loglik0),
model = list(type = type,
measure.inv = measure.inv,
coeff.inv = coeff.inv,
N = N, Ng = Ng, G = G,
C = C, W = W, M = M, R = R,
P = P, Q = Q,
npar = npar,
df = min(N - 1L, fulldf - 1L) - npar),
vname = list(y.names = y.names,
g.names = g.names,
r.names = r.names,
X.names = X.names,
Z.names = Z.names,
x.names = x.names,
z.names = z.names,
resp.name = resp.name))
)
}
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.