4.2-1

alexanderrobitzsch · Feb 6, 2024 · a4bff44 · a4bff44
1 parent cfec9d6
commit a4bff44
Show file tree

Hide file tree

Showing 33 changed files with 239 additions and 156 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 4.1-1
-Date: 2024-01-18 10:36:46
+Version: 4.2-1
+Date: 2024-02-06 09:12:35
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 4.001001
+## File Version: 4.002001
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 

diff --git a/R/conf.detect.R b/R/conf.detect.R
@@ -1,5 +1,5 @@
 ## File Name: conf.detect.R
-## File Version: 1.208
+## File Version: 1.211
 
 
 # Confirmatory DETECT analysis
@@ -8,21 +8,21 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE,
                 bias_corr=TRUE)
 {
     CALL <- match.call()
-    cat("-----------------------------------------------------------\n" )
-    cat("Confirmatory DETECT Analysis \n" )
+    cat('-----------------------------------------------------------\n' )
+    cat('Confirmatory DETECT Analysis \n' )
     utils::flush.console()
     h1 <- is.matrix(score)
     if (h1){
         PP <- ncol(score)
     }
     is_one_score <- TRUE
     if (! h1 ){
-        cat("Conditioning on 1 Score\n" )
+        cat('Conditioning on 1 Score\n' )
     } else {
-        cat(paste("Conditioning on ",PP, " Scores\n", sep="") )
+        cat(paste('Conditioning on ',PP, ' Scores\n', sep='') )
         is_one_score <- FALSE
     }
-    cat(paste("Bandwidth Scale:", bwscale, "\n" ) )
+    cat(paste('Bandwidth Scale:', bwscale, '\n' ) )
     utils::flush.console()
     scale_score <- TRUE
     if (!smooth){
@@ -39,7 +39,7 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE,
         ccovtable.list <- list()
         args_ccov_np$progress <- FALSE
         for (pp in 1:PP){
-            cat( paste( "DETECT Calculation Score ", pp, "\n", sep="") ) ;
+            cat( paste( 'DETECT Calculation Score ', pp, '\n', sep='') ) ;
             utils::flush.console()
             args_ccov_np$score <- score[,pp]
             ccovtable.list[[pp]] <- do.call( what=ccov.np, args=args_ccov_np)
@@ -48,21 +48,20 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE,
                     detect.index( ccovtable, itemcluster=itemcluster ) } )
         detect.matrix <- matrix( unlist( lapply( detect.list, FUN=function( ll){
                                 c( ll[1,], ll[2,], ll[3,] ) } ) ), nrow=PP, byrow=TRUE)
-        detect.summary <- data.frame( "NScores"=PP, "Mean"=colMeans( detect.matrix ),
-                                "SD"=apply( detect.matrix, 2, stats::sd ),
-                                "Min"=apply( detect.matrix, 2, min ),
-                                "Max"=apply( detect.matrix, 2, max ) )
-        rownames(detect.summary) <- c("DETECT Unweighted", "DETECT Weighted",
-                                        "ASSI Unweighted", "ASSI Weighted",
-                                        "RATIO Unweighted", "RATIO Weighted" )
+        detect.summary <- data.frame( NScores=PP, Mean=colMeans( detect.matrix ),
+                                SD=apply( detect.matrix, 2, stats::sd ),
+                                Min=apply( detect.matrix, 2, min ),
+                                Max=apply( detect.matrix, 2, max ) )
+        rownames(detect.summary) <- c('DETECT Unweighted', 'DETECT Weighted',
+                                        'ASSI Unweighted', 'ASSI Weighted',
+                                        'RATIO Unweighted', 'RATIO Weighted' )
     }
-    cat("-----------------------------------------------------------\n" )
+    cat('-----------------------------------------------------------\n' )
     if ( ! h1){
-        res <- list( "detect"=res, "ccovtable"=ccovtable,
-                    "detect.summary"=res )
+        res <- list( detect=res, ccovtable=ccovtable, detect.summary=res )
     } else {
-        res <- list( "detect"=detect.list, "ccovtable"=ccovtable.list,
-            "detect.summary"=detect.summary )
+        res <- list( detect=detect.list, ccovtable=ccovtable.list,
+                        detect.summary=detect.summary )
     }
     res$is_one_score <- is_one_score
     res$CALL <- CALL
@@ -72,7 +71,7 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE,
     #--- print
     print(round(res$detect.summary,3))
     #--- return
-    class(res) <- "conf.detect"
+    class(res) <- 'conf.detect'
     return(res)
 }
 
diff --git a/R/invariance.alignment.R b/R/invariance.alignment.R
@@ -1,5 +1,5 @@
 ## File Name: invariance.alignment.R
-## File Version: 3.820
+## File Version: 3.962
 
 
 invariance.alignment <- function( lambda, nu, wgt=NULL,
@@ -26,9 +26,22 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
         wgt <- 1+0*nu
     }
 
+    #- reparametrization of meth argument
+    meth0 <- meth
+    if (meth0==4){ meth <- 0}     # Mplus FREE
+    if (meth0==3){ meth <- 0.5}      # Mplus FIXED
+
     #- choose fixed value
     fixed <- invariance_alignment_choose_fixed(fixed=fixed, G=G, Gmax=999)
     reparam <- ! fixed
+    if (meth%in%c(0,0.5)){
+        constraint <- 'prod'
+        reparam <- TRUE
+        num_deriv <- TRUE
+    }
+    if (overparam){
+        num_deriv <- TRUE
+    }
 
     W1 <- dim(wgt)
     wgtM <- matrix( colSums(wgt,na.rm=TRUE), nrow=W1[1], ncol=W1[2], byrow=TRUE )
@@ -66,34 +79,45 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
 
     group_combis <- group.combis-1
     G1 <- G-1
+    if (overparam){
+        G1 <- G
+    }
     ind_alpha <- seq_len(G1)
     ind_psi <- G1 + ind_alpha
+    if (meth==0){
+        ind_alpha <- seq_len(G)
+        ind_psi <- G + seq_len(G1)
+    }
 
     #-- define optimization functions
-    ia_fct_optim <- function(x, lambda, nu, overparam, eps){
+    ia_fct_optim <- function(x, lambda, nu, overparam, eps, meth_ ){
         res <- invariance_alignment_define_parameters(x=x, ind_alpha=ind_alpha,
-                        ind_psi=ind_psi, reparam=reparam)
+                        ind_psi=ind_psi, reparam=reparam, meth=meth_)
+        alpha0 <- res$alpha0
+        psi0 <- res$psi0
         val <- sirt_rcpp_invariance_alignment_opt_fct( nu=nu, lambda=lambda,
-                        alpha0=res$alpha0, psi0=res$psi0, group_combis=group_combis,
+                        alpha0=alpha0, psi0=psi0, group_combis=group_combis,
                         wgt=wgt, align_scale=align.scale,
                         align_pow=align.pow, eps=eps, wgt_combi=wgt_combi, type=type,
-                        reparam=FALSE, meth=meth)
+                        reparam=FALSE, meth=meth_)
         val <- val$fopt
-        if (overparam){
+        if (overparam | meth==0 ){
             G <- nrow(lambda)
-            val <- val+sum(wgt[,1]*x[1:G]^2 + wgt[,1]*x[G+(1:G)]^2 )
+            fac <- sum(wgt[,1]) / 1000
+            val <- val+fac*sum(x[1:G]^2)
         }
         return(val)
     }
-    ia_grad_optim <- function(x, lambda, nu, overparam, eps){
+
+    ia_grad_optim <- function(x, lambda, nu, overparam, eps, meth_){
         res <- invariance_alignment_define_parameters(x=x, ind_alpha=ind_alpha,
                             ind_psi=ind_psi, reparam=reparam)
         alpha0 <- res$alpha0
         psi0 <- res$psi0
         grad <- sirt_rcpp_invariance_alignment_opt_grad( nu=nu, lambda=lambda,
                         alpha0=alpha0, psi0=psi0, group_combis=group_combis, wgt=wgt,
                         align_scale=align.scale, align_pow=align.pow, eps=eps,
-                        wgt_combi=wgt_combi, type=type, reparam=reparam, meth=meth)
+                        wgt_combi=wgt_combi, type=type, reparam=reparam, meth=meth_)
         grad <- grad[-c(1,G+1)]
         return(grad)
     }
@@ -106,8 +130,14 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
 
     #** evaluate optimization function at initial solution
     x0 <- c( alpha0[-1], psi0[-1] )
+    if (overparam){
+        x0 <- c(alpha0, psi0)
+    }
+    if (meth==0){
+        x0 <- c( alpha0, psi0[-1] )
+    }
     fct_optim_inits <- ia_fct_optim(x=x0, lambda=lambda, nu=nu,
-                            overparam=overparam, eps=eps)
+                            overparam=overparam, eps=eps, meth=meth )
 
     #* estimate alignment parameters
     min_val <- .01
@@ -117,11 +147,16 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
         GL <- G
         par <- c( alpha0, psi0 )
     }
+    if (meth==0){
+        par <- c( alpha0, psi0[-1] )
+    }
     lower <- c(rep(-Inf,GL), rep(min_val, GL))
     if (reparam){
         grad_optim <- NULL
     }
-
+    if (meth==0){
+        lower <- c(rep(-Inf,G), rep(min_val, GL))
+    }
     #* define sequence of epsilon values
     eps_vec <- 10^eps_grid
     eps_vec <- sirt_define_eps_sequence(eps=eps, eps_vec=eps_vec)
@@ -131,19 +166,21 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
     nu1 <- nu
     est_loop <- 1
     psi_list <- list()
+
     while(est_loop>=1){
         for (eps in eps_vec){
             res_optim <- sirt_optimizer(optimizer=optimizer, par=par, fn=ia_fct_optim,
                                 grad=ia_grad_optim, lower=lower, hessian=FALSE,
                                 lambda=lambda1, nu=nu1, overparam=overparam,
-                                eps=eps, ...)
+                                eps=eps, meth_=meth, ...)
             par <- res_optim$par
             res <- invariance_alignment_define_parameters(x=res_optim$par,
-                            ind_alpha=ind_alpha, ind_psi=ind_psi, reparam=reparam)
+                            ind_alpha=ind_alpha, ind_psi=ind_psi, reparam=reparam,
+                            meth=meth)
             alpha0 <- res$alpha0
             psi0 <- res$psi0
-        }
-        if (meth%in%c(1,2)){
+        } # end eps grid loop
+        if (meth%in%c(0,0.5,1,2)){
             est_loop <- 0
         }
         if (meth>=3){
@@ -153,10 +190,12 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
             }
             if (est_loop==1){
                 est_loop <- est_loop + 1
-                nu1 <- nu/lambda*psi_list[[1]]
+                # nu1 <- nu/sqrt(mod1)*psi_list[[1]]
+                #=> likely dead code
             }
         }
-    }
+    }  # end while loop
+
     if (meth==3){
         psi0 <- psi_list[[1]]
     }
@@ -168,7 +207,7 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
         NP <- length(par)
 
         #- gradient computation
-        ia_grad_optim_num <- function(x, lambda, nu, overparam, eps, h=1e-4){
+        ia_grad_optim_num <- function(x, lambda, nu, overparam, eps, meth=1, h=1e-4){
             NP <- length(x)
             par <- x
             grad <- rep(0,NP)
@@ -188,7 +227,8 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
         h <- 1e-4
         hess_par <- matrix(NA, nrow=NP, ncol=NP)
         rownames(hess_par) <- colnames(hess_par) <- names(par)
-        args <- list(x=par, lambda=lambda, nu=nu, overparam=overparam, eps=eps)
+        args <- list(x=par, lambda=lambda, nu=nu, overparam=overparam, eps=eps,
+                        meth=meth)
         pp <- 1
         for (pp in 1:NP){
             args$x <- mgsem_add_increment(x=par, h=h, i1=pp)
@@ -230,10 +270,9 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
 
     }
 
-
     # center parameters
     res <- invariance_alignment_center_parameters(alpha0=alpha0, psi0=psi0,
-                    center=center)
+                    center=center, meth=meth )
     alpha0 <- res$alpha0
     psi0 <- res$psi0
 
@@ -286,8 +325,8 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
             nu=nu0, nu.resid=nu.resid, fopt=fopt, align.scale=align.scale,
             align.pow=align.pow0, res_optim=res_optim, eps=eps, wgt=wgt,
             miss_items=missM, numb_items=numb_items, vcov=vcov,
-            fct_optim_inits=fct_optim_inits, fixed=fixed, meth=meth,
-            s1=s1, s2=s2, time_diff=time_diff, CALL=CALL)
+            fct_optim_inits=fct_optim_inits, fixed=fixed, meth=meth0,
+            meth_internal=meth, s1=s1, s2=s2, time_diff=time_diff, CALL=CALL)
     class(res) <- 'invariance.alignment'
     return(res)
 }

diff --git a/R/invariance_alignment_center_parameters.R b/R/invariance_alignment_center_parameters.R
@@ -1,8 +1,8 @@
 ## File Name: invariance_alignment_center_parameters.R
-## File Version: 0.10
+## File Version: 0.130
 
 invariance_alignment_center_parameters <- function(alpha0, psi0, center,
-    reparam=FALSE, convert=FALSE)
+    reparam=FALSE, convert=FALSE, meth=1)
 {
     if (reparam & convert){
         alpha0 <- alpha0 * psi0
@@ -12,6 +12,14 @@ invariance_alignment_center_parameters <- function(alpha0, psi0, center,
         log_psi <- log(psi0)
         psi0 <- exp(log_psi - mean(log_psi))
     }
+    if (meth %in% c(0,0.5) ){
+        fac <- 1/psi0[1]
+        psi0 <- fac*psi0
+        alpha0 <- fac*alpha0
+        if (meth==0){
+            alpha0 <- alpha0 - alpha0[1]
+        }
+    }
     #--- output
     res <- list(alpha0=alpha0, psi0=psi0)
     return(res)

diff --git a/R/invariance_alignment_define_parameters.R b/R/invariance_alignment_define_parameters.R
@@ -1,23 +1,35 @@
 ## File Name: invariance_alignment_define_parameters.R
-## File Version: 0.193
+## File Version: 0.215
 
 invariance_alignment_define_parameters <- function(x, ind_alpha, ind_psi,
-        fix_first_psi=TRUE, reparam=FALSE, constraint="prod", overparam=FALSE)
+        fix_first_psi=TRUE, reparam=FALSE, constraint="prod", overparam=FALSE,
+        meth=1)
 {
+    # meth=0=> free optimization with ridge regularization
     alpha0 <- c(0, x[ind_alpha])
+    if (meth==0){
+        alpha0 <- alpha0[-c(1)]
+    }
     if (fix_first_psi){
         psi0 <- c(1, x[ind_psi])
     } else {
         psi0 <- c(x[ind_psi])
     }
-    if (reparam & ( constraint=='prod') ){
+    if (( reparam & ( constraint=='prod') ) | (meth %in% c(0,0.5) ) ){
         prod_psi <- prod(psi0)
         NX <- length(psi0)
-        psi0 <- psi0 / ( prod_psi^(1/NX) )
+        log_psi0 <- log(psi0)
+        psi0 <- exp( log_psi0 - mean(log_psi0) )
+        # psi0 <- psi0 / ( prod_psi^(1/NX) )
     }
     if (reparam & ( constraint=='sum') ){
         psi0 <- psi0 - mean(psi0) + 1
     }
+    if (meth==0){
+        # alpha0[1] <- -sum(alpha0[-c(1)])
+        # alpha0 <- alpha0 - mean(alpha0)
+        # Implementing the constraints leads to bias
+    }
     #--- output
     res <- list(alpha0=alpha0, psi0=psi0)
     return(res)