3.13-1

alexanderrobitzsch · May 17, 2022 · 3d95b68 · 3d95b68
1 parent 5952cba
commit 3d95b68
Show file tree

Hide file tree

Showing 219 changed files with 1,579 additions and 766 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.12-43
-Date: 2022-04-04 11:31:55
+Version: 3.13-1
+Date: 2022-05-17 11:06:20
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 
@@ -30,8 +30,6 @@ Suggests:
     coda, igraph, lavaan, lavaan.survey, MASS, Matrix,
     miceadds, minqa, mirt, mvtnorm, nloptr, optimx, pbivnorm, 
     pbv, psych, sfsmisc, sm, survey
-Enhances:
-    immer
 LinkingTo: 
     pbv, Rcpp, RcppArmadillo
 URL: 

diff --git a/NAMESPACE b/NAMESPACE
@@ -192,6 +192,7 @@ export(btm)
 export(btm_sim)
 export(categorize)
 export(ccov.np)
+export(cfa_meas_inv)
 export(class.accuracy.rasch)
 export(colCumsums.sirt)
 export(conf.detect)
@@ -331,12 +332,14 @@ export(rinvgamma2)
 export(rm.facets)
 export(rm.sdt)
 export(rm_proc_data)
+export(rmvn)
 export(rowCumsums.sirt)
 export(rowIntervalIndex.sirt)
 export(rowKSmallest.sirt)
 export(rowKSmallest2.sirt)
 export(rowMaxs.sirt)
 export(rowMins.sirt)
+export(ruvn)
 export(scale_group_means)
 export(soft_thresholding)
 export(sia.sirt)

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.012043
+## File Version: 3.013001
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
@@ -391,8 +391,8 @@ sirt_rcpp_rm_sdt_calc_probs_grm_item <- function(tau_item, a_item, theta_k, VV,
     .Call('_sirt_sirt_rcpp_rm_sdt_calc_probs_grm_item', PACKAGE='sirt', tau_item, a_item, theta_k, VV, K, TP, eps, use_log)
 }
 
-sirt_rcpp_xxirt_compute_posterior_expected_counts <- function(dat1_resp_gg, p_aj_xi_gg) {
-    .Call('_sirt_sirt_rcpp_xxirt_compute_posterior_expected_counts', PACKAGE='sirt', dat1_resp_gg, p_aj_xi_gg)
+sirt_rcpp_xxirt_compute_posterior_expected_counts <- function(dat1_resp_gg, p_aj_xi_gg, weights_gg) {
+    .Call('_sirt_sirt_rcpp_xxirt_compute_posterior_expected_counts', PACKAGE='sirt', dat1_resp_gg, p_aj_xi_gg, weights_gg)
 }
 
 sirt_rcpp_xxirt_compute_likelihood <- function(dat, dat_resp_bool, probs, TP, maxK) {

diff --git a/R/Rhat_sirt.R b/R/Rhat_sirt.R
@@ -1,5 +1,5 @@
 ## File Name: Rhat_sirt.R
-## File Version: 1.09
+## File Version: 1.12
 
 
 ####################################
@@ -12,7 +12,8 @@
     # Source: Statistical Science, Vol. 7, No. 4 (Nov., 1992), pp. 457-472
     # Stable URL: http://www.jstor.org/stable/2246093
     ## Matches gelman.diag() from package "coda", but not WinBUGS() "summary" component.
-    ## Better than gelman.diag() because multivariate stat is not bothered to be calculated
+    ## Better than gelman.diag() because multivariate stat is not
+    ## bothered to be calculated
 Rhat1 <- function(mat)
 {
     m <- ncol(mat)

diff --git a/R/btm_fit_statistics.R b/R/btm_fit_statistics.R
@@ -1,5 +1,5 @@
 ## File Name: btm_fit_statistics.R
-## File Version: 0.2071
+## File Version: 0.209
 
 
 #**** item outfit and infit statistic
@@ -72,7 +72,8 @@ btm_fit_statistics <- function( probs, dat0, ind1, ind2, TP, judge=NULL,
             a1$N_dyad <- stats::aggregate( 1+0*dat1$result, list(dat1$dyad), sum )[,2]
             a1 <- a1[ ( a1$N_dyad > 2 ) & ( a1$mode %in% c(0,1) ), ]
             dat1$mode <- a1[ match(dat1$dyad, a1$dyad), "mode" ]
-            a2 <- stats::aggregate( dat1$result==dat1$mode, list(dat1$judge), mean, na.rm=TRUE)
+            a2 <- stats::aggregate( dat1$result==dat1$mode, list(dat1$judge),
+                                        mean, na.rm=TRUE)
             fit_judges$agree <- a2[ match(fit_judges$judge, a2[,1]), 2]
         }
     }

diff --git a/R/ccov.np.R b/R/ccov.np.R
@@ -1,5 +1,5 @@
 ## File Name: ccov.np.R
-## File Version: 1.216
+## File Version: 1.218
 
 
 #---- nonparametric estimation of conditional covariance
@@ -52,7 +52,8 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200),
 
     sirt_progress_cat(progress=progress)
     #-- weights thetagrid
-    wgt_thetagrid <- ccov_np_score_density(score=score, thetagrid=thetagrid, smooth=smooth)
+    wgt_thetagrid <- ccov_np_score_density(score=score, thetagrid=thetagrid,
+                            smooth=smooth)
 
     #-- display progress
     if (progress){

diff --git a/R/cfa_meas_inv.R b/R/cfa_meas_inv.R
@@ -0,0 +1,69 @@
+## File Name: cfa_meas_inv.R
+## File Version: 0.157
+
+
+cfa_meas_inv <- function(dat, group, weights=NULL, alpha=0.01, verbose=FALSE,
+                    op=c("~1","=~") )
+{
+    requireNamespace("lavaan")
+    #--- define data input
+    items <- colnames(dat)
+    if (is.null(weights)){
+        weights <- rep(1, nrow(dat))
+    }
+    dat <- data.frame(group=group, weight=weights, dat)
+    args <- list(data=dat, group="group")
+    args$sampling.weights <- "weight"
+
+    #--- define model under scalar invariance
+    lavmodel <- invariance_alignment_cfa_config_estimate_define_lavaan_model(
+                        items_gg=items, label_F="F", model="2PM")
+    args$model <- lavmodel
+    args$meanstructure <- TRUE
+    args$std.lv <- TRUE
+    args$group.equal <- c("loadings","intercepts")
+
+    mod1 <- lavaan::cfa(data=dat, model=lavmodel, std.lv=TRUE, meanstructure=TRUE,
+                    group.equal=c("loadings","intercepts"), group="group",
+                    sampling.weights="weight")
+    partable1 <- lavaan::parameterTable(object=mod1)
+    mimod1 <- mi_inv_lavaan_modification_indices(mod=mod1, op=op)
+    partable0 <- partable1
+    partable <- partable1
+    mod0 <- mod1
+
+    pars_mi <- meas_inv_cfa_proc_partable(partable=partable1, items=items)
+
+    #--- change model by subsequently freeing parameters
+    nfp <- 0
+    critval <- stats::qchisq(1-alpha, df=1)
+    free <- TRUE
+    while(free){
+
+        res <- meas_inv_cfa_modify_partable(partable=partable1, mimod=mimod1,
+                            critval=critval)
+        free <- res$free_parameter
+
+        if (free){
+            nfp <- nfp+1
+            partable2 <- res$partable
+            mod2 <- lavaan::cfa(data=dat, model=partable2, group="group",
+                        sampling.weights="weight")
+            partable1 <- lavaan::parameterTable(object=mod2)
+            mimod1 <- mi_inv_lavaan_modification_indices(mod=mod2, op=op)
+            if (verbose){
+                cat(paste0("freed ", nfp, " parameters\n" ) )
+                utils::flush.console()
+            }
+
+        }
+    }
+
+    pars_pi <- meas_inv_cfa_proc_partable(partable=partable1, items=items)
+
+    #--- output
+    res <- list(pars_mi=pars_mi, pars_pi=pars_pi, alpha=alpha, critval=critval,
+                    nfp=nfp, partable=partable1, dat=dat, items=items,
+                    mod_mi=mod0, mod_pi=mod1)
+    return(res)
+}
diff --git a/R/conf.detect.R b/R/conf.detect.R
@@ -1,5 +1,5 @@
 ## File Name: conf.detect.R
-## File Version: 1.206
+## File Version: 1.208
 
 
 # Confirmatory DETECT analysis
@@ -52,8 +52,9 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE,
                                 "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" )
+        rownames(detect.summary) <- c("DETECT Unweighted", "DETECT Weighted",
+                                        "ASSI Unweighted", "ASSI Weighted",
+                                        "RATIO Unweighted", "RATIO Weighted" )
     }
     cat("-----------------------------------------------------------\n" )
     if ( ! h1){

diff --git a/R/data.prep.R b/R/data.prep.R
@@ -1,5 +1,5 @@
 ## File Name: data.prep.R
-## File Version: 1.145
+## File Version: 1.147
 
 #----- data preparations for rasch.jml and rasch.mml
 data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE,
@@ -43,10 +43,11 @@ data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE,
     }
     # item pattern corresponding to frequency pattern
     if ( use.freqpatt){
-        dat2 <- matrix( as.numeric( unlist( strsplit( paste(dat1[,1]), "" ) ) ), ncol=ncol(dat), byrow=T)
+        dat2 <- matrix( as.numeric( unlist( strsplit( paste(dat1[,1]), "" ) ) ),
+                            ncol=ncol(dat), byrow=TRUE)
     } else {
         dat2 <- dat.9 }
-        dat2.resp <- 1 * ( dat2 !=9 )
+        dat2.resp <- 1*( dat2 !=9 )
         dat2[ dat2==9 ] <- 0
         # mean right
         dat1$mean <- rowSums( dat2 * dat2.resp )  / rowSums( dat2.resp )

diff --git a/R/equating.rasch.R b/R/equating.rasch.R
@@ -1,5 +1,5 @@
 ## File Name: equating.rasch.R
-## File Version: 0.242
+## File Version: 0.244
 
 
 #---- Equating (linking) in the Rasch model
@@ -18,15 +18,19 @@ equating.rasch <- function( x, y, theta=seq( -4, 4, len=100),
     opt_interval <- 10*c(-1,1)
     #-- Haebara function
     ha <- function(B){
-                fct1 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,2], alpha1=alpha1, alpha2=alpha2 )
-                fct2 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,3] - B, alpha1=alpha1, alpha2=alpha2 )
+                fct1 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,2], alpha1=alpha1,
+                                    alpha2=alpha2 )
+                fct2 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,3] - B,
+                                    alpha1=alpha1, alpha2=alpha2 )
                 sum( (fct1 - fct2)^2 )
             }
     B.ha <- stats::optimize( f=ha, interval=opt_interval )$minimum
     # Stocking and Lord Approach
     sl <- function(B){
-                fct1 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,2], alpha1=alpha1, alpha2=alpha2 )
-                fct2 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,3] - B, alpha1=alpha1, alpha2=alpha2 )
+                fct1 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,2],
+                                alpha1=alpha1, alpha2=alpha2 )
+                fct2 <- .prob.raschtype.genlogis( theta=theta, b=b.xy[,3] - B,
+                                alpha1=alpha1, alpha2=alpha2 )
                 sum( (rowSums( fct1 - fct2 ) )^2 )
             }
     B.sl <- stats::optimize( f=sl, interval=opt_interval )$minimum
@@ -41,7 +45,8 @@ equating.rasch <- function( x, y, theta=seq( -4, 4, len=100),
     colnames(transf.par) <- c("item", "TransfItempar.Gr1", "Itempar.Gr2"  )
     transf.par <- transf.par[ order( paste(transf.par$item ) ), ]
     # calculate variance and linking error
-    des <- data.frame( "N.Items"=nrow(b.xy), "SD"=stats::sd( b.xy$TransfItempar.Gr1 - b.xy$Itempar.Gr2 ) )
+    des <- data.frame( N.Items=nrow(b.xy),
+                            SD=stats::sd( b.xy$TransfItempar.Gr1 - b.xy$Itempar.Gr2 ) )
     des$Var <- des$SD^2
     des$linkerror <- sqrt( des["SD"]^2 / des["N.Items"] )[1,1]
     #--- output

diff --git a/R/gom.em.R b/R/gom.em.R
@@ -1,5 +1,5 @@
 ## File Name: gom.em.R
-## File Version: 5.411
+## File Version: 5.420
 
 
 #-- gom EM algorithm
@@ -197,7 +197,8 @@ gom.em <- function( dat, K=NULL, problevels=NULL, weights=NULL, model="GOM",
     #--- Newton-Raphson steps
     if (newton_raphson){
         lambda_logit <- as.vector( stats::qlogis(lambda) )
-        lambda_logit <- stats::aggregate(lambda_logit, list(lambda_partable$par_index), mean)[,2]
+        lambda_logit <- stats::aggregate(lambda_logit,
+                                list(lambda_partable$par_index), mean)[,2]
         ind_lambda <- seq(1, max(lambda_partable$par_index))
         names(lambda_logit) <- lambda_partable$par_name[ lambda_partable$free ]
         pi_k_logit <- sirt_probs_to_logit(y=pi.k)
@@ -299,7 +300,6 @@ gom.em <- function( dat, K=NULL, problevels=NULL, weights=NULL, model="GOM",
         mu <- res$mu
         Sigma <- res$Sigma
         pi.k <- res$pi.k
-
         #--- distribution parameters
         c1 <- stats::cov2cor(Sigma)
         if (progress){
@@ -384,14 +384,14 @@ gom.em <- function( dat, K=NULL, problevels=NULL, weights=NULL, model="GOM",
     #--- output
     s2 <- Sys.time()
     res <- list(deviance=dev, ic=ic, item=item, person=person, EAP.rel=EAP.rel,
-        MAP=MAP, EAP=EAP, classdesc=classdesc, lambda=lambda, se.lambda=se.lambda,
-        mu=mu, Sigma=Sigma, b=b, se.b=se.b, f.yi.qk=f.yi.qk, f.qk.yi=f.qk.yi,
-        probs=probs, n.ik=n.ik, iter=iter, dat=dat0, dat2=dat2, dat2.resp=dat2.resp,
-        I=I, K=K, TP=TP, G=1, theta.k=theta.k, pi.k=pi.k, problevels=problevels,
-        model=model, plmat=plmat, mu=mu, Sigma=Sigma, weights=weights,
-        lambda.index=lambda.index, lambda_partable=lambda_partable,
-        optimization=optimization, newton_raphson=newton_raphson,
-        s1=s1, s2=s2, time_diff=s2-s1, CALL=CALL)
+                MAP=MAP, EAP=EAP, classdesc=classdesc, lambda=lambda, se.lambda=se.lambda,
+                mu=mu, Sigma=Sigma, b=b, se.b=se.b, f.yi.qk=f.yi.qk, f.qk.yi=f.qk.yi,
+                probs=probs, n.ik=n.ik, iter=iter, dat=dat0, dat2=dat2, dat2.resp=dat2.resp,
+                I=I, K=K, TP=TP, G=1, theta.k=theta.k, pi.k=pi.k, problevels=problevels,
+                model=model, plmat=plmat, mu=mu, Sigma=Sigma, weights=weights,
+                lambda.index=lambda.index, lambda_partable=lambda_partable,
+                optimization=optimization, newton_raphson=newton_raphson,
+                s1=s1, s2=s2, time_diff=s2-s1, CALL=CALL)
     class(res) <- "gom"
     return(res)
 }

diff --git a/R/gom_em_loglike_opt_fun.R b/R/gom_em_loglike_opt_fun.R
@@ -1,5 +1,5 @@
 ## File Name: gom_em_loglike_opt_fun.R
-## File Version: 0.267
+## File Version: 0.271
 
 
 
@@ -9,8 +9,8 @@ gom_em_loglike_opt_fun <- function(x, ind_lambda, ind_pi, I, K, theta.k, theta0.
 {
     ncat <- 2
     res <- gom_em_loglike_parameter_conversion( x=x, ind_lambda=ind_lambda,
-        ind_pi=ind_pi, I=I, K=K, ind_mu=ind_mu, ind_sigma=ind_sigma, model=model,
-        theta_grid=theta_grid, lambda_partable=lambda_partable )
+                ind_pi=ind_pi, I=I, K=K, ind_mu=ind_mu, ind_sigma=ind_sigma, model=model,
+                theta_grid=theta_grid, lambda_partable=lambda_partable )
     lambda <- res$lambda
     pi.k <- res$pi.k
     mu <- res$mu

diff --git a/R/gom_em_loglike_parameter_conversion.R b/R/gom_em_loglike_parameter_conversion.R
@@ -1,9 +1,9 @@
 ## File Name: gom_em_loglike_parameter_conversion.R
-## File Version: 0.18
+## File Version: 0.192
 
 
 gom_em_loglike_parameter_conversion <- function(x, ind_lambda, ind_pi, I, K,
-    ind_mu, ind_sigma, model, theta_grid, lambda_partable)
+            ind_mu, ind_sigma, model, theta_grid, lambda_partable)
 {
     mu <- NULL
     Sigma <- NULL