mgsem: Estimation of Multiple-Group Structural Equation Models

View source: R/mgsem.R

mgsemR Documentation

Estimation of Multiple-Group Structural Equation Models

Description

Estimates a multiple-group structural equation model. The function allows arbitrary prior distributions on model parameters and regularized estimation with the SCAD and the LASSO penalty. Moreover, it can also conduct robust moment estimation using the L_p loss function \rho(x)=|x|^p for p \ge 0. See Robitzsch (2023) for more details.

Usage

mgsem(suffstat, model, data=NULL, group=NULL, weights=NULL, estimator="ML",
     p_me=2, p_pen=1, pen_type="scad", diffpar_pen=NULL, pen_sample_size=TRUE,
     a_scad=3.7, eps_approx=0.001, comp_se=TRUE, se_delta_formula=FALSE,
     prior_list=NULL, hessian=TRUE, fixed_parms=FALSE, cd=FALSE,
     cd_control=list(maxiter=20, tol=5*1e-04, interval_length=0.05, method="exact"),
     partable_start=NULL, num_approx=FALSE, technical=NULL, control=list())

Arguments

suffstat

List containing sufficient statistics

model

Model specification, see examples. Can have entries est, index, lower, upper, prior, pen_l2, pen_lp, pen_difflp. Each entry can be defined for model matrices ALPHA, NU, LAM, PHI, and PSI.

data

Optional data frame

group

Optional vector of group identifiers

weights

Optional vector of sampling weights

estimator

Character. Can be either "ML" for maximum likelihood fitting function or "ME" for robust moment estimation.

p_me

Power in $L_p$ loss function for robust moment estimation

p_pen

Power for penalty in regularized estimation. For regular LASSO and SCAD penalty functions, it is $p=1$.

pen_type

Penalty type. Can be either "scad" or "lasso".

diffpar_pen

List containing values of regularization parameters in fused lasso estimation

pen_sample_size

List containing values for sample sizes for regularized estimation

a_scad

Parameter $a$ used in SCAD penalty

eps_approx

Approximation value for nondifferentiable robust moment fitting function or penalty function

comp_se

Logical indicating whether standard errors should be computed

se_delta_formula

Logical indicating whether standard errors should be computed according to the delta formula

prior_list

List containing specifications of the prior distributions

hessian

Logical indicating whether the Hessian matrix should be computed

fixed_parms

Logical indicating whether all model parameters should be fixed

cd

Logical indicating whether coordinate descent should be used for estimation

cd_control

Control parameters for coordinate descent estimation

partable_start

Starting values for parameter estimation

num_approx

Logical indicating whether derivatives should be computed based on numerical differentiation

technical

Parameters used for optimization in sirt_optimizer

control

Control paramaters for optimization

Details

[MORE INFORMATION TO BE ADDED]

Value

A list with following values

coef

Coeffients

vcov

Variance matrix

se

Vector of standard errors

partable

Parameter table

model

Specified model

opt_res

Result from optimization

opt_value

Value of fitting function

residuals

Residuals of sufficient statistics

ic

Information criteria

technical

Specifications of optimizer

suffstat_vcov

Variance matrix of sufficient statistics

me_delta_method

Input and output matrices for delta method if estimator="ME"

data_proc

Processed data

case_ll

Casewise log-likelihood function

...

Further values

References

Robitzsch, A. (2023). Model-robust estimation of multiple-group structural equation models. Algorithms, 16(4), 210. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3390/a16040210")}

Examples

## Not run: 
#############################################################################
# EXAMPLE 1: Noninvariant item intercepts in a multiple-group SEM
#############################################################################

#---- simulate data
set.seed(65)
G <- 3  # number of groups
I <- 5  # number of items
# define lambda and nu parameters
lambda <- matrix(1, nrow=G, ncol=I)
nu <- matrix(0, nrow=G, ncol=I)
err_var <- matrix(1, nrow=G, ncol=I)

# define extent of noninvariance
dif_int <- .5

#- 1st group: N(0,1)
nu[1,4] <- dif_int
#- 2nd group: N(0.3,1.5)
gg <- 2 ;
nu[gg,1] <- -dif_int
#- 3nd group: N(.8,1.2)
gg <- 3
nu[gg,2] <- -dif_int
#- define distributions of groups
mu <- c(0,.3,.8)
sigma <- sqrt(c(1,1.5,1.2))
N <- rep(1000,3) # sample sizes per group

exact <- FALSE
suffstat <- sirt::invariance_alignment_simulate(nu, lambda, err_var, mu, sigma, N,
                output="suffstat", groupwise=TRUE, exact=exact)

#---- model specification

# model specifications joint group
est <- list(
        ALPHA=matrix( c(0), ncol=1),
        NU=matrix( 0, nrow=I, ncol=1),
        LAM=matrix(1, nrow=I, ncol=1),
        PHI=matrix(0,nrow=1,ncol=1),
        PSI=diag(rep(1,I))
        )

# parameter index
index <- list(
        ALPHA=0*est$ALPHA,
        NU=1+0*est$NU,
        LAM=1+0*est$LAM,
        PHI=0*est$PHI,
        PSI=diag(1,I)
        )

# lower bounds
lower <- list(
        PSI=diag(rep(0.01,I)), PHI=matrix(0.01,1,1)
        )

#*** joint parameters
group0 <- list(est=est, index=index, lower=lower)

#*** group1
est <- list(
        ALPHA=matrix( c(0), ncol=1),
        NU=matrix( 0, nrow=I, ncol=1),
        LAM=matrix(0, nrow=I, ncol=1),
        PHI=matrix(1,nrow=1,ncol=1)
            )

# parameter index
index <- list(
        ALPHA=0*est$ALPHA,
        NU=0*est$NU,
        LAM=1*est$LAM,
        PHI=0*est$PHI
        )

group1 <- list(est=est, index=index, lower=lower)

#*** group 2 and group 3

# modify parameter index
index$ALPHA <- 1+0*est$ALPHA
index$PHI <- 1+0*est$PHI
group3 <- group2 <- list(est=est, index=index, lower=lower)

#*** define model
model <- list(group0=group0, group1=group1, group2=group2, group3=group3)

#-- estimate model with ML
res1 <- sirt::mgsem( suffstat=suffstat, model=model2, eps_approx=1e-4, estimator="ML",
                    technical=list(maxiter=500, optimizer="optim"),
                    hessian=FALSE, comp_se=FALSE, control=list(trace=1) )
str(res1)

#-- robust moment estimation with p=0.5

optimizer <- "optim"
technical <- list(maxiter=500, optimizer=optimizer)
eps_approx <- 1e-3

res2 <- sirt::mgsem( suffstat=suffstat, model=res1$model, p_me=0.5,
                    eps_approx=eps_approx, estimator="ME", technical=technical,
                    hessian=FALSE, comp_se=FALSE, control=list(trace=1) )

#---- regularized estimation

nu_lam <- 0.1    # regularization parameter

# redefine model
define_model <- function(model, nu_lam)
{
    pen_lp <- list( NU=nu_lam+0*model$group1$est$NU)
    ee <- "group1"
    for (ee in c("group1","group2","group3"))
    {
        model[[ee]]$index$NU <- 1+0*index$NU
        model[[ee]]$pen_lp <- pen_lp
    }
    return(model)
}

model3 <- define_model(model=model, nu_lam=nu_lam)
pen_type <- "scad"

res3 <- sirt::mgsem( suffstat=suffstat, model=model3, p_pen=1, pen_type=pen_type,
                    eps_approx=eps_approx, estimator="ML",
                    technical=list(maxiter=500, optimizer="optim"),
                    hessian=FALSE, comp_se=FALSE, control=list(trace=1) )
str(res3)

## End(Not run)

sirt documentation built on May 29, 2024, 8:43 a.m.