Package {TSQLEM}

Title:	Two Stage Estimation for Generalized Structural Equation Models
Version:	0.1.0
Description:	Provides a framework to estimate high dimensional generalized structural equation models using two stage quasi-likelihood expectation-maximization. The structural model supports binomial (logit and probit), Poisson, negative binomial, and gamma distributions for the outcome variable. Hattab (2026) "A Two Stage Quasi-Likelihood Estimation Method for High Dimensional Generalized Structural Equation Models" <doi:10.48550/arXiv.2608.16017>.
License:	GPL (≥ 3)
Encoding:	UTF-8
RoxygenNote:	8.0.0
Depends:	R (≥ 4.6.0)
Imports:	lavaan, Matrix, MASS, mvtnorm
Suggests:	psych
NeedsCompilation:	no
Packaged:	2026-09-17 13:57:55 UTC; Admin
Author:	Mohammad Hattab [aut, cre, cph]
Maintainer:	Mohammad Hattab <mwhattab@unm.edu>
Repository:	CRAN
Date/Publication:	2026-09-28 08:10:08 UTC

Sample the Full Conditional Distribution of the Latent Space

Description

Extracts the full conditional distribution of the latent variables given the observed indicators and the structural outcome via Markov Chain Monte Carlo (MCMC).

Usage

CondLatent(
  df,
  y_name,
  betahat,
  phi = 1,
  family,
  link,
  scores,
  cov_scores,
  nsim = 10000,
  individual_index = NULL,
  method = "MH",
  mh_burnin = 2000,
  mh_tune_c = 1
)

Arguments

Value

A list of individual-level data containing the conditional distribution of the latent space given the observed data, MCMC draws, and the MH acceptance rate.

See Also

TSQLEM

Examples

# 1. PREP THE DATA
# ==========================================
library(psych)
data(bfi)
df <- na.omit(bfi[, c(1:25, 26)]) # Must be complete cases

# 2. DEFINE MEASUREMENT MODEL
# ==========================================
model_cfa <-  '
  Agree =~ A1 + A2 + A3 + A4 + A5
  Cons  =~ C1 + C2 + C3 + C4 + C5
  Extra =~ E1 + E2 + E3 + E4 + E5
  Neuro =~ N1 + N2 + N3 + N4 + N5
  Open  =~ O1 + O2 + O3 + O4 + O5'

# 3. RUN TSQLEM ALGORITHM
# ==========================================
Fit <- TSQLEM(df = df,
              y_name = "gender",
              lav_syn_cfa = model_cfa,
              B = 5000,
              num.it = 20,
              tol = 10^-5,
              cfa_estimator = "ML",
              family = "binomial",
              link = "probit",
              DebiasBern = TRUE,
              probit_exact = TRUE)

print(Fit$beta_hat)

# 4. SAMPLE FULL CONDITIONAL LATENT DISTRIBUTIONS
# ==========================================
# Extract counterfactual/conditional densities for individuals 10, 376, and 1057
mcmc_results <- CondLatent(df = df,
                               y_name = "gender",
                               betahat = Fit$beta_hat,
                               phi = Fit$dispersion_hat, # Inherits 1 for binomial
                               family = "binomial",
                               link = "probit",
                               scores = Fit$scores,
                               cov_scores = Fit$cov_scores,
                               nsim = 10^4,
                               individual_index = c(10, 376, 1057),
                               method = "MH",
                               mh_burnin =3000,
                               mh_tune_c = 1.25)

# Review acceptance rates and the mean of the conditional distributions 
print(mcmc_results[["10"]]$mh_accept_rate)
print(mcmc_results[["10"]]$conditional_means)

Two-Stage Quasi-Likelihood Estimation for Generalized Structural Equation Models

Description

Fits high-dimensional generalized structural equation models (GSEMs) utilizing a Two-Stage Quasi-Likelihood Expectation-Maximization (TSQLEM) framework.

Usage

TSQLEM(
  df,
  y_name,
  lav_syn_cfa,
  B,
  start_beta = NA,
  num.it = 10,
  tol = 10^-4,
  cfa_estimator = "ML",
  family = "binomial",
  link = "logit",
  DebiasBern = TRUE,
  probit_exact = TRUE
)

Arguments

Details

The algorithm operates in two modular stages:

1. Stage 1: Fits a measurement model using lavaan to extract latent factor scores and their conditional covariance matrix.

2. Stage 2: Employs marginal quasi-likelihood estimating equations via an iteratively reweighted least squares (IRLS) scheme, utilizing either exact analytical integration or Monte Carlo simulations.

For Gamma and Negative Binomial models, extra dispersion parameters are estimated internally via profile working models using standard variance parameterizations:

• Gamma: Modeled with Var(y) = dispersion * mean^2.

• Negative Binomial: Modeled with Var(y) = mean + (mean^2 / dispersion).

Supported outcome distributions and link functions include:

• family = "binomial": paired with link = "logit" or link = "probit"

• family = "poisson": paired with link = "log"

• family = "negbin": paired with link = "log"

• family = "gamma": paired with link = "log"

Value

A list containing:

• beta_hat: The final estimated structural parameters.

• dispersion_hat: The estimated dispersion parameter (applicable to Gamma/Negative Binomial).

• BREAK: An integer indicating if the algorithm did not coverage (1) or ran normally (0).

• Num.it: The number of iterations completed.

• eb_beta: The empirical Bayes starting values extracted during initialization.

• beta_history: A matrix tracking parameter estimates across iterations.

• gradient_history: A vector tracking the gradient norm across iterations.

• scores: A matrix containing the Stage 1 empirical Bayes expected latent factor scores for each complete case.

• cov_scores: A matrix representing the approximate conditional covariance of the latent space, extracted from the measurement model.

Author(s)

Mohammad W. Hattab

References

Hattab, M. W. (2026). A Two Stage Quasi-Likelihood Estimation Method for High Dimensional Generalized Structural Equation Models. arXiv:2608.16017. doi:10.48550/arXiv.2608.16017

See Also

CondLatent

Examples

# 1. PREP THE DATA
# ==========================================
library(psych)
data(bfi)
df <- bfi[, c(1:25, 26)] 
df <- na.omit(df) 
model_cfa <- '
   Agree =~ A1 + A2 + A3 + A4 + A5
   Cons  =~ C1 + C2 + C3 + C4 + C5
   Extra =~ E1 + E2 + E3 + E4 + E5
   Neuro =~ N1 + N2 + N3 + N4 + N5
   Open  =~ O1 + O2 + O3 + O4 + O5'

Fit <- TSQLEM(df = df,
               y_name = "gender", 
               lav_syn_cfa = model_cfa,
               B = 5000, 
               num.it = 20, 
               tol = 10^-5, 
               cfa_estimator = "ML", 
               family = "binomial", 
               link = "probit",  
               DebiasBern = TRUE,
               probit_exact = TRUE)

print(Fit$beta_hat)