## File Name: mgsem.R
## File Version: 0.554
mgsem <- function(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=1e-3, 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-4, interval_length=0.05, method="exact"),
partable_start=NULL, num_approx=FALSE, technical=NULL, control=list() )
{
#- pen_type: lasso, scad or none
#*** preliminaries
CALL <- match.call()
s1 <- Sys.time()
test <- FALSE
# test <- TRUE
#** coordinate decent
if (cd){
fixed_parms <- TRUE
}
#*** process data if suffstat not available
if (missing(suffstat)){
data_proc <- mgsem_proc_data(data=data, group=group, weights=weights)
groups <- data_proc$groups
suffstat <- data_proc$suffstat
weights <- data_proc$weights
is_data <- TRUE
G <- length(groups)
} else {
groups <- names(suffstat)
G <- length(groups)
if (is.null(groups)){
groups <- paste0('Group',1L:G)
}
data_proc <- NULL
is_data <- FALSE
}
technical$is_data <- is_data
#*** compute covariance matrix of sufficient statistics
suffstat_vcov <- mgsem_suffstat_covariance_matrix(suffstat=suffstat)
#*** process technical defaults
res <- mgsem_proc_technical(technical=technical, control=control, p_me=p_me,
p_pen=p_pen, eps_approx=eps_approx, suffstat=suffstat,
estimator=estimator, diffpar_pen=diffpar_pen, cd_control=cd_control)
technical <- res$technical
technical$estimator <- estimator
control <- res$control
eps_approx <- res$eps_approx
p <- res$p
cd_control <- res$cd_control
#*** process sufficient statistics
res <- mgsem_proc_suffstat(suffstat=suffstat)
suffstat <- res$suffstat
G <- res$G
I <- res$I
N <- res$N
N_group <- res$N_group
W <- res$W
random_sd <- -9
if (test){
random_sd <- 1e-1
}
#*** process model specification
res <- mgsem_proc_model(model=model, G=G, prior_list=prior_list,
technical=technical, N_group=N_group, random_sd=random_sd,
pen_type=pen_type, fixed_parms=fixed_parms,
partable_start=partable_start, diffpar_pen=diffpar_pen,
pen_sample_size=pen_sample_size, W=W)
model <- res$model
partable <- res$partable
NP <- res$NP
coef <- res$coef
ND <- res$ND
D <- res$D
is_B <- res$is_B
technical <- res$technical
types <- res$types
difflp_info <- res$difflp_info
loop_parms <- res$loop_parms
if (estimator=='ML'){
eval_fun <- 'mgsem_loglike_suffstat'
grad_param_fun <- 'mgsem_loglike_suffstat_derivative_parameter'
grad_suffstat_fun <- 'mgsem_loglike_suffstat_derivative'
}
if (estimator=='ME'){
eval_fun <- 'mgsem_loss_function_suffstat'
grad_param_fun <- 'mgsem_loss_function_suffstat_derivative_parameter'
grad_suffstat_fun <- 'mgsem_loss_function_suffstat'
}
if (technical$use_deriv){
grad_fun <- mgsem_grad_fun
} else {
grad_fun <- mgsem_grad_fun_numeric_approx
}
x <- coef
opt_fun_args <- list(NP=NP, ND=ND, suffstat=suffstat, model=model,
partable=partable, G=G, I=I, D=D, is_B=is_B,
N=N, N_group=N_group,
estimator=estimator, eval_fun=eval_fun,
grad_param_fun=grad_param_fun,
grad_suffstat_fun=grad_suffstat_fun,
technical=technical, p_pen=p_pen, p_me=p_me,
p=p, eps_approx=eps_approx, num_approx=num_approx,
prior_list=prior_list, use_penalty=TRUE, types=types,
difflp_info=difflp_info, loop_parms=loop_parms,
pen_type=pen_type, a_scad=a_scad )
#- test function
res <- mgsem_test_fun(test=test, coef=coef, opt_fun_args=opt_fun_args)
#**** estimation
if (! fixed_parms){
#- define lower and upper bounds
ind <- which(partable$unique>0)
lower <- partable[ ind, 'lower' ]
upper <- partable[ ind, 'upper' ]
#- use optimizer
opt_res <- sirt_optimizer(optimizer=technical$optimizer, par=x, fn=mgsem_opt_fun,
grad=grad_fun, opt_fun_args=opt_fun_args,
method='L-BFGS-B', lower=lower, upper=upper,
hessian=hessian, control=control )
} else {
#**** no estimation
ll <- mgsem_opt_fun(x=x, opt_fun_args=opt_fun_args)
opt_res <- list(par=x, opt_value=ll)
#--- coordinate descent
if (cd){
opt_res <- mgsem_cd_opt(x=x, opt_fun_args=opt_fun_args, tol=cd_control$tol,
maxiter=cd_control$maxiter, method=cd_control$method,
interval_length=cd_control$interval_length)
}
}
#-- collect results
coef <- opt_res$par
opt_value <- opt_res$value
partable <- mgsem_coef2partable(coef=coef, partable=partable)
model <- mgsem_partable2model(partable=partable, model=model)
#-- evaluate optimization function and gradient
opt_fun_args$partable <- partable
opt_fun_args$model <- model
opt_fun_output <- mgsem_opt_fun(x=coef, opt_fun_args=opt_fun_args, output_all=TRUE)
implied <- opt_fun_output$implied
est_tot <- opt_fun_output$est_tot
grad_fun_output <- mgsem_grad_fun(x=coef, opt_fun_args=opt_fun_args, output_all=TRUE)
#-- vcov for estimator='ME'
res <- mgsem_vcov_me(coef=coef, opt_fun_args=opt_fun_args,
suffstat_vcov=suffstat_vcov, comp_se=comp_se,
se_delta_formula=se_delta_formula)
vcov <- res$vcov
se <- res$se
comp_se_me <- res$comp_se_me
me_delta_method <- res
#-- residual statistics
residuals <- mgsem_output_proc_residuals(implied=implied, suffstat=suffstat)
#-- compute case-wise log-likelihood
case_ll <- mgsem_output_proc_casewise_likelihood(data_proc=data_proc,
implied=implied, estimator=estimator)
#-- computation of standard errors
res <- mgsem_observed_information(coef=coef, opt_fun_args=opt_fun_args,
technical=technical, comp_se=comp_se, comp_se_me=comp_se_me)
if (!comp_se_me){
vcov <- res$vcov
comp_se <- res$comp_se
se <- res$se
}
info_loglike <- res$info_loglike
info_loglike_pen <- res$info_loglike_pen
#-- information criteria
ic <- mgsem_ic(opt_fun_output=opt_fun_output, opt_fun_args=opt_fun_args,
partable=partable, technical=technical)
#--- output
s2 <- Sys.time()
res <- list(coef=coef, vcov=vcov, se=se,
partable=partable, model=model, suffstat=suffstat,
opt_res=opt_res, opt_value=opt_value, implied=implied,
est_tot=est_tot, residuals=residuals, opt_fun_output=opt_fun_output,
ic=ic, info_loglike=info_loglike, info_loglike_pen=info_loglike_pen,
estimator=estimator, p_pen=p_pen, p_me=p_me,
pen_type=pen_type, eps_approx=eps_approx,
technical=technical, comp_se=comp_se, groups=groups, group=group,
data=data, data_proc=data_proc, case_ll=case_ll,
suffstat_vcov=suffstat_vcov, me_delta_method=me_delta_method,
CALL=CALL, s1=s1, s2=s2)
class(res) <- 'mgsem'
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.