test.R

In kim0sun/catlvm: Models for Multiple Latent Classs Variables

debug(slcm)
{
   library(devtools)
   library(dplyr)
   load_all()
   lca   = slcm(L1[2] ~ X1 + X2 + X3)
   lcas  = slcm(L1[2] ~ X1 + X2 + X3,
                L2[2] ~ Y1 + Y2 + Y3)
   jlca  = slcm(L1[2] ~ X1 + X2 + X3,
                  L2[3] ~ Y1 + Y2 + Y3,
                  L3[3] ~ Z1 + Z2 + Z3,
                  JC[5] ~ L1 + L2 + L3)
   lcpa  = slcm(L1[2] ~ X1 + X2 + X3,
                  L2[2] ~ Y1 + Y2 + Y3,
                  L3[2] ~ Z1 + Z2 + Z3,
                  PF[2] ~ L1 + L2 + L3,
                  constraints = list(c("L1", "L2", "L3")))

   jlcpa = slcm(L1[3] ~ X11 + X21 + X31,
                  M1[3] ~ Y11 + Y21 + Y31,
                  N1[3] ~ Z11 + Z21 + Z31,
                  L2[3] ~ X12 + X22 + X32,
                  M2[3] ~ Y12 + Y22 + Y32,
                  N2[3] ~ Z12 + Z22 + Z32,
                  L3[3] ~ X13 + X23 + X33,
                  M3[3] ~ Y13 + Y23 + Y33,
                  N3[3] ~ Z13 + Z23 + Z33,
                  J1[3] ~ L1 + M1 + N1,
                  J2[3] ~ L2 + M2 + N2,
                  J3[3] ~ L3 + M3 + N3,
                  JP[3] ~ J1 + J2 + J3,
                  constraints = list(
                     c("L1", "L2", "L3"), c("M1", "M2", "M3"), c("N1", "N2", "N3"),
                     c("J1 ~ L1", "J2 ~ L2", "J3 ~ L3"),
                     c("J1 ~ M1", "J2 ~ M2", "J3 ~ M3"),
                     c("J1 ~ N1", "J2 ~ N2", "J3 ~ N3")))
   jlcpa1 = slcm(L1[3] ~ X11 + X21 + X31,
                M1[3] ~ Y11 + Y21 + Y31,
                N1[3] ~ Z11 + Z21 + Z31,
                L2[3] ~ X12 + X22 + X32,
                M2[3] ~ Y12 + Y22 + Y32,
                N2[3] ~ Z12 + Z22 + Z32,
                L3[3] ~ X13 + X23 + X33,
                M3[3] ~ Y13 + Y23 + Y33,
                N3[3] ~ Z13 + Z23 + Z33,
                J1[4] ~ L1 + M1 + N1,
                J2[4] ~ L2 + M2 + N2,
                J3[4] ~ L3 + M3 + N3,
                JP[3] ~ J1 + J2 + J3,
                constraints = list(
                   c("L1", "L2", "L3"), c("M1", "M2", "M3"), c("N1", "N2", "N3"),
                   c("J1 ~ L1", "J2 ~ L2", "J3 ~ L3"),
                   c("J1 ~ M1", "J2 ~ M2", "J3 ~ M3"),
                   c("J1 ~ N1", "J2 ~ N2", "J3 ~ N3")))
   jjcpa = slcm(L1[3] ~ X11 + X21 + X31,
                  J2[4] ~ L2 + M2 + N2,
                  J3[4] ~ L3 + M3 + N3,
                  M1[3] ~ Y11 + Y21 + Y31,
                  N1[3] ~ Z11 + Z21 + Z31,
                  L2[3] ~ X12 + X22 + X32,
                  M2[3] ~ Y12 + Y22 + Y32,
                  J1[5] ~ L1 + M1 + N1,
                  N2[3] ~ Z12 + Z22 + Z32,
                  L3[3] ~ X13 + X23 + X33,
                  M3[3] ~ Y13 + Y23 + Y33,
                  N3[3] ~ Z13 + Z23 + Z33,
                  JP[4] ~ J1 + J2 + J3)
   lta = slcm(L1[3] ~ X11 + X21 + X31,
                L2[3] ~ X12 + X22 + X32,
                L3[3] ~ X13 + X23 + X33,
                L1 ~ L2, L2 ~ L3,
                constraints = list(c("L1", "L2", "L3")))
   jlta = slcm(L1[3] ~ X11 + X21 + X31,
                 M1[3] ~ Y11 + Y21 + Y31,
                 N1[3] ~ Z11 + Z21 + Z31,
                 L2[3] ~ X12 + X22 + X32,
                 M2[3] ~ Y12 + Y22 + Y32,
                 N2[3] ~ Z12 + Z22 + Z32,
                 L3[3] ~ X13 + X23 + X33,
                 M3[3] ~ Y13 + Y23 + Y33,
                 N3[3] ~ Z13 + Z23 + Z33,
                 J1[3] ~ L1 + M1 + N1,
                 J2[3] ~ L2 + M2 + N2,
                 J3[3] ~ L3 + M3 + N3,
                 J1 ~ J2, J2 ~ J3,
                 constraints = list(
                    c("L1", "L2", "L3"), c("M1", "M2", "N3"), c("N1", "N2", "N3"),
                    c("J1 ~ L1", "J2 ~ L2", "J3 ~ L3"),
                    c("J1 ~ M1", "J2 ~ M2", "J3 ~ M3"),
                    c("J1 ~ N1", "J2 ~ N2", "J3 ~ N3")))
   mlcpa = slcm(L1[3] ~ X11 + X21 + X31,
                 M1[3] ~ Y11 + Y21 + Y31,
                 L2[3] ~ X12 + X22 + X32,
                 M2[3] ~ Y12 + Y22 + Y32,
                 L3[3] ~ X13 + X23 + X33,
                 M3[3] ~ Y13 + Y23 + Y33,
                 LP[3] ~ L1 + M1 + L2 + M2 + L3 + M3,
                 constraints = list(
                    c("L1", "L2", "L3"), c("M1", "M2", "M3"),
                    c("LP->L1", "LP->L2", "LP->L3"),
                    c("LP->M1", "LP->M2", "LP->M3")))
   lcawg = slcm(LG[2] ~ Z1 + Z2 + Z3,
                  LC[2] ~ X1 + X2 + X3,
                  LG ~ LC)
   lcpawg = slcm(LG[2] ~ Z1 + Z2 + Z3,
                   LG ~ P1,
                   L1[2] ~ X11 + X12 + X13,
                   L2[2] ~ X21 + X22 + X23,
                   L3[2] ~ X31 + X32 + X33,
                   P1[2] ~ L1 + L2 + L3,
                   constraints = list(c("L1", "L2", "L3")))
}

# lcas; plot(lcas)
# lca; plot(lca)
# lcpa; plot(lcpa)
# jlca; plot(jlca, abbreviation = TRUE)
# jlcpa; plot(jlcpa, abbreviation = TRUE)
# lta; plot(lta)
debug(fit)
library(tidyverse)
object <- jlcpa1
jlcpa1$args
sim <- object %>% simulate(1000, ncat = 2)
fit <- object %>% fit(data = sim$response, method = "em",
                              control = list(em.iterlim = 5000))
transitionPlot(t(matrix(colSums(exp(fit$joint[[2]])), 3, 3)))
transitionPlot(t(round(fit$estimates$param$tau$B*100)))
str(fit$score)
fit$joint[[1]]

fit %>% parameter_numbers()
undebug(estimate)
fff <- fit %>% fit(data = sim$response, restriction = c("rho1", "rho7"))
fff$estimates$se$param
fff$estimates$param
fit$estimates$param
logLik(fff)
fit$args$npar
fff$args$npar
colSums(fff$estimates$param$rho$a)

newdata <- sim$response[sample(1:1000, 300),]
fit_nlm <- object %>% fit(data = sim$response, method = "nlm", restriction = c("tau3"))
rapply(fit_hb$estimates$param, round, how = "list", digits = 3)
fit_hb <- object %>% fit(data = sim$response, method = "hybrid", restriction = c("tau3"))
fit_hb$estimates$param
fit_hb$estimates$se

library(randomLCA)
library(tidyverse)
data(symptoms)
dat = symptoms[rep(seq(symptoms$Freq), symptoms$Freq),]

lcpa_sym <- slcm(
   LC1[3] ~ Nightcough.13 + Wheeze.13 + Itchyrash.13 + FlexDerma.13,
   LC2[3] ~ Nightcough.45 + Wheeze.45 + Itchyrash.45 + FlexDerma.45,
   LC3[3] ~ Nightcough.6  + Wheeze.6  + Itchyrash.6  + FlexDerma.6,
   LC4[3] ~ Nightcough.7  + Wheeze.7  + Itchyrash.7  + FlexDerma.7,
   LCP[2] ~ LC1 + LC2 + LC3 + LC4,
   constraints = list(c("LC1", "LC2", "LC3", "LC4"))
)

lta_sym <- slcm(
   LC1[2] ~ Nightcough.13 + Wheeze.13 + Itchyrash.13 + FlexDerma.13,
   LC2[2] ~ Nightcough.45 + Wheeze.45 + Itchyrash.45 + FlexDerma.45,
   LC3[2] ~ Nightcough.6  + Wheeze.6  + Itchyrash.6  + FlexDerma.6,
   LC4[2] ~ Nightcough.7  + Wheeze.7  + Itchyrash.7  + FlexDerma.7,
   LC1 ~ LC2,
   LC2 ~ LC3,
   LC3 ~ LC4,
   constraints = list(c("LC1", "LC2", "LC3", "LC4"))
)

lcpa_sym = lcpa_sym %>% fit(dat)
object <- lcpa_sym
lta_sym = lta_sym %>% fit(data = dat, method = "em")
lta_sym %>% parameter_numbers()
lta_sym %>% estimates()
debug(estimate)
f <- lcpa_sym %>% fit()

f %>% estimates()
f %>% parameter_numbers()
ff <- f %>% fit()



library(glca)
data("gss08")
names(gss08)
ff <- slcm(L1[3] ~ DEFECT + HLTH + RAPE + POOR + SINGLE)
fit <- ff %>% fit(data = gss08)
debug(regress)
fit %>% regress(SEX ~ L1, data = gss08)