if (FALSE) {
#############################################################################
# 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)
}
Run the code above in your browser using DataLab