From 0e25144d927cee586537e31f89392764d062690a Mon Sep 17 00:00:00 2001
From: Robitzsch <robitzsch@ipn.uni-kiel.de>
Date: Sun, 2 Apr 2023 14:59:21 +0200
Subject: [PATCH] 3.13-128
---
DESCRIPTION | 4 +-
R/L0_polish.R | 8 +-
R/RcppExports.R | 2 +-
R/dexppow.R | 6 +-
R/lq_fit.R | 4 +-
R/lq_fit_estimate_power.R | 4 +-
R/lsdm_est_logist_quant.R | 8 +-
R/mgsem.R | 50 ++++++++----
R/mgsem_bdiag.R | 14 ++++
R/mgsem_compute_model_implied_moments.R | 5 +-
R/mgsem_duplication_matrix.R | 26 +++++++
R/mgsem_eval_lp_penalty_matrix.R | 8 +-
R/mgsem_eval_lp_penalty_vector.R | 14 ++--
R/mgsem_evaluate_penalties.R | 44 ++++++++++-
R/mgsem_evaluate_penalties_evaluate_entry.R | 12 +--
R/mgsem_ginv.R | 10 +++
R/mgsem_grad_fun.R | 4 +-
R/mgsem_loglike_data.R | 4 +-
R/mgsem_loglike_suffstat.R | 5 +-
R/mgsem_loglike_suffstat_derivative.R | 5 +-
R/mgsem_moments_derivative_parameter.R | 30 ++++----
R/mgsem_observed_information.R | 19 +++--
R/mgsem_opt_fun.R | 6 +-
R/mgsem_output_proc_casewise_likelihood.R | 6 +-
R/mgsem_partable2coef.R | 6 +-
R/mgsem_partable2model.R | 16 ++--
R/mgsem_power_fun_differentiable_approx.R | 10 +--
R/mgsem_proc_data.R | 9 ++-
R/mgsem_proc_model.R | 55 ++++++++-----
R/mgsem_proc_model_add_specs_all.R | 30 ++++----
R/mgsem_proc_model_difflp_information.R | 16 ++--
R/mgsem_proc_model_include_missing_entries.R | 14 ++--
R/mgsem_proc_model_is_B.R | 4 +-
R/mgsem_proc_model_partable_define_index.R | 20 ++---
R/mgsem_proc_model_single_group.R | 4 +-
R/mgsem_proc_model_update_penalties_matrix.R | 4 +-
R/mgsem_proc_suffstat.R | 14 ++--
R/mgsem_proc_technical.R | 26 +++----
R/mgsem_sqrt_diag.R | 7 ++
R/mgsem_suffstat_covariance_matrix.R | 60 +++++++++++++++
R/mgsem_test_fun.R | 8 +-
R/mgsem_vcov_me.R | 81 ++++++++++++++++++++
R/mgsem_vec.R | 7 ++
R/mgsem_vech.R | 4 +-
R/mi_inv_lavaan_modification_indices.R | 4 +-
R/stratified.cronbach.alpha.R | 6 +-
README.md | 4 +-
docs/404.html | 2 +-
docs/authors.html | 8 +-
docs/index.html | 2 +-
docs/pkgdown.yml | 2 +-
inst/NEWS | 4 +-
src/RcppExports.cpp | 2 +-
53 files changed, 508 insertions(+), 219 deletions(-)
create mode 100644 R/mgsem_bdiag.R
create mode 100644 R/mgsem_duplication_matrix.R
create mode 100644 R/mgsem_ginv.R
create mode 100644 R/mgsem_sqrt_diag.R
create mode 100644 R/mgsem_suffstat_covariance_matrix.R
create mode 100644 R/mgsem_vcov_me.R
create mode 100644 R/mgsem_vec.R
diff --git a/DESCRIPTION b/DESCRIPTION
index bd8bd435..8b86caef 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 4b1bf372..58171b45 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 8a03ac4c..691cb05f 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 c2500c4a..2ea66f1e 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 883cbdd8..58a981f8 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 4aa1500d..95141b45 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 e1385684..b769fe18 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 e48e4f92..673cd0ac 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 00000000..a5621a89
--- /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 5711cd8c..a94beb1c 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 00000000..4180bbfa
--- /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 e76d0e67..3ca424c2 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 7a22bdc4..354e553e 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 1e396040..da8e5fdf 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 fdff3cc8..9bf37449 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 00000000..4c04fe97
--- /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 753b035a..6f97e86c 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 1a6779f1..fc42de87 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 f186d343..6eb54685 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 db7b18b3..30c6946f 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 9fe839af..55ef1ac4 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 ef27cd17..28e183fe 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 835a4f57..6c1854fd 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 544ae027..2e2ea77b 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 54e19517..f833672a 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 fca60b8b..fc37656d 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 b17b29a1..41c5bb19 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 13680fd6..c80914c3 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 caf80fb9..66980e14 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 64f90a46..9bf4a69a 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 76cec9a1..0478b476 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 3759a32d..9037638f 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 ef8ca46a..7de57235 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 7bbd1ba4..d9f53500 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 3b63a322..511ca426 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 4dd1e6c1..cc31e14c 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 fd83b106..bd1163ff 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 b08d3933..8bb16266 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 00000000..e29cdb5e
--- /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 00000000..b2bef27b
--- /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 86822d86..4399937d 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 00000000..f9609c0a
--- /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 00000000..acba4e05
--- /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 a482c168..eafea972 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 d6509efa..7ae80a50 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 dac0cd0a..f0b845c2 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 cc20cdd6..9c7789ff 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://github.com/alexanderrobitzsch/sirt)  
+[](https://github.com/alexanderrobitzsch/sirt)  
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 61d6e445..a50f4a8f 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 312eab0c..17ca3b44 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 a9f299e9..2748435f 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 e08c2367..13245d9e 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 9b91bc88..22df49dc 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 1940306a..3dba00dc 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