3.13-128

alexanderrobitzsch · Apr 2, 2023 · 0e25144 · 0e25144
1 parent 9dd6b9f
commit 0e25144
Show file tree

Hide file tree

Showing 53 changed files with 508 additions and 219 deletions.
diff --git 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 <[email protected]>
 Description: 

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -1,31 +1,31 @@
 ## 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,
     pen_type="lasso")
 {
 
     # 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)
         }
     }