R/estModel.R

In slca: Structural Modeling for Multiple Latent Class Variables

estModel <- function(method, control, par, mf, arg) {
   llf <- function(
      logit, fix0, fix1, ref, id, y, nobs, nvar, nlev, nlv, nrl, nlf,
      npi, ntau, nrho, ul, vl, lf, tr, rt, eqrl, eqlf,
      nc, nk, nl, ncl, nc_pi, nk_tau, nl_tau, nc_rho, nr_rho
   ) {
      logits <- numeric(length(id))
      logits[-c(fix0, fix1, ref)] <- logit
      logits[fix0] <- -Inf
      logits[fix1] <- Inf
      logits[ref] <- 0
      par <- unlist(tapply(logits, id, norm2))

      fll(y, par, nobs, nvar, nlev, nlv, nrl, nlf,
          npi, ntau, nrho, ul, vl, lf, tr, rt, eqrl, eqlf,
          nc, nk, nl, ncl, nc_pi, nk_tau, nl_tau, nc_rho, nr_rho)
   }
   while (any(cond <- arg$ref_idx %in% which(arg$fix0))) {
      arg$ref[cond] <- arg$ref[cond] - 1
      arg$ref_idx[cond] <- arg$ref_idx[cond] - 1
   }

   if (method == "em") {
      if (control$verbose) cat("EM iteration begin.\n")
      em <- em_est(
         attr(mf, "y"), arg$nobs, arg$nvar, unlist(arg$nlev), par, arg$fix0,
         arg$nlv, arg$nrl, arg$nlf, arg$npi, arg$ntau, arg$nrho,
         arg$ul, arg$vl, arg$lf, arg$tr, arg$rt, arg$eqrl, arg$eqlf,
         arg$nc, arg$nk, arg$nl, arg$ncl,
         arg$nc_pi, arg$nk_tau, arg$nl_tau, arg$nc_rho, arg$nr_rho,
         control$em.iterlim, control$em.tol, control$verbose
      )
      par <- em$param
      logit <- par - par[arg$ref_idx[arg$id]]
      fix0 <- which(logit == -Inf)
      fix1 <- which(logit == Inf)
      em.conv <- em$converged
      em.niter <- em$niter
      nlm.conv <- NA
      if (control$verbose) cat(".. done.\n")
   } else if (method == "nlm") {
      if (control$verbose) cat("nlm iteration begin.\n")
      logit <- par - par[arg$ref_idx[arg$id]]
      fix0 <- which(logit == -Inf)
      fix1 <- which(logit == Inf)

      nonlm <- stats::nlm(
         llf, logit[-c(fix0, fix1, arg$ref_idx)],
         fix0, fix1, arg$ref_idx, arg$id, y = attr(mf, "y"),
         nobs = arg$nobs, nvar = arg$nvar, nlev = unlist(arg$nlev),
         nlv = arg$nlv, nrl = arg$nrl, nlf = arg$nlf,
         npi = arg$npi, ntau = arg$ntau, nrho = arg$nrho,
         ul = arg$ul, vl = arg$vl, lf = arg$lf, tr = arg$tr, rt = arg$rt,
         eqrl = arg$eqrl, eqlf = arg$eqlf,
         nc = arg$nc, nk = arg$nk, nl = arg$nl, ncl = arg$ncl,
         nc_pi = arg$nc_pi, nk_tau = arg$nk_tau, nl_tau = arg$nl_tau,
         nc_rho = arg$nc_rho, nr_rho = arg$nr_rho,
         iterlim = control$nlm.iterlim,
         gradtol = control$nlm.tol, steptol = control$nlm.tol
      )
      logit[-c(fix0, fix1, arg$ref_idx)] <- nonlm$estimate
      par <- unlist(tapply(logit, arg$id, norm2))
      em.conv <- NA
      nlm.conv <- nonlm$code < 3
      if (control$verbose) cat(".. done.\n")
   } else if (method == "hybrid") {
      if (control$verbose) cat("EM iteration begin.\n")
      em <- em_est(
         attr(mf, "y"), arg$nobs, arg$nvar, unlist(arg$nlev), par, arg$fix0,
         arg$nlv, arg$nrl, arg$nlf, arg$npi, arg$ntau, arg$nrho,
         arg$ul, arg$vl, arg$lf, arg$tr, arg$rt, arg$eqrl, arg$eqlf,
         arg$nc, arg$nk, arg$nl, arg$ncl,
         arg$nc_pi, arg$nk_tau, arg$nl_tau, arg$nc_rho, arg$nr_rho,
         control$em.iterlim, control$em.tol, control$verbose
      )
      par <- em$param
      if (control$verbose) cat(".. done. \nnlm iteration begin.\n")
      logit <- par - par[arg$ref_idx[arg$id]]
      fix0 <- which(logit == -Inf)
      fix1 <- which(logit == Inf)

      nonlm <- stats::nlm(
         llf, logit[-c(fix0, fix1, arg$ref_idx)],
         fix0 = fix0, fix1 = fix1, ref = arg$ref_idx,
         id = arg$id, y = attr(mf, "y"),
         nobs = arg$nobs, nvar = arg$nvar, nlev = unlist(arg$nlev),
         nlv = arg$nlv, nrl = arg$nrl, nlf = arg$nlf,
         npi = arg$npi, ntau = arg$ntau, nrho = arg$nrho,
         ul = arg$ul, vl = arg$vl, lf = arg$lf, tr = arg$tr, rt = arg$rt,
         eqrl = arg$eqrl, eqlf = arg$eqlf,
         nc = arg$nc, nk = arg$nk, nl = arg$nl, ncl = arg$ncl,
         nc_pi = arg$nc_pi, nk_tau = arg$nk_tau, nl_tau = arg$nl_tau,
         nc_rho = arg$nc_rho, nr_rho = arg$nr_rho,
         iterlim = control$nlm.iterlim,
         gradtol = control$nlm.tol, steptol = control$nlm.tol
      )
      logit[-c(fix0, fix1, arg$ref_idx)] <- nonlm$estimate
      par <- unlist(tapply(logit, arg$id, norm2))

      em.conv  <- em$converged
      nlm.conv <- nonlm$code < 3
      if (control$verbose) cat(".. done.\n")
   }

   hess <- matrix(NA, length(par), length(par))
   if (control$hessian) {
      nonlm <- stats::nlm(
         llf, logit[-c(fix0, fix1, arg$ref_idx)],
         fix0, fix1, arg$ref_idx, arg$id, y = attr(mf, "y"),
         nobs = arg$nobs, nvar = arg$nvar, nlev = unlist(arg$nlev),
         nlv = arg$nlv, nrl = arg$nrl, nlf = arg$nlf,
         npi = arg$npi, ntau = arg$ntau, nrho = arg$nrho,
         ul = arg$ul, vl = arg$vl, lf = arg$lf, tr = arg$tr, rt = arg$rt,
         eqrl = arg$eqrl, eqlf = arg$eqlf,
         nc = arg$nc, nk = arg$nk, nl = arg$nl, ncl = arg$ncl,
         nc_pi = arg$nc_pi, nk_tau = arg$nk_tau, nl_tau = arg$nl_tau,
         nc_rho = arg$nc_rho, nr_rho = arg$nr_rho,
         iterlim = 1, hessian = TRUE
      )
      ind <- c(fix0, fix1, arg$ref_idx)
      hess[-ind, -ind] <- nonlm$hessian
   }

   list(par = par, logit = logit, hess = hess,
        conv = c(EM = em.conv, nlm = nlm.conv))
}