diff --git a/DESCRIPTION b/DESCRIPTION
index bd8bd43..8b86cae 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.13-105
-Date: 2023-03-19 12:31:21
+Version: 3.13-128
+Date: 2023-04-02 12:30:22
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <robitzsch@ipn.uni-kiel.de>
 Description: 
diff --git a/R/L0_polish.R b/R/L0_polish.R
index 4b1bf37..58171b4 100644
--- a/R/L0_polish.R
+++ b/R/L0_polish.R
@@ -1,5 +1,5 @@
 ## File Name: L0_polish.R
-## File Version: 0.12
+## File Version: 0.131
 
 
 L0_polish <- function(x, tol, conv=0.01, maxiter=30, type=1, verbose=TRUE)
@@ -9,9 +9,9 @@ L0_polish <- function(x, tol, conv=0.01, maxiter=30, type=1, verbose=TRUE)
     while(res$iterate_further){
         res <- L0_polish_one_iteration(x=res$x_update, tol=tol, type=type, eps=conv)
         if (verbose){
-            v1 <- paste0("Interactions detected: ", res$N_elim)
-            v2 <- paste0(" | Absolute value residual: ", round(res$max_resid,3) )
-            cat(v1, v2, "\n")
+            v1 <- paste0('Interactions detected: ', res$N_elim)
+            v2 <- paste0(' | Absolute value residual: ', round(res$max_resid,3) )
+            cat(v1, v2, '\n')
             utils::flush.console()
         }
     }
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 8a03ac4..691cb05 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.013105
+## File Version: 3.013128
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
diff --git a/R/dexppow.R b/R/dexppow.R
index c2500c4..2ea66f1 100644
--- a/R/dexppow.R
+++ b/R/dexppow.R
@@ -1,5 +1,5 @@
 ## File Name: dexppow.R
-## File Version: 0.04
+## File Version: 0.052
 
 
 ## copied from normalp::dnormp
@@ -10,6 +10,8 @@ dexppow <- function (x, mu=0, sigmap=1, pow=2, log=FALSE)
     expon1 <- (abs(x - mu))^p
     expon2 <- p * sigmap^p
     dsty <- (1/cost) * exp(-expon1/expon2)
-    if (log){ dsty <- log(dsty) }
+    if (log){
+        dsty <- log(dsty)
+    }
     return(dsty)
 }
diff --git a/R/lq_fit.R b/R/lq_fit.R
index 883cbdd..58a981f 100644
--- a/R/lq_fit.R
+++ b/R/lq_fit.R
@@ -1,5 +1,5 @@
 ## File Name: lq_fit.R
-## File Version: 0.152
+## File Version: 0.154
 
 lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL,
         est_pow=FALSE, optimizer="optim", eps_vec=10^seq(0,-10, by=-.5),
@@ -71,7 +71,7 @@ lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL,
                 iterate_powers <- FALSE
             }
             iter <- iter + 1
-            # cat( paste0("eps=",eps, " | pow=",pow,"\n"))
+            # cat( paste0('eps=',eps, ' | pow=',pow,'\n'))
         }
     }
 
diff --git a/R/lq_fit_estimate_power.R b/R/lq_fit_estimate_power.R
index 4aa1500..95141b4 100644
--- a/R/lq_fit_estimate_power.R
+++ b/R/lq_fit_estimate_power.R
@@ -1,5 +1,5 @@
 ## File Name: lq_fit_estimate_power.R
-## File Version: 0.14
+## File Version: 0.15
 
 
 lq_fit_estimate_power <- function(e, pow_init=2, lower_pow=.1, upper_pow=10)
@@ -20,7 +20,7 @@ lq_fit_estimate_power <- function(e, pow_init=2, lower_pow=.1, upper_pow=10)
             (vi - sqrt(gamma(1/p) * gamma(3/p))/gamma(2/p))^2
         }
         res <- stats::optim(par=pow0, fn=fvi, lower=lower_pow, upper=upper_pow,
-                            method="L-BFGS-B")
+                            method='L-BFGS-B')
         p <- res$par
     }
     return(p)
diff --git a/R/lsdm_est_logist_quant.R b/R/lsdm_est_logist_quant.R
index e138568..b769fe1 100644
--- a/R/lsdm_est_logist_quant.R
+++ b/R/lsdm_est_logist_quant.R
@@ -1,5 +1,5 @@
 ## File Name: lsdm_est_logist_quant.R
-## File Version: 0.192
+## File Version: 0.193
 
 
 #--- Function for calculating logistic functions and probability quantiles
@@ -11,8 +11,8 @@ lsdm_est_logist_quant <- function( probcurves, theta, quantiles, wgt_theta,
     b0 <- NULL
     if (est.icc){
         pars.probcurves <- matrix( 0, nrow=I, ncol=5 )
-        colnames(pars.probcurves) <- c("b.2PL", "a.2PL", "sigma.2PL", "b.1PL",
-                                            "sigma.1PL")
+        colnames(pars.probcurves) <- c('b.2PL', 'a.2PL', 'sigma.2PL', 'b.1PL',
+                                            'sigma.1PL')
         rownames(pars.probcurves) <- rownames(probcurves)
         for (kk in 1:I){
             if (!is.null(b)){
@@ -33,7 +33,7 @@ lsdm_est_logist_quant <- function( probcurves, theta, quantiles, wgt_theta,
                 } )
             } )
     probcurves.quant <- as.data.frame(probcurves.quant)
-    colnames(probcurves.quant) <- paste( "Q", 100*quantiles, sep="")
+    colnames(probcurves.quant) <- paste( 'Q', 100*quantiles, sep='')
     rownames(probcurves.quant) <- rownames(probcurves)
     if (est.icc){
         pars.probcurves <- cbind( probcurves.quant, pars.probcurves )
diff --git a/R/mgsem.R b/R/mgsem.R
index e48e4f9..673cd0a 100644
--- a/R/mgsem.R
+++ b/R/mgsem.R
@@ -1,9 +1,11 @@
 ## File Name: mgsem.R
-## File Version: 0.423
+## File Version: 0.497
 
 mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
         estimator="ML", p_me=2, p_pen=1, pen_type="scad",
-        a_scad=3.7, eps_approx=1e-3, comp_se=TRUE, prior_list=NULL, hessian=TRUE,
+        diffpar_pen=NULL, a_scad=3.7, eps_approx=1e-3, comp_se=TRUE,
+        se_delta_formula=FALSE,
+        prior_list=NULL, hessian=TRUE,
         fixed_parms=FALSE, partable_start=NULL,
         num_approx=FALSE, technical=NULL, control=list() )
 {
@@ -23,17 +25,25 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
         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',1: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)
+                estimator=estimator, diffpar_pen=diffpar_pen)
     technical <- res$technical
     technical$estimator <- estimator
     control <- res$control
@@ -56,7 +66,7 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
     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)
+                    partable_start=partable_start, diffpar_pen=diffpar_pen)
     model <- res$model
     partable <- res$partable
     NP <- res$NP
@@ -69,12 +79,12 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
     difflp_info <- res$difflp_info
     loop_parms <- res$loop_parms
 
-    if (estimator=="ML"){
+    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"){
+    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'
@@ -107,13 +117,13 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
 
         #- define lower and upper bounds
         ind <- which(partable$unique>0)
-        lower <- partable[ ind, "lower" ]
-        upper <- partable[ ind, "upper" ]
+        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,
+                        method='L-BFGS-B', lower=lower, upper=upper,
                         hessian=hessian, control=control )
     } else {
         #**** no estimation
@@ -135,6 +145,14 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
     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
+
     #-- residual statistics
     residuals <- mgsem_output_proc_residuals(implied=implied, suffstat=suffstat)
 
@@ -144,10 +162,12 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
 
     #-- computation of standard errors
     res <- mgsem_observed_information(coef=coef, opt_fun_args=opt_fun_args,
-                    technical=technical, comp_se=comp_se)
-    vcov <- res$vcov
-    comp_se <- res$comp_se
-    se <- res$se
+                    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
 
@@ -166,7 +186,7 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
                     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,
-                    CALL=CALL, s1=s1, s2=s2)
-    class(res) <- "mgsem"
+                    suffstat_vcov=suffstat_vcov, CALL=CALL, s1=s1, s2=s2)
+    class(res) <- 'mgsem'
     return(res)
 }
diff --git a/R/mgsem_bdiag.R b/R/mgsem_bdiag.R
new file mode 100644
index 0000000..a5621a8
--- /dev/null
+++ b/R/mgsem_bdiag.R
@@ -0,0 +1,14 @@
+## File Name: mgsem_bdiag.R
+## File Version: 0.05
+
+mgsem_bdiag <- function(x1, x2)
+{
+    vars <- c(rownames(x1),rownames(x2))
+    n1 <- ncol(x1)
+    n2 <- ncol(x2)
+    mat <- matrix(0,nrow=n1+n2,ncol=n1+n2)
+    rownames(mat) <- colnames(mat)
+    mat[1:n1,1:n1] <- x1
+    mat[n1+1:n2,n1+1:n2] <- x2
+    return(mat)
+}
diff --git a/R/mgsem_compute_model_implied_moments.R b/R/mgsem_compute_model_implied_moments.R
index 5711cd8..a94beb1 100644
--- a/R/mgsem_compute_model_implied_moments.R
+++ b/R/mgsem_compute_model_implied_moments.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_compute_model_implied_moments.R
-## File Version: 0.168
+## File Version: 0.169
 
 
 mgsem_compute_model_implied_moments <- function(est, is_B=FALSE, calc_Sigma=TRUE,
@@ -8,11 +8,10 @@ mgsem_compute_model_implied_moments <- function(est, is_B=FALSE, calc_Sigma=TRUE
     Mu <- NULL
     Sigma <- NULL
     if (is_B){
-        requireNamespace("MASS")
         D <- ncol(est$PHI)
         ID <- diag(rep(1,D))
         B <- res$B
-        B1 <- MASS::ginv(ID-B)
+        B1 <- mgsem_ginv(ID-B)
         LAMB <- est$LAM %*% B1
         Mu <- LAMB %*% est$ALPHA + est$NU
         # Sigma <- LAMB %*% tcrossprod(est$PHI, LAMB) + est$PSI
diff --git a/R/mgsem_duplication_matrix.R b/R/mgsem_duplication_matrix.R
new file mode 100644
index 0000000..4180bbf
--- /dev/null
+++ b/R/mgsem_duplication_matrix.R
@@ -0,0 +1,26 @@
+## File Name: mgsem_duplication_matrix.R
+## File Version: 0.11
+
+mgsem_duplication_matrix <- function(x)
+{
+    NV <- ncol(x)
+    dfr1 <- data.frame( index1=rep(1:NV, NV), index2=rep(1:NV, each=NV))
+    N1 <- nrow(dfr1)
+    dfr2 <- dfr1[ dfr1[,1] >=dfr1[,2], ]
+    ND <- nrow(dfr2)
+    dupl <- matrix(0, nrow=ND, ncol=N1)
+    dd <- 1
+    for (dd in 1:ND){
+        h1 <- dfr2$index1[dd]
+        h2 <- dfr2$index2[dd]
+        i1 <- which( ( dfr1$index1==h1 ) & ( dfr1$index2==h2 ) )
+        i2 <- which( ( dfr1$index2==h1 ) & ( dfr1$index1==h2 ) )
+        if (i1==i2){
+            dupl[dd,i1] <- 1
+        } else {
+            dupl[dd,i1] <- 0.5
+            dupl[dd,i2] <- 0.5
+        }
+    }
+    return(dupl)
+}
diff --git a/R/mgsem_eval_lp_penalty_matrix.R b/R/mgsem_eval_lp_penalty_matrix.R
index e76d0e6..3ca424c 100644
--- a/R/mgsem_eval_lp_penalty_matrix.R
+++ b/R/mgsem_eval_lp_penalty_matrix.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_eval_lp_penalty_matrix.R
-## File Version: 0.083
+## File Version: 0.084
 
 mgsem_eval_lp_penalty_matrix <- function(x, fac, p, n, h, eps_approx,
     pen_type="lasso", a_scad=3.7)
@@ -8,11 +8,11 @@ mgsem_eval_lp_penalty_matrix <- function(x, fac, p, n, h, eps_approx,
     I1 <- length(x1)
     y <- matrix(x1, nrow=I1, ncol=I1)-sirt_matrix2(x=x1, nrow=I1)
     y <- mgsem_power_fun_differentiable_approx(x=y, p=p,
-                            eps=eps_approx, deriv=FALSE, approx_method="lp")
-    if (pen_type=="lasso"){
+                            eps=eps_approx, deriv=FALSE, approx_method='lp')
+    if (pen_type=='lasso'){
         val <- fac*y
     }
-    if (pen_type=="scad"){
+    if (pen_type=='scad'){
         val <- mgsem_scad_penalty(x=y, lambda=fac, a=a_scad)
     }
     res <- sum(n*val)
diff --git a/R/mgsem_eval_lp_penalty_vector.R b/R/mgsem_eval_lp_penalty_vector.R
index 7a22bdc..354e553 100644
--- a/R/mgsem_eval_lp_penalty_vector.R
+++ b/R/mgsem_eval_lp_penalty_vector.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_eval_lp_penalty_vector.R
-## File Version: 0.198
+## File Version: 0.199
 
 
 mgsem_eval_lp_penalty_vector <- function(x, fac, n, p, eps_approx, deriv, h, a=3.7,
@@ -7,25 +7,25 @@ mgsem_eval_lp_penalty_vector <- function(x, fac, n, p, eps_approx, deriv, h, a=3
 {
 
     # Lasso penalty
-    if (pen_type=="lasso"){
+    if (pen_type=='lasso'){
         val <- mgsem_power_fun_differentiable_approx(x=x, p=p,
-                            eps=eps_approx, deriv=deriv, approx_method="lp")
+                            eps=eps_approx, deriv=deriv, approx_method='lp')
         val <- n*fac*val
     }
 
     # SCAD penalty
-    if (pen_type=="scad"){
+    if (pen_type=='scad'){
         if (deriv){
             val <- mgsem_power_fun_differentiable_approx(x=x+h, p=p,
-                                eps=eps_approx, deriv=FALSE, approx_method="lp")
+                                eps=eps_approx, deriv=FALSE, approx_method='lp')
             val1 <- mgsem_scad_penalty(x=val, lambda=fac, a=a)
             val <- mgsem_power_fun_differentiable_approx(x=x-h, p=p,
-                                eps=eps_approx, deriv=FALSE, approx_method="lp")
+                                eps=eps_approx, deriv=FALSE, approx_method='lp')
             val2 <- mgsem_scad_penalty(x=val, lambda=fac, a=a)
             val <- n*(val1-val2)/(2*h)
         } else {
             val <- mgsem_power_fun_differentiable_approx(x=x, p=p,
-                                eps=eps_approx, deriv=FALSE, approx_method="lp")
+                                eps=eps_approx, deriv=FALSE, approx_method='lp')
             val <- n*mgsem_scad_penalty(x=val, lambda=fac, a=a)
         }
     }
diff --git a/R/mgsem_evaluate_penalties.R b/R/mgsem_evaluate_penalties.R
index 1e39604..da8e5fd 100644
--- a/R/mgsem_evaluate_penalties.R
+++ b/R/mgsem_evaluate_penalties.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_evaluate_penalties.R
-## File Version: 0.294
+## File Version: 0.329
 
 
 mgsem_evaluate_penalties <- function(x, partable, prior_list, technical,
@@ -22,7 +22,7 @@ mgsem_evaluate_penalties <- function(x, partable, prior_list, technical,
 
     #- loop over parameters
     for (dd in loop_parms){
-        index <- partable[dd,"index"]
+        index <- partable[dd,'index']
         res <- mgsem_evaluate_penalties_evaluate_entry(
                     x=x, res=res, dd=dd, index=index,
                     partable=partable, technical=technical,
@@ -108,7 +108,45 @@ mgsem_evaluate_penalties <- function(x, partable, prior_list, technical,
         res$pen_difflp <- z
     }
 
-    res$pen_all <- res$pen_prior - res$pen_l2 - res$pen_lp - res$pen_difflp
+    #*** penalty for diffpar
+    res$pen_diffpar_lp <- 0
+    if (technical$is_diffpar_pen){
+        diffpar_pen_list_entries <- technical$diffpar_pen$diffpar_pen_list_entries
+        NDP <- nrow(diffpar_pen_list_entries)
+        p <- technical$diffpar_pen$p
+        n <- partable2$N_group
+        # vector of differences of parameter
+        z <- x[ diffpar_pen_list_entries$index1 ] - x[ diffpar_pen_list_entries$index2 ]
+        n2 <- sqrt(n[ diffpar_pen_list_entries$index1 ] *
+                            n[ diffpar_pen_list_entries$index2 ] )
+        args_pen <- list(x=z, fac=diffpar_pen_list_entries$W, n=n2, p=p,
+                        eps_approx=eps_approx, pen_type=pen_type, h=h, deriv=deriv)
+        fun_pen <- 'mgsem_eval_lp_penalty_vector'
+        val <- do.call(what=fun_pen, args=args_pen)
+
+        #* no derivative
+        if (!deriv){
+            val <- sum(val)
+        }
+
+        #* derivative
+        if (deriv){
+            der_z <- val
+            NP <- length(x)
+            val <- rep(0,NP)
+            for (hh in 1:NDP){
+                i1 <- diffpar_pen_list_entries$index1[hh]
+                val[i1] <- val[i1] + der_z[hh]
+                i2 <- diffpar_pen_list_entries$index2[hh]
+                val[i2] <- val[i2] - der_z[hh]
+            }
+        }
+        res$pen_diffpar_lp <- val
+    }
+
+    #*** sum all penalties
+    res$pen_all <- res$pen_prior - res$pen_l2 - res$pen_lp - res$pen_difflp -
+                            res$pen_diffpar_lp
 
     #--- output
     return(res)
diff --git a/R/mgsem_evaluate_penalties_evaluate_entry.R b/R/mgsem_evaluate_penalties_evaluate_entry.R
index fdff3cc..9bf3744 100644
--- a/R/mgsem_evaluate_penalties_evaluate_entry.R
+++ b/R/mgsem_evaluate_penalties_evaluate_entry.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_evaluate_penalties_evaluate_entry.R
-## File Version: 0.188
+## File Version: 0.189
 
 
 mgsem_evaluate_penalties_evaluate_entry <- function(x, res, dd, index, partable,
@@ -16,15 +16,15 @@ mgsem_evaluate_penalties_evaluate_entry <- function(x, res, dd, index, partable,
         null_vec <- 0
     }
 
-    # res[["pen_l2"]] <- null_vec
+    # res[['pen_l2']] <- null_vec
 
-    # fun_eval="none" produces a value of zero
+    # fun_eval='none' produces a value of zero
 
     #*** priors
     if (technical$is_prior){
-        pen_entry <- "pen_prior"
-        prior <- partable[dd,"prior"]
-        if (prior!="none"){
+        pen_entry <- 'pen_prior'
+        prior <- partable[dd,'prior']
+        if (prior!='none'){
             fun_eval <- prior_list[[prior]]
             args_eval <- list()
             val <- mgsem_evaluate_penalties_evaluate_entry_fun_eval(x=x,
diff --git a/R/mgsem_ginv.R b/R/mgsem_ginv.R
new file mode 100644
index 0000000..4c04fe9
--- /dev/null
+++ b/R/mgsem_ginv.R
@@ -0,0 +1,10 @@
+## File Name: mgsem_ginv.R
+## File Version: 0.04
+
+mgsem_ginv <- function(X)
+{
+    requireNamespace('MASS')
+    res <- MASS::ginv(X=X)
+    rownames(res) <- colnames(res) <- colnames(X)
+    return(res)
+}
diff --git a/R/mgsem_grad_fun.R b/R/mgsem_grad_fun.R
index 753b035..6f97e86 100644
--- a/R/mgsem_grad_fun.R
+++ b/R/mgsem_grad_fun.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_grad_fun.R
-## File Version: 0.171
+## File Version: 0.172
 
 
 mgsem_grad_fun <- function(x, opt_fun_args, output_all=FALSE)
@@ -25,7 +25,7 @@ mgsem_grad_fun <- function(x, opt_fun_args, output_all=FALSE)
     for (gg in 1:G){
         grad_suffstat_fun_args <- list(suffstat=opt_fun_args$suffstat[[gg]],
                                         Mu=implied0[[gg]])
-        if (estimator=="ME"){
+        if (estimator=='ME'){
             grad_suffstat_fun_args$p <- opt_fun_args$p
             grad_suffstat_fun_args$eps <- opt_fun_args$eps_approx
             grad_suffstat_fun_args$deriv <- TRUE
diff --git a/R/mgsem_loglike_data.R b/R/mgsem_loglike_data.R
index 1a6779f..fc42de8 100644
--- a/R/mgsem_loglike_data.R
+++ b/R/mgsem_loglike_data.R
@@ -1,10 +1,10 @@
 ## File Name: mgsem_loglike_data.R
-## File Version: 0.02
+## File Version: 0.03
 
 
 mgsem_loglike_data <- function(dat, Mu, Sigma)
 {
-    requireNamespace("mvtnorm")
+    requireNamespace('mvtnorm')
     ll <- sum(mvtnorm::dmvnorm(dat, mean=Mu, sigma=Sigma, log=TRUE))
     return(ll)
 }
diff --git a/R/mgsem_loglike_suffstat.R b/R/mgsem_loglike_suffstat.R
index f186d34..6eb5468 100644
--- a/R/mgsem_loglike_suffstat.R
+++ b/R/mgsem_loglike_suffstat.R
@@ -1,9 +1,8 @@
 ## File Name: mgsem_loglike_suffstat.R
-## File Version: 0.135
+## File Version: 0.136
 
 mgsem_loglike_suffstat <- function(suffstat, Mu, Sigma, output_all=FALSE )
 {
-    requireNamespace("MASS")
     N <- suffstat$N
     M <- suffstat$M
     S <- suffstat$S
@@ -12,7 +11,7 @@ mgsem_loglike_suffstat <- function(suffstat, Mu, Sigma, output_all=FALSE )
         Mu <- res$Mu
         Sigma <- res$Sigma
     }
-    S1 <- MASS::ginv(Sigma)
+    S1 <- mgsem_ginv(X=Sigma)
     p <- length(Mu)
     m1 <- M-Mu
 
diff --git a/R/mgsem_loglike_suffstat_derivative.R b/R/mgsem_loglike_suffstat_derivative.R
index db7b18b..30c6946 100644
--- a/R/mgsem_loglike_suffstat_derivative.R
+++ b/R/mgsem_loglike_suffstat_derivative.R
@@ -1,10 +1,9 @@
 ## File Name: mgsem_loglike_suffstat_derivative.R
-## File Version: 0.16
+## File Version: 0.171
 
 
 mgsem_loglike_suffstat_derivative <- function(suffstat, Mu, Sigma )
 {
-    requireNamespace("MASS")
     N <- suffstat$N
     M <- suffstat$M
     S <- suffstat$S
@@ -13,7 +12,7 @@ mgsem_loglike_suffstat_derivative <- function(suffstat, Mu, Sigma )
         Mu <- res$Mu
         Sigma <- res$Sigma
     }
-    S1 <- MASS::ginv(Sigma)
+    S1 <- mgsem_ginv(X=Sigma)
     p <- length(Mu)
     m1 <- M-Mu
 
diff --git a/R/mgsem_moments_derivative_parameter.R b/R/mgsem_moments_derivative_parameter.R
index 9fe839a..55ef1ac 100644
--- a/R/mgsem_moments_derivative_parameter.R
+++ b/R/mgsem_moments_derivative_parameter.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_moments_derivative_parameter.R
-## File Version: 0.393
+## File Version: 0.395
 
 mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
         i1, i2, h, is_B, eps=1e-12, num_approx=FALSE)
@@ -11,11 +11,11 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
     I <- nrow(est$LAM)
     log0 <- Sigma_der_logical <- matrix(FALSE, nrow=I, ncol=I)
 
-    if (type %in% c("PHI","PSI")){
+    if (type %in% c('PHI','PSI')){
         symm <- TRUE
         calc_Mu <- FALSE
     }
-    if (type %in% c("NU","ALPHA")){
+    if (type %in% c('NU','ALPHA')){
         calc_Sigma <- FALSE
     }
 
@@ -23,7 +23,7 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
     if (is_B){
         num_approx <- TRUE
     }
-    if (type %in% c("LAM","PHI")){
+    if (type %in% c('LAM','PHI')){
         num_approx <- FALSE
     }
 
@@ -32,19 +32,19 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
 
     num_approx <- TRUE
 
-    if (type %in% c("ALPHA","NU","PSI") ) {
+    if (type %in% c('ALPHA','NU','PSI') ) {
         num_approx <- FALSE
     }
 
 
     # issue is for LAM and PHI
-    if (type %in% c("PHI") ) {
+    if (type %in% c('PHI') ) {
         if (i1==i2){
             num_approx <- FALSE
         }
     }
 
-    if (type %in% c("LAM") ) {
+    if (type %in% c('LAM') ) {
         # num_approx <- FALSE
     }
 
@@ -76,14 +76,14 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
         }
         if (calc_Sigma){
             Sigma_der <- (implied1$Sigma-implied2$Sigma)/(2*h)
-            if (type %in% c("B","PHI")){
+            if (type %in% c('B','PHI')){
                 Sigma_der_logical <- mgsem_differ_from_zero(x=Sigma_der, eps=eps)
             }
-            if (type %in% c("PSI")){
+            if (type %in% c('PSI')){
                 Sigma_der_logical[i1,i2] <- TRUE
                 Sigma_der_logical[i2,i1] <- TRUE
             }
-            if (type %in% c("LAM")){
+            if (type %in% c('LAM')){
                 Sigma_der_logical[i1,] <- TRUE
                 Sigma_der_logical[,i1] <- TRUE
             }
@@ -100,20 +100,20 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
         Sigma_der <- mat0
         Sigma_der_logical <- log0
 
-        if (type=="ALPHA"){
+        if (type=='ALPHA'){
             Mu_der <- est$LAM[,i1,drop=FALSE]
         }
-        if (type=="NU"){
+        if (type=='NU'){
             Mu_der[i1,1] <- 1
         }
 
-        if (type=="PSI"){
+        if (type=='PSI'){
             Sigma_der <- mgsem_add_increment(x=mat0,h=1,i1=i1,i2=i2, symm=symm)
             Sigma_der_logical[i1,i2] <- TRUE
             Sigma_der_logical[i2,i1] <- TRUE
         }
 
-        if (type=="LAM"){
+        if (type=='LAM'){
             Mu_der[i1,1] <- est$ALPHA[i2,1]
             # D1 <- matrix(0, nrow=I, ncol=ncol(est$LAM))
             # D1[i1,i2] <- 1
@@ -137,7 +137,7 @@ mgsem_moments_derivative_parameter <- function(est, est_add=NULL, type,
             Sigma_der_logical[i1,] <- TRUE
             Sigma_der_logical[,i1] <- TRUE
         }
-        if (type=="PHI"){
+        if (type=='PHI'){
             # D <- nrow(est$PHI)
             # D1 <- matrix(0, nrow=D, ncol=D)
             # D1[i1,i2] <- 1
diff --git a/R/mgsem_observed_information.R b/R/mgsem_observed_information.R
index ef27cd1..28e183f 100644
--- a/R/mgsem_observed_information.R
+++ b/R/mgsem_observed_information.R
@@ -1,17 +1,20 @@
 ## File Name: mgsem_observed_information.R
-## File Version: 0.13
+## File Version: 0.148
 
 
-mgsem_observed_information <- function(coef, opt_fun_args, technical, comp_se=TRUE)
+mgsem_observed_information <- function(coef, opt_fun_args, technical, comp_se=TRUE,
+        comp_se_me=FALSE)
 {
-    if (technical$estimator!="ML"){
+    is_ml <- technical$estimator=='ML'
+    if (! is_ml){
         comp_se <- FALSE
     }
     res <- list(coef=coef, comp_se=comp_se, se=NA*coef,
-                    vcov=NULL, info_loglike=NULL, info_loglike_pen=NULL)
+                    vcov=NULL, info_loglike=NULL, info_loglike_pen=NULL,
+                    is_ml=is_ml, comp_se_me=comp_se_me)
 
     if (comp_se){
-        requireNamespace("MASS")
+
         grad_fun <- 'mgsem_numerical_gradient'
         h <- technical$h
 
@@ -19,7 +22,7 @@ mgsem_observed_information <- function(coef, opt_fun_args, technical, comp_se=TR
         opt_fun_args1 <- opt_fun_args
         opt_fun_args1$use_penalty <- FALSE
         args_grad <- list(par=coef, FUN=mgsem_grad_fun, symmetrize=TRUE,
-                    h=h, opt_fun_args=opt_fun_args1, output_all=FALSE)
+                                h=h, opt_fun_args=opt_fun_args1, output_all=FALSE)
         res$info_loglike <- do.call(what=grad_fun, args=args_grad)
 
         #- observed information with penalty
@@ -27,8 +30,8 @@ mgsem_observed_information <- function(coef, opt_fun_args, technical, comp_se=TR
         res$info_loglike_pen <- do.call(what=grad_fun, args=args_grad)
 
         #- standard errors
-        res$vcov <- MASS::ginv(X=res$info_loglike)
-        res$se <- sqrt(diag(res$vcov))
+        res$vcov <- mgsem_ginv(X=res$info_loglike)
+        res$se <- mgsem_sqrt_diag(x=res$vcov)
 
     }
 
diff --git a/R/mgsem_opt_fun.R b/R/mgsem_opt_fun.R
index 835a4f5..6c1854f 100644
--- a/R/mgsem_opt_fun.R
+++ b/R/mgsem_opt_fun.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_opt_fun.R
-## File Version: 0.222
+## File Version: 0.224
 
 
 mgsem_opt_fun <- function(x, opt_fun_args, output_all=FALSE)
@@ -35,14 +35,14 @@ mgsem_opt_fun <- function(x, opt_fun_args, output_all=FALSE)
         eval_args <- list(suffstat=opt_fun_args$suffstat[[gg]],
                             Mu=implied$Mu, Sigma=implied$Sigma,
                             output_all=output_all)
-        if (estimator=="ME"){
+        if (estimator=='ME'){
             eval_args$p <- opt_fun_args$p_me
             eval_args$eps <- opt_fun_args$eps_approx
             eval_args$deriv <- FALSE
             eval_args$approx_method <- opt_fun_args$technical$approx_method
         }
         ll_gg <- do.call(what=opt_fun_args$eval_fun, args=eval_args)
-        if (output_all & (estimator=="ML") ){
+        if (output_all & (estimator=='ML') ){
             S1_list[[gg]] <- ll_gg$S1
             mean_residual_list[[gg]] <- ll_gg$mean_residual
             ll_gg <- ll_gg$loglike
diff --git a/R/mgsem_output_proc_casewise_likelihood.R b/R/mgsem_output_proc_casewise_likelihood.R
index 544ae02..2e2ea77 100644
--- a/R/mgsem_output_proc_casewise_likelihood.R
+++ b/R/mgsem_output_proc_casewise_likelihood.R
@@ -1,12 +1,12 @@
 ## File Name: mgsem_output_proc_casewise_likelihood.R
-## File Version: 0.06
+## File Version: 0.07
 
 
 mgsem_output_proc_casewise_likelihood <- function(data_proc, implied, estimator="ML")
 {
     case_ll <- NULL
-    if (estimator=="ML" & ( ! is.null(data_proc) ) ){
-        requireNamespace("mvtnorm")
+    if (estimator=='ML' & ( ! is.null(data_proc) ) ){
+        requireNamespace('mvtnorm')
         N <- data_proc$N
         G <- data_proc$G
         data <- data_proc$data
diff --git a/R/mgsem_partable2coef.R b/R/mgsem_partable2coef.R
index 54e1951..f833672 100644
--- a/R/mgsem_partable2coef.R
+++ b/R/mgsem_partable2coef.R
@@ -1,11 +1,11 @@
 ## File Name: mgsem_partable2coef.R
-## File Version: 0.04
+## File Version: 0.05
 
 mgsem_partable2coef <- function(partable)
 {
     dfr <- partable
     ind <- which(dfr$unique==1)
-    coef <- dfr[ ind,"est"]
-    names(coef) <- dfr[ind,"name"]
+    coef <- dfr[ ind,'est']
+    names(coef) <- dfr[ind,'name']
     return(coef)
 }
diff --git a/R/mgsem_partable2model.R b/R/mgsem_partable2model.R
index fca60b8..fc37656 100644
--- a/R/mgsem_partable2model.R
+++ b/R/mgsem_partable2model.R
@@ -1,22 +1,22 @@
 ## File Name: mgsem_partable2model.R
-## File Version: 0.14
+## File Version: 0.151
 
 
 mgsem_partable2model <- function(partable, model, index=FALSE)
 {
     ND <- nrow(partable)
-    entries <- c("est")
+    entries <- c('est')
     if (index){
-        entries <- c("est","index")
+        entries <- c('est','index')
     }
     for (entry in entries){
         for (dd in 1:ND){
-            hh <- partable[dd,"group"]+1
-            type <- paste(partable[dd,"type"])
+            hh <- partable[dd,'group']+1
+            type <- paste(partable[dd,'type'])
             mat <- model[[hh]][[entry]][[type]]
-            mat[ partable[dd,"i1"], partable[dd,"i2"] ] <- partable[dd,entry]
-            if (type %in% c("PHI","PSI") ){
-                mat[ partable[dd,"i2"], partable[dd,"i1"] ] <- partable[dd,entry]
+            mat[ partable[dd,'i1'], partable[dd,'i2'] ] <- partable[dd,entry]
+            if (type %in% c('PHI','PSI') ){
+                mat[ partable[dd,'i2'], partable[dd,'i1'] ] <- partable[dd,entry]
             }
             model[[hh]][[entry]][[ type ]] <- mat
         }
diff --git a/R/mgsem_power_fun_differentiable_approx.R b/R/mgsem_power_fun_differentiable_approx.R
index b17b29a..41c5bb1 100644
--- a/R/mgsem_power_fun_differentiable_approx.R
+++ b/R/mgsem_power_fun_differentiable_approx.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_power_fun_differentiable_approx.R
-## File Version: 0.08
+## File Version: 0.091
 
 mgsem_power_fun_differentiable_approx <- function(x, p, eps, deriv=FALSE,
         approx_method="lp")
@@ -7,7 +7,7 @@ mgsem_power_fun_differentiable_approx <- function(x, p, eps, deriv=FALSE,
     # logcomp <- TRUE
     logcomp <- FALSE
     if (deriv){
-        if (approx_method=="lp"){
+        if (approx_method=='lp'){
             if (!logcomp){
                 res <- p*((x^2+eps)^(p/2-1))*x
             } else {
@@ -16,14 +16,14 @@ mgsem_power_fun_differentiable_approx <- function(x, p, eps, deriv=FALSE,
                 res <- x*exp( p2*log(x^2+eps) + logp )
             }
         }
-        if (approx_method=="l2"){
+        if (approx_method=='l2'){
             res <- 2*x
         }
     } else {
-        if (approx_method=="lp"){
+        if (approx_method=='lp'){
             res <- (x^2+eps)^(p/2)
         }
-        if (approx_method=="l2"){
+        if (approx_method=='l2'){
             res <- x^2
         }
     }
diff --git a/R/mgsem_proc_data.R b/R/mgsem_proc_data.R
index 13680fd..c80914c 100644
--- a/R/mgsem_proc_data.R
+++ b/R/mgsem_proc_data.R
@@ -1,11 +1,14 @@
 ## File Name: mgsem_proc_data.R
-## File Version: 0.04
+## File Version: 0.055
 
 mgsem_proc_data <- function(data, group, weights)
 {
+    N <- nrow(data)
+    if (is.null(group)){
+        group <- rep(1,N)
+    }
     groups <- unique(group)
     G <- length(groups)
-    N <- nrow(data)
     if (is.null(weights)){
         weights <- rep(1,N)
     }
@@ -15,7 +18,7 @@ mgsem_proc_data <- function(data, group, weights)
         ind_gg <- which(group==groups[gg])
         dat_gg <- data[ ind_gg, ]
         w <- weights[ind_gg]
-        res <- stats::cov.wt( x=dat_gg, w=w, method="ML")
+        res <- stats::cov.wt( x=dat_gg, w=w, method='ML')
         s_gg <- list(M=res$center, S=res$cov, N=sum(w) )
         suffstat[[gg]] <- s_gg
     }
diff --git a/R/mgsem_proc_model.R b/R/mgsem_proc_model.R
index caf80fb..66980e1 100644
--- a/R/mgsem_proc_model.R
+++ b/R/mgsem_proc_model.R
@@ -1,27 +1,27 @@
 ## File Name: mgsem_proc_model.R
-## File Version: 0.286
+## File Version: 0.298
 
 mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
         prior_list=NULL, pen_type="lasso", fixed_parms=FALSE,
-        partable_start=NULL)
+        partable_start=NULL, diffpar_pen=NULL)
 {
 
     dfr <- NULL
-    types <- c("ALPHA", "NU", "LAM", "PHI", "PSI")
-    symm_types <- c( "PHI", "PSI")
+    types <- c('ALPHA', 'NU', 'LAM', 'PHI', 'PSI')
+    symm_types <- c( 'PHI', 'PSI')
     N <- sum(N_group)
 
     names_prior_list <- names(prior_list)
     is_B <- mgsem_proc_model_is_B(model=model)
     technical$is_B <- is_B
     if (is_B){
-        types <- c(types, "B")
+        types <- c(types, 'B')
     }
 
-    I <- mgsem_proc_model_extract_dimension(model=model, entry="est",
-                            type="LAM", nrow=TRUE)
-    D <- mgsem_proc_model_extract_dimension(model=model, entry="est",
-                            type="LAM", nrow=FALSE)
+    I <- mgsem_proc_model_extract_dimension(model=model, entry='est',
+                            type='LAM', nrow=TRUE)
+    D <- mgsem_proc_model_extract_dimension(model=model, entry='est',
+                            type='LAM', nrow=FALSE)
 
     #** process case of single model
     model <- mgsem_proc_model_single_group(model=model)
@@ -34,7 +34,7 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
 
         #-- include missing entries
         model[[hh]] <- mgsem_proc_model_include_missing_entries(model_hh=model[[hh]],
-                            types=types, entries=c("est","index"), I=I, D=D)
+                            types=types, entries=c('est','index'), I=I, D=D)
         model_hh <- model[[hh]]
         est <- model_hh$est
         index <- model_hh$index
@@ -65,9 +65,9 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
                             dfr1 <- data.frame( type=type, i1=ii, i2=jj,
                                                 group=group)
                             dfr1$name <- paste0(dfr1$type, dfr1$i1, dfr1$i2,
-                                                "_G", dfr1$group)
-                            dfr1$name2 <- paste0(dfr1$type, dfr1$i1, "-", dfr1$i2,
-                                                "_G", dfr1$group)
+                                                '_G', dfr1$group)
+                            dfr1$name2 <- paste0(dfr1$type, dfr1$i1, '-', dfr1$i2,
+                                                '_G', dfr1$group)
                             symm <- symm0
                             if (ii==jj){
                                 symm <- FALSE
@@ -85,8 +85,8 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
 
                             #-- check for entries
                             #-- model specifications
-                            entries <- c("lower", "upper","prior", "pen_l2",
-                                            "pen_lp", "pen_difflp")
+                            entries <- c('lower', 'upper','prior', 'pen_l2',
+                                            'pen_lp', 'pen_difflp')
                             dfr1 <- mgsem_proc_model_add_specs_all(model=model_hh,
                                             entries=entries, type=type, ii=ii, jj=jj,
                                             dfr1=dfr1, names_prior_list=names_prior_list,
@@ -122,15 +122,34 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
     technical$is_pen_difflp <- difflp_info$is_pen_difflp
 
     #*** specifications in technical
-    technical$is_prior <- sum(dfr$prior!="none") > 0
+    technical$is_prior <- sum(dfr$prior!='none') > 0
     technical$is_pen_l2 <- sum(dfr$pen_l2>0) > 0
     technical$is_pen_lp <- sum(dfr$pen_lp>0) > 0
 
+    technical$is_diffpar_pen <- ! is.null( diffpar_pen )
+    if (technical$is_diffpar_pen){
+        coef_indices <- which( rowSums( diffpar_pen$W ) > 0 )
+        diffpar_pen$coef_indices <- coef_indices
+        dp1 <- NULL
+        NW <- ncol(W)
+        for (ww in 1:NW){
+            for (uu in 1:NW){
+                val <- W[ww,uu]
+                if (abs(val) > 1e-14){
+                    dp2 <- data.frame(index1=ww, index2=uu, W=val)
+                    dp1 <- rbind(dp1, dp2)
+                }
+            }
+        }
+        diffpar_pen$diffpar_pen_list_entries <- dp1
+        technical$diffpar_pen <- diffpar_pen
+    }
+
     dfr <- as.data.frame(dfr)
 
     #** coefficient vector of estimated parameters
     if ( ! is.null(partable_start) ){
-        dfr$start <- dfr$est <- partable_start
+        dfr$start <- dfr$est <- partable_start$est
     }
     coef <- mgsem_partable2coef(partable=dfr)
 
@@ -153,7 +172,7 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group,
     loop_parms <- (1:ND)[ dfr$unique==1]
 
     #- rewrite penalty parameters into model matrices
-    entries <- c("pen_l2", "pen_lp", "pen_difflp")
+    entries <- c('pen_l2', 'pen_lp', 'pen_difflp')
     model <- mgsem_proc_model_update_penalties_matrix(partable=dfr,
                     entries=entries, model=model)
 
diff --git a/R/mgsem_proc_model_add_specs_all.R b/R/mgsem_proc_model_add_specs_all.R
index 64f90a4..9bf4a69 100644
--- a/R/mgsem_proc_model_add_specs_all.R
+++ b/R/mgsem_proc_model_add_specs_all.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_add_specs_all.R
-## File Version: 0.096
+## File Version: 0.097
 
 
 mgsem_proc_model_add_specs_all <- function(model, entries, type, ii, jj, dfr1,
@@ -12,32 +12,32 @@ mgsem_proc_model_add_specs_all <- function(model, entries, type, ii, jj, dfr1,
     N1 <- 1
     NE <- length(entries)
     for (entry in entries){
-        if (entry=="lower"){ default <- -Inf }
-        if (entry=="upper"){ default <- Inf }
-        if (entry=="prior"){ default <- "none" }
-        if (entry=="pen_l2"){ default <- 0 }
-        if (entry=="pen_lp"){ default <- 0 }
-        if (entry=="pen_difflp"){ default <- 0 }
+        if (entry=='lower'){ default <- -Inf }
+        if (entry=='upper'){ default <- Inf }
+        if (entry=='prior'){ default <- 'none' }
+        if (entry=='pen_l2'){ default <- 0 }
+        if (entry=='pen_lp'){ default <- 0 }
+        if (entry=='pen_difflp'){ default <- 0 }
 
         val <- mgsem_proc_model_add_specs(model=model, entry=entry, type=type,
                     ii=ii, jj=jj, default=default)
-        if (entry=="prior"){ val <- paste(val) }
-        if (entry %in% c("pen_l2","pen_lp","pen_difflp")){
+        if (entry=='prior'){ val <- paste(val) }
+        if (entry %in% c('pen_l2','pen_lp','pen_difflp')){
             val <- N1*val
-            if ( (entry=="pen_difflp") & (group==0)){
+            if ( (entry=='pen_difflp') & (group==0)){
                 val <- 0
             }
             #- no penalty
-            if (pen_type=="none"){
+            if (pen_type=='none'){
                 val <- 0
             }
         }
 
         #- correctness checks
-        if (entry=="prior"){
-            if ( ( val!="none") & ( ! val %in% names_prior_list ) ){
-                v1 <- paste0( "Specified prior in ", type, "[",
-                        ii, ",", jj, "] in group ", group, " not in 'prior_list'!\n")
+        if (entry=='prior'){
+            if ( ( val!='none') & ( ! val %in% names_prior_list ) ){
+                v1 <- paste0( 'Specified prior in ', type, '[',
+                        ii, ',', jj, '] in group ', group, ' not in \'prior_list\'!\n')
                 stop(v1)
             }
         }
diff --git a/R/mgsem_proc_model_difflp_information.R b/R/mgsem_proc_model_difflp_information.R
index 76cec9a..0478b47 100644
--- a/R/mgsem_proc_model_difflp_information.R
+++ b/R/mgsem_proc_model_difflp_information.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_difflp_information.R
-## File Version: 0.186
+## File Version: 0.187
 
 
 mgsem_proc_model_difflp_information <- function(partable)
@@ -10,11 +10,11 @@ mgsem_proc_model_difflp_information <- function(partable)
     lpdiff_diff_indices <- list()
     pen_difflp_factor <- list()
 
-    pen_entry <- "pen_difflp"
+    pen_entry <- 'pen_difflp'
     loop_parms <- which(partable$unique==1 & partable$group>0)
 
     for (dd in loop_parms){
-        index_dd <- partable[dd,"index"]
+        index_dd <- partable[dd,'index']
         ind <- which( partable$type==partable$type[dd] & partable$i1==partable$i1[dd] &
                 partable$i2==partable$i2[dd] & (partable$group > 0) & (partable$rid!=dd) &
                 partable$unique==1)
@@ -23,7 +23,7 @@ mgsem_proc_model_difflp_information <- function(partable)
         }
         if ( (length(ind)>0) & (partable$pen_difflp[dd]>0)  ){
             is_lpdiff_entry[dd] <- TRUE
-            lpdiff_diff_indices[[dd]] <- partable[ind, "index"]
+            lpdiff_diff_indices[[dd]] <- partable[ind, 'index']
             fac <- sqrt(partable[dd,pen_entry]*partable[ind,pen_entry])
             pen_difflp_factor[[dd]] <- fac
         } else {
@@ -32,17 +32,17 @@ mgsem_proc_model_difflp_information <- function(partable)
     }
 
     #--- design matrices
-    cp <- partable[ partable$unique==1, "name"]
+    cp <- partable[ partable$unique==1, 'name']
     lpdiff_facmat <- matrix(0, nrow=NP, ncol=NP)
     rownames(lpdiff_facmat) <- colnames(lpdiff_facmat) <- cp
     lpdiff_n <- lpdiff_facmat
     for (dd in loop_parms){
         if (is_lpdiff_entry[dd]){
-            i1 <- partable[dd,"index"]
+            i1 <- partable[dd,'index']
             i2 <- lpdiff_diff_indices[[dd]]
             lpdiff_facmat[i1,i2] <- pen_difflp_factor[[dd]]
-            n1 <- partable[partable$index==i1,"N_group"]
-            n2 <- partable[partable$index%in%i2,"N_group"]
+            n1 <- partable[partable$index==i1,'N_group']
+            n2 <- partable[partable$index%in%i2,'N_group']
             lpdiff_n[i1,i2] <- 0.5*sqrt(n1*n2)
         }
     }
diff --git a/R/mgsem_proc_model_include_missing_entries.R b/R/mgsem_proc_model_include_missing_entries.R
index 3759a32..9037638 100644
--- a/R/mgsem_proc_model_include_missing_entries.R
+++ b/R/mgsem_proc_model_include_missing_entries.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_include_missing_entries.R
-## File Version: 0.09
+## File Version: 0.11
 
 
 mgsem_proc_model_include_missing_entries <- function(model_hh, types,
@@ -7,12 +7,12 @@ mgsem_proc_model_include_missing_entries <- function(model_hh, types,
 {
     for (tt in types){
         val <- 0
-        if (tt=="NU"){ NR <- I; NC <- 1 }
-        if (tt=="ALPHA"){ NR <- D; NC <- 1 }
-        if (tt=="LAM"){ NR <- I; NC <- D }
-        if (tt=="PHI"){ NR <- D; NC <- D }
-        if (tt=="PSI"){ NR <- I; NC <- I }
-        if (tt=="B"){ NR <- D; NC <- D }
+        if (tt=='NU'){ NR <- I; NC <- 1 }
+        if (tt=='ALPHA'){ NR <- D; NC <- 1 }
+        if (tt=='LAM'){ NR <- I; NC <- D }
+        if (tt=='PHI'){ NR <- D; NC <- D }
+        if (tt=='PSI'){ NR <- I; NC <- I }
+        if (tt=='B'){ NR <- D; NC <- D }
         for (ee in entries){
             if (is.null(model_hh[[ee]][[tt]])){
                 model_hh[[ee]][[tt]] <- matrix(val, nrow=NR, ncol=NC)
diff --git a/R/mgsem_proc_model_is_B.R b/R/mgsem_proc_model_is_B.R
index ef8ca46..7de5723 100644
--- a/R/mgsem_proc_model_is_B.R
+++ b/R/mgsem_proc_model_is_B.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_is_B.R
-## File Version: 0.02
+## File Version: 0.03
 
 
 mgsem_proc_model_is_B <- function(model)
@@ -7,7 +7,7 @@ mgsem_proc_model_is_B <- function(model)
     H <- length(model)
     is_B <- 0
     for (hh in 1:H){
-        is_B <- is_B + ( ! is.null( model[[hh]][["est"]][["B"]] ) )
+        is_B <- is_B + ( ! is.null( model[[hh]][['est']][['B']] ) )
     }
     is_B <- ( is_B > 0 )
     return(is_B)
diff --git a/R/mgsem_proc_model_partable_define_index.R b/R/mgsem_proc_model_partable_define_index.R
index 7bbd1ba..d9f5350 100644
--- a/R/mgsem_proc_model_partable_define_index.R
+++ b/R/mgsem_proc_model_partable_define_index.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_partable_define_index.R
-## File Version: 0.11
+## File Version: 0.121
 
 
 mgsem_proc_model_partable_define_index <- function(partable)
@@ -16,26 +16,26 @@ mgsem_proc_model_partable_define_index <- function(partable)
     #** find duplicated parameters
     ND <- nrow(dfr)
     for (dd in 1:ND){
-        if (dfr[dd,"unique"]==0){
-            i1 <- which( dfr[,"index"]==dfr[dd,"index"])
-            i2 <- which( dfr[,"unique"]==1)
+        if (dfr[dd,'unique']==0){
+            i1 <- which( dfr[,'index']==dfr[dd,'index'])
+            i2 <- which( dfr[,'unique']==1)
             i3 <- intersect(i1,i2)
-            dfr[dd,"name"] <- dfr[i3,"name"]
+            dfr[dd,'name'] <- dfr[i3,'name']
         }
         #- recycled parameters
-        group_dd <- dfr[dd,"group"]
+        group_dd <- dfr[dd,'group']
         if (group_dd>0){
-            name_dd <- dfr[dd,"name"]
-            name0_dd <- gsub(paste0("_G",group_dd),"_G0", name_dd)
+            name_dd <- dfr[dd,'name']
+            name0_dd <- gsub(paste0('_G',group_dd),'_G0', name_dd)
             ind <- match(name0_dd, dfr$name)
             if (length(ind)==1){
-                dfr[dd,"recycle"] <- ind
+                dfr[dd,'recycle'] <- ind
             }
         }
     }
     dfr$recycle <- ifelse(is.na(dfr$recycle), 0, dfr$recycle )
 
-    dfr <- data.frame("rid"=1:ND, dfr)
+    dfr <- data.frame(rid=1:ND, dfr)
     NP <- max(dfr$index)
     #--- output
     res <- list(partable=dfr, ND=ND, NP=NP)
diff --git a/R/mgsem_proc_model_single_group.R b/R/mgsem_proc_model_single_group.R
index 3b63a32..511ca42 100644
--- a/R/mgsem_proc_model_single_group.R
+++ b/R/mgsem_proc_model_single_group.R
@@ -1,11 +1,11 @@
 ## File Name: mgsem_proc_model_single_group.R
-## File Version: 0.07
+## File Version: 0.081
 
 mgsem_proc_model_single_group <- function(model)
 {
     H <- length(model)
     if (H==1){
-        entries <- c("est","index")
+        entries <- c('est','index')
         group1 <- list()
         group0 <- model[[1]]
         for (ee in entries){
diff --git a/R/mgsem_proc_model_update_penalties_matrix.R b/R/mgsem_proc_model_update_penalties_matrix.R
index 4dd1e6c..cc31e14 100644
--- a/R/mgsem_proc_model_update_penalties_matrix.R
+++ b/R/mgsem_proc_model_update_penalties_matrix.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_model_update_penalties_matrix.R
-## File Version: 0.071
+## File Version: 0.072
 
 mgsem_proc_model_update_penalties_matrix <- function(partable, entries, model)
 {
@@ -7,7 +7,7 @@ mgsem_proc_model_update_penalties_matrix <- function(partable, entries, model)
     ND <- nrow(partable)
     for (entry in entries){
         for (dd in 1:ND){
-            if (partable[dd,"unique"]==1){
+            if (partable[dd,'unique']==1){
                 group <- partable$group[dd]+1
                 type_dd <- paste(partable$type[dd])
                 mat_gg_ee <- model[[group]][[entry]]
diff --git a/R/mgsem_proc_suffstat.R b/R/mgsem_proc_suffstat.R
index fd83b10..bd1163f 100644
--- a/R/mgsem_proc_suffstat.R
+++ b/R/mgsem_proc_suffstat.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_proc_suffstat.R
-## File Version: 0.097
+## File Version: 0.098
 
 
 mgsem_proc_suffstat <- function(suffstat)
@@ -11,22 +11,22 @@ mgsem_proc_suffstat <- function(suffstat)
         N_group[gg] <- suffstat_gg$N
         # rename entries
         ns_gg <- names(suffstat[[gg]])
-        if( ! ( "M" %in% ns_gg ) ){
+        if( ! ( 'M' %in% ns_gg ) ){
             ind <- min( which( unlist( lapply(suffstat_gg, is.vector) ) ))
-            ns_gg[ind] <- "M"
+            ns_gg[ind] <- 'M'
         }
-        if( ! ( "S" %in% ns_gg ) ){
+        if( ! ( 'S' %in% ns_gg ) ){
             ind <- min( which( unlist( lapply(suffstat_gg, is.matrix) ) ))
-            ns_gg[ind] <- "S"
+            ns_gg[ind] <- 'S'
         }
         names(suffstat_gg) <- ns_gg
         I <- length(suffstat_gg$M)
 
         #*** weights for moment-based estimation
-        if ( ! ( "weights_M" %in% ns_gg ) ){
+        if ( ! ( 'weights_M' %in% ns_gg ) ){
             suffstat_gg$weights_M <- 1+0*suffstat_gg$M
         }
-        if ( ! ( "weights_S" %in% ns_gg ) ){
+        if ( ! ( 'weights_S' %in% ns_gg ) ){
             suffstat_gg$weights_S <- 1+0*suffstat_gg$S
         }
         suffstat_gg$vech_weights_S <- mgsem_vech(x=suffstat_gg$weights_S)
diff --git a/R/mgsem_proc_technical.R b/R/mgsem_proc_technical.R
index b08d393..8bb1626 100644
--- a/R/mgsem_proc_technical.R
+++ b/R/mgsem_proc_technical.R
@@ -1,46 +1,46 @@
 ## File Name: mgsem_proc_technical.R
-## File Version: 0.172
+## File Version: 0.174
 
 mgsem_proc_technical <- function(technical, control, p_me, p_pen, eps_approx, suffstat,
-            estimator)
+            estimator, diffpar_pen=NULL)
 {
     #*** defaults technical; update list
-    technical0 <- list(maxiter=1000, h=1e-5, eps_zero=1e-12, optimizer="nlminb",
+    technical0 <- list(maxiter=1000, h=1e-5, eps_zero=1e-12, optimizer='nlminb',
                         use_deriv=TRUE,    num_approx=FALSE,
                         eps_count_penal=1e-2,
-                        approx_method="lp", use_rcpp_penalty=TRUE)
+                        approx_method='lp', use_rcpp_penalty=TRUE)
     use_approx_method <- ! is.null(technical$approx_method)
     technical <- mgsem_update_list_entries(add_list=technical, output_list=technical0)
 
     if ((p_pen==2) & (p_me==2) & ( ! use_approx_method) ){
-        technical$approx_method    <- "l2"
+        technical$approx_method    <- 'l2'
     }
 
-    if (technical$approx_method=="l2" & ( (p_pen!=2) | ( p_me!=2) ) ){
-        technical$approx_method    <- "lp"
+    if (technical$approx_method=='l2' & ( (p_pen!=2) | ( p_me!=2) ) ){
+        technical$approx_method    <- 'lp'
     }
-    # if (technical$approx_method=="l2"){
+    # if (technical$approx_method=='l2'){
     #     eps_approx <- 0
     # }
 
     #- optimizer
-    if (technical$optimizer=="nlminb"){
+    if (technical$optimizer=='nlminb'){
         control$iter.max <- technical$maxiter
     }
-    if (technical$optimizer=="optim"){
+    if (technical$optimizer=='optim'){
         control$maxit <- technical$maxiter
     }
-    if (technical$optimizer=="bobyqa"){
+    if (technical$optimizer=='bobyqa'){
         fac <- max(15, 3*nrow(suffstat[[1]][[2]]))
         technical$maxiter <- technical$maxiter*fac
         control$maxfun <- technical$maxiter
     }
-    if (technical$optimizer=="nloptr"){
+    if (technical$optimizer=='nloptr'){
         control$maxeval <- technical$maxiter
     }
 
     p <- p_pen
-    if (estimator=="ME"){
+    if (estimator=='ME'){
         p <- p_me
     }
 
diff --git a/R/mgsem_sqrt_diag.R b/R/mgsem_sqrt_diag.R
new file mode 100644
index 0000000..e29cdb5
--- /dev/null
+++ b/R/mgsem_sqrt_diag.R
@@ -0,0 +1,7 @@
+## File Name: mgsem_sqrt_diag.R
+## File Version: 0.01
+
+mgsem_sqrt_diag <- function(x)
+{
+    sqrt(diag(x))
+}
diff --git a/R/mgsem_suffstat_covariance_matrix.R b/R/mgsem_suffstat_covariance_matrix.R
new file mode 100644
index 0000000..b2bef27
--- /dev/null
+++ b/R/mgsem_suffstat_covariance_matrix.R
@@ -0,0 +1,60 @@
+## File Name: mgsem_suffstat_covariance_matrix.R
+## File Version: 0.099
+
+mgsem_suffstat_covariance_matrix <- function(suffstat)
+{
+    G <- length(suffstat)
+    dfr <- NULL
+    # gg <- 1
+
+    for (gg in 1:G){
+        suffstat_gg <- suffstat[[gg]]
+        names_gg <- names(suffstat_gg)
+        mu <- suffstat_gg[[ intersect( c('mu','M'), names_gg) ]]
+        Sigma <- suffstat_gg[[ intersect( c('S','Sigma'), names_gg) ]]
+        N <- suffstat_gg$N
+        vars <- names(mu)
+        I1 <- length(mu)
+        if (is.null(vars)){
+            vars <- paste0('V',1:I1)
+        }
+        #** mu
+        label <- paste0('mu_G', gg, '_', vars)
+        se <- sqrt( diag(Sigma) / N )
+        dfr1 <- data.frame(type='mu', group=gg, index1=1:I1, index2=NA, label=label,
+                    par=mu, se=se)
+        V1 <- Sigma / N
+        rownames(V1) <- colnames(V1) <- label
+
+        #** sigma
+        c1 <- rbind( 1:I1, 1:I1)
+        c2 <- utils::combn(x=I1, m=2)
+        c3 <- rbind( t(c1), t(c2))
+        c3 <- c3[ order(c3[,1]), ]
+        label <- paste0('Sigma_G', gg, '_', vars[c3[,1]], '_', vars[c3[,2]] )
+        par <- Sigma[ c3[,1:2] ]
+
+        #* compute duplication matrix
+        K <- mgsem_duplication_matrix(x=Sigma)
+
+        dfr2 <- data.frame(type='sigma', group=gg, index1=c3[,1], index2=c3[,2],
+                        label=label, par=par)
+        N2 <- nrow(dfr2)
+        V2 <- matrix(0, nrow=N2, ncol=N2)
+        V2 <- 2* ( K %*% ( kronecker(X=Sigma, Y=Sigma) ) %*% t(K) ) / N
+        rownames(V2) <- colnames(V2) <- label
+        dfr2$se <- sqrt(diag(V2))
+        V3 <- mgsem_bdiag(x1=V1, x2=V2)
+        dfr <- rbind(dfr, dfr1, dfr2)
+        if (gg==1){
+            V <- V3
+        } else {
+            V <- mgsem_bdiag(x1=V, x2=V3)
+        }
+
+    }    # end gg
+
+    #--- output
+    res <- list(suffstat_pars=dfr, suffstat_vcov=V)
+    return(res)
+}
diff --git a/R/mgsem_test_fun.R b/R/mgsem_test_fun.R
index 86822d8..4399937 100644
--- a/R/mgsem_test_fun.R
+++ b/R/mgsem_test_fun.R
@@ -1,11 +1,11 @@
 ## File Name: mgsem_test_fun.R
-## File Version: 0.15
+## File Version: 0.161
 
 
 mgsem_test_fun <- function(test, coef, opt_fun_args)
 {
     if (test){
-        requireNamespace("miceadds")
+        requireNamespace('miceadds')
         #- function evaluation
         ll <- mgsem_opt_fun(x=coef, opt_fun_args=opt_fun_args)
         #- numerical gradient
@@ -17,8 +17,8 @@ mgsem_test_fun <- function(test, coef, opt_fun_args)
         # dfr <- cbind( grad, grad1)
 
         #- print
-        miceadds::Revalpr("ll")
-        miceadds::Revalpr_maxabs("grad","grad1")
+        miceadds::Revalpr('ll')
+        miceadds::Revalpr_maxabs('grad','grad1')
         stop()
     }
 }
diff --git a/R/mgsem_vcov_me.R b/R/mgsem_vcov_me.R
new file mode 100644
index 0000000..f9609c0
--- /dev/null
+++ b/R/mgsem_vcov_me.R
@@ -0,0 +1,81 @@
+## File Name: mgsem_vcov_me.R
+## File Version: 0.162
+
+mgsem_vcov_me <- function(coef, opt_fun_args, suffstat_vcov, comp_se,
+            se_delta_formula=FALSE)
+{
+    estimator <- opt_fun_args$technical$estimator
+    comp_se_me <- (estimator=='ME') & comp_se
+    if ( (estimator=='ML') & comp_se & se_delta_formula ){
+        comp_se_me <- TRUE
+    }
+    if (comp_se_me){
+        NP <- length(coef)
+        parnames <- names(coef)
+        #-- derivative with respect to coef
+        grad_der_coef <- matrix(0, nrow=NP, ncol=NP)
+        rownames(grad_der_coef) <- colnames(grad_der_coef) <- parnames
+        h <- opt_fun_args$technical$h
+        for (pp in 1:NP){
+            coef1 <- mgsem_add_increment(x=coef, h=h, i1=pp)
+            res1 <- mgsem_grad_fun(x=coef1, opt_fun_args=opt_fun_args, output_all=FALSE)
+            coef2 <- mgsem_add_increment(x=coef, h=-h, i1=pp)
+            res2 <- mgsem_grad_fun(x=coef2, opt_fun_args=opt_fun_args, output_all=FALSE)
+            grad_der_coef[,pp] <- ( res1 - res2 )/(2*h)
+        }
+        grad_der_coef <- ( grad_der_coef + t(grad_der_coef) ) / 2
+
+        #-- derivative with respect input parameters of sufficient statistics
+        suffstat_pars <- suffstat_vcov$suffstat_pars
+        SP <- nrow(suffstat_pars)
+        V <- suffstat_vcov$suffstat_vcov
+        grad_der_suffstat <- matrix(0, nrow=NP, ncol=SP)
+        rownames(grad_der_suffstat) <- parnames
+        colnames(grad_der_suffstat) <- suffstat_pars$label
+        opt_fun_args1 <- opt_fun_args
+        suffstat <- opt_fun_args$suffstat
+        for (pp in 1:SP){
+            suffstat_pars_pp <- suffstat_pars[pp,]
+            group_pp <- suffstat_pars_pp$group
+            val_pp <- list(NA,2)
+            for (oo in 1:2){
+                suffstat1 <- suffstat
+                u <- h
+                if (oo==2){
+                    u <- -h
+                }
+                if (suffstat_pars_pp$type=='mu'){
+                    entry1 <- mgsem_add_increment(x=suffstat1[[group_pp]]$M, h=u,
+                                        i1=suffstat_pars_pp$index1)
+                    suffstat1[[group_pp]]$M <- entry1
+                } else {
+                    entry1 <- mgsem_add_increment(x=suffstat1[[group_pp]]$S, h=u,
+                                        i1=suffstat_pars_pp$index1,
+                                        i2=suffstat_pars_pp$index2, symm=TRUE)
+                    suffstat1[[group_pp]]$S <- entry1
+                }
+                opt_fun_args1$suffstat <- suffstat1
+                val_pp[[oo]] <- mgsem_grad_fun(x=coef, opt_fun_args=opt_fun_args1,
+                                    output_all=FALSE)
+            }
+            der_est <- (val_pp[[1]]-val_pp[[2]])/(2*h)
+
+            grad_der_suffstat[,pp] <- der_est
+        }
+
+        #-- compute transformation matrix A
+        W1 <- mgsem_ginv(X=grad_der_coef)
+        A <- - ( W1 %*% grad_der_suffstat )
+        vcov <- A %*% V %*% t(A)
+        rownames(vcov) <- colnames(vcov) <- names(coef)
+        se <- mgsem_sqrt_diag(x=vcov)
+        names(se) <- names(coef)
+    } else {
+        vcov <- NULL
+        se <- NULL
+    }
+
+    #--- output
+    res <- list(vcov=vcov, se=se, comp_se_me=comp_se_me)
+    return(res)
+}
diff --git a/R/mgsem_vec.R b/R/mgsem_vec.R
new file mode 100644
index 0000000..acba4e0
--- /dev/null
+++ b/R/mgsem_vec.R
@@ -0,0 +1,7 @@
+## File Name: mgsem_vec.R
+## File Version: 0.01
+
+mgsem_vec <- function(x)
+{
+    as.vector(x)
+}
diff --git a/R/mgsem_vech.R b/R/mgsem_vech.R
index a482c16..eafea97 100644
--- a/R/mgsem_vech.R
+++ b/R/mgsem_vech.R
@@ -1,8 +1,8 @@
 ## File Name: mgsem_vech.R
-## File Version: 0.01
+## File Version: 0.02
 
 mgsem_vech <- function(x)
 {
-    res <- x[ ! upper.tri(x)]
+    res <- x[ ! upper.tri(x) ]
     return(res)
 }
diff --git a/R/mi_inv_lavaan_modification_indices.R b/R/mi_inv_lavaan_modification_indices.R
index d6509ef..7ae80a5 100644
--- a/R/mi_inv_lavaan_modification_indices.R
+++ b/R/mi_inv_lavaan_modification_indices.R
@@ -1,9 +1,9 @@
 ## File Name: mi_inv_lavaan_modification_indices.R
-## File Version: 0.05
+## File Version: 0.06
 
 mi_inv_lavaan_modification_indices <- function(mod, op=c("~1","=~"))
 {
-    requireNamespace("lavaan")
+    requireNamespace('lavaan')
     res <- lavaan::modificationIndices(object=mod, free.remove=FALSE,
                     op=op, sort=TRUE)
     return(res)
diff --git a/R/stratified.cronbach.alpha.R b/R/stratified.cronbach.alpha.R
index dac0cd0..f0b845c 100644
--- a/R/stratified.cronbach.alpha.R
+++ b/R/stratified.cronbach.alpha.R
@@ -1,5 +1,5 @@
 ## File Name: stratified.cronbach.alpha.R
-## File Version: 0.253
+## File Version: 0.257
 
 
 # stratified Cronbach's Alpha
@@ -23,7 +23,9 @@ stratified.cronbach.alpha <- function( data, itemstrata=NULL )
     # stratified alpha
     dfr$alpha.stratified <- NA
     var_tot <- dfr[ -1, "var.tot" ]
-    dfr$alpha.stratified[1] <- 1 - sum (( 1 - dfr[ -1, "alpha" ] )*var_tot ) / var_tot
+    dfr_alpha <- dfr[ -1, "alpha" ]
+    # dfr$alpha.stratified[1] <- 1 - sum (( 1 - dfr_alpha )*var_tot ) / var_tot
+    dfr$alpha.stratified[1] <- 1 - sum (( 1 - dfr_alpha )*var_tot ) / dfr[ 1, "var.tot"]
     obji <- dfr
     obji[, -c(1:2)] <- round( obji[,-c(1:2) ], 3 )
     if ( ! stratcomp ){
diff --git a/README.md b/README.md
index cc20cdd..9c7789f 100644
--- a/README.md
+++ b/README.md
@@ -25,9 +25,9 @@ The CRAN version can be installed from within R using:
 utils::install.packages("sirt")
 ```
 
-#### GitHub version `sirt` 3.13-105 (2023-03-19)
+#### GitHub version `sirt` 3.13-128 (2023-04-02)
 
-[![](https://img.shields.io/badge/github%20version-3.13--105-orange.svg)](https://github.com/alexanderrobitzsch/sirt)&#160;&#160;
+[![](https://img.shields.io/badge/github%20version-3.13--128-orange.svg)](https://github.com/alexanderrobitzsch/sirt)&#160;&#160;
 
 The version hosted [here](https://github.com/alexanderrobitzsch/sirt) is the development version of `sirt`. 
 The GitHub version can be installed using `devtools` as:
diff --git a/docs/404.html b/docs/404.html
index 61d6e44..a50f4a8 100644
--- a/docs/404.html
+++ b/docs/404.html
@@ -71,7 +71,7 @@
       </button>
       <span class="navbar-brand">
         <a class="navbar-link" href="index.html">sirt</a>
-        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-105</span>
+        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-128</span>
       </span>
     </div>
 
diff --git a/docs/authors.html b/docs/authors.html
index 312eab0..17ca3b4 100644
--- a/docs/authors.html
+++ b/docs/authors.html
@@ -71,7 +71,7 @@
       </button>
       <span class="navbar-brand">
         <a class="navbar-link" href="index.html">sirt</a>
-        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-105</span>
+        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-128</span>
       </span>
     </div>
 
@@ -111,12 +111,12 @@ <h1>Citation</h1>
       <small class="dont-index">Source: <a href='https://github.com/alexanderrobitzsch/sirt/blob/master/inst/CITATION'><code>inst/CITATION</code></a></small>
     </div>
 
-    <p>Robitzsch, A. (2023). sirt: Supplementary Item Response Theory Models. R package version 3.13-105. https://CRAN.R-project.org/package=sirt</p>
-    <pre>@Manual{sirt_3.13-105,
+    <p>Robitzsch, A. (2023). sirt: Supplementary Item Response Theory Models. R package version 3.13-128. https://CRAN.R-project.org/package=sirt</p>
+    <pre>@Manual{sirt_3.13-128,
   title = {sirt: Supplementary Item Response Theory Models},
   author = {Alexander Robitzsch},
   year = {2023},
-  note = {R package version 3.13-105},
+  note = {R package version 3.13-128},
   url = {https://CRAN.R-project.org/package=sirt},
 }</pre>
 
diff --git a/docs/index.html b/docs/index.html
index a9f299e..2748435 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -47,7 +47,7 @@
       </button>
       <span class="navbar-brand">
         <a class="navbar-link" href="index.html">sirt</a>
-        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-105</span>
+        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">3.13-128</span>
       </span>
     </div>
 
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index e08c236..13245d9 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -2,5 +2,5 @@ pandoc: 1.13.1
 pkgdown: 1.5.1
 pkgdown_sha: ~
 articles: []
-last_built: 2023-03-19T11:58Z
+last_built: 2023-04-02T10:58Z
 
diff --git a/inst/NEWS b/inst/NEWS
index 9b91bc8..22df49d 100644
--- a/inst/NEWS
+++ b/inst/NEWS
@@ -35,7 +35,7 @@ CHANGELOG sirt
 
 
 ------------------------------------------------------------------------
-VERSIONS sirt 3.13  |  2023-03-19  |  Last: sirt 3.13-105
+VERSIONS sirt 3.13  |  2023-04-02  |  Last: sirt 3.13-128
 ------------------------------------------------------------------------
 
 NOTE  * included arguments 'pw_linear', 'pw_quadratic' that allow 
@@ -66,7 +66,7 @@ ADDED * included function lsem.test() for testing of no parameter change
         based on a fitted bootstrap object
 NOTE  * included utility function print_digits() that prints a data frame
         with a fixed number of digits
- 
+FIXED * corrected a bug in stratified.cronbach.alpha()
          
 
 DATA  * included/modified datasets: ---
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 1940306..3dba00d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,5 +1,5 @@
 //// File Name: RcppExports.cpp
-//// File Version: 3.013105
+//// File Version: 3.013128
 // Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393