4.2-51

alexanderrobitzsch · Apr 12, 2024 · 8e16305 · 8e16305
1 parent a9dea4f
commit 8e16305
Show file tree

Hide file tree

Showing 103 changed files with 1,468 additions and 351 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 4.2-40
-Date: 2024-03-19 13:56:57
+Version: 4.2-51
+Date: 2024-04-12 16:00:09
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

diff --git a/NAMESPACE b/NAMESPACE
@@ -398,6 +398,7 @@ export(xxirt_createParTable)
 export(xxirt_createThetaDistribution)
 export(xxirt_hessian)
 export(xxirt_modifyParTable)
+export(xxirt_sandwich_pml)
 export(yen.q3)		# defunct
 
 

diff --git a/R/IRT.anova.sirt.R b/R/IRT.anova.sirt.R
@@ -0,0 +1,41 @@
+## File Name: IRT.anova.sirt.R
+## File Version: 0.02
+
+
+
+IRT.anova.sirt <- function (object, ...)
+{
+    cl <- match.call()
+    cl2 <- paste(cl)[-1]
+    if (length(list(object, ...)) !=2) {
+        stop('anova method can only be applied for comparison of two models.\n')
+    }
+    objects <- list(object, ...)
+    model1a <- objects[[1]]
+    model2a <- objects[[2]]
+    model1 <- IRT.IC(model1a)
+    model2 <- IRT.IC(model2a)
+    dfr1 <- data.frame(Model=cl2[1], loglike=model1['loglike'],
+                Deviance=-2 * model1['loglike'])
+    dfr1$Npars <- model1['Npars']
+    dfr1$AIC <- model1['AIC']
+    dfr1$BIC <- model1['BIC']
+    dfr2 <- data.frame(Model=cl2[2], loglike=model2['loglike'],
+                Deviance=-2 * model2['loglike'])
+    dfr2$Npars <- model2['Npars']
+    dfr2$AIC <- model2['AIC']
+    dfr2$BIC <- model2['BIC']
+    dfr <- rbind(dfr1, dfr2)
+    dfr <- dfr[ order(dfr$Npars), ]
+    dfr$Chisq <- NA
+    dfr$df <- NA
+    dfr$p <- NA
+    dfr[1, 'Chisq'] <- dfr[1, 'Deviance'] - dfr[2, 'Deviance']
+    dfr[1, 'df'] <- abs(dfr[1, 'Npars'] - dfr[2, 'Npars'])
+    dfr[1, 'p'] <- round(1 - stats::pchisq(dfr[1, 'Chisq'], df=dfr[1,'df']), 5)
+    for (vv in 2L:(ncol(dfr))){
+        dfr[, vv] <- round(dfr[, vv], 5)
+    }
+    rownames(dfr) <- NULL
+    invisible(dfr)
+}
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 4.002040
+## File Version: 4.002051
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
@@ -407,3 +407,15 @@ sirt_rcpp_xxirt_newton_raphson_derivative_par <- function(dat, dat_resp_bool, ra
     .Call('_sirt_sirt_rcpp_xxirt_newton_raphson_derivative_par', PACKAGE='sirt', dat, dat_resp_bool, ratio, p_xi_aj, item, prior_Theta, group0, weights, ll_case0, eps)
 }
 
+sirt_rcpp_xxirt_nr_pml_opt_fun <- function(prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps) {
+    .Call('_sirt_sirt_rcpp_xxirt_nr_pml_opt_fun', PACKAGE='sirt', prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps)
+}
+
+sirt_rcpp_xxirt_nr_pml_grad_fun_eval <- function(prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps, NP, der_prior_Theta, val1, val2, pp_Theta, der_probs_items, index_freq1, index_freq2) {
+    .Call('_sirt_sirt_rcpp_xxirt_nr_pml_grad_fun_eval', PACKAGE='sirt', prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps, NP, der_prior_Theta, val1, val2, pp_Theta, der_probs_items, index_freq1, index_freq2)
+}
+
+sirt_rcpp_xxirt_nr_pml_casewise_opt_fun <- function(prior_Theta, probs_items, W1, W2_long, G, K, I, TP, NI2, eps, group0, weights, dat1, dat_resp) {
+    .Call('_sirt_sirt_rcpp_xxirt_nr_pml_casewise_opt_fun', PACKAGE='sirt', prior_Theta, probs_items, W1, W2_long, G, K, I, TP, NI2, eps, group0, weights, dat1, dat_resp)
+}
+
diff --git a/R/anova_sirt.R b/R/anova_sirt.R
@@ -1,5 +1,5 @@
 ## File Name: anova_sirt.R
-## File Version: 0.242
+## File Version: 0.244
 
 ##############################################################
 # anova rasch.mml
@@ -39,44 +39,6 @@ anova.xxirt <- anova.rasch.mml
 ##############################################################
 
 
-##################################################################
-IRT.anova.sirt <- function (object, ...)
-{
-    cl <- match.call()
-    cl2 <- paste(cl)[-1]
-    if (length(list(object, ...)) !=2) {
-        stop('anova method can only be applied for comparison of two models.\n')
-    }
-    objects <- list(object, ...)
-    model1a <- objects[[1]]
-    model2a <- objects[[2]]
-    model1 <- IRT.IC(model1a)
-    model2 <- IRT.IC(model2a)
-    dfr1 <- data.frame(Model=cl2[1], loglike=model1['loglike'],
-                Deviance=-2 * model1['loglike'])
-    dfr1$Npars <- model1['Npars']
-    dfr1$AIC <- model1['AIC']
-    dfr1$BIC <- model1['BIC']
-    dfr2 <- data.frame(Model=cl2[2], loglike=model2['loglike'],
-                Deviance=-2 * model2['loglike'])
-    dfr2$Npars <- model2['Npars']
-    dfr2$AIC <- model2['AIC']
-    dfr2$BIC <- model2['BIC']
-    dfr <- rbind(dfr1, dfr2)
-    dfr <- dfr[ order(dfr$Npars), ]
-    dfr$Chisq <- NA
-    dfr$df <- NA
-    dfr$p <- NA
-    dfr[1, 'Chisq'] <- dfr[1, 'Deviance'] - dfr[2, 'Deviance']
-    dfr[1, 'df'] <- abs(dfr[1, 'Npars'] - dfr[2, 'Npars'])
-    dfr[1, 'p'] <- round(1 - stats::pchisq(dfr[1, 'Chisq'], df=dfr[1,'df']), 5)
-    for (vv in 2:(ncol(dfr))){
-        dfr[, vv] <- round(dfr[, vv], 5)
-    }
-    rownames(dfr) <- NULL
-    invisible(dfr)
-}
-
 ##############################################################
 # Likelihood ratio test for rasch.copula2 objects
 anova.rasch.copula2 <- function( object, ... )
@@ -116,7 +78,7 @@ anova.rasch.copula2 <- function( object, ... )
     dfr[1,'Chisq'] <- dfr[1,'Deviance'] - dfr[2,'Deviance']
     dfr[1,'df'] <- abs( dfr[1,'Npars'] - dfr[2,'Npars'] )
     dfr[1, 'p' ] <- round( 1 - stats::pchisq( dfr[1,'Chisq'], df=dfr[1,'df'] ), 5 )
-    for ( vv in 2:( ncol(dfr))){
+    for ( vv in 2L:( ncol(dfr))){
         dfr[,vv] <- round( dfr[,vv], 5 )
     }
     print(dfr)

diff --git a/R/btm.R b/R/btm.R
@@ -1,5 +1,5 @@
 ## File Name: btm.R
-## File Version: 1.537
+## File Version: 1.539
 
 
 #--- Bradley-Terry model in sirt
@@ -20,7 +20,7 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU
     }
 
     if ( ignore.ties ){
-        admiss <- admiss[1:2]
+        admiss <- admiss[c(1,2)]
         delta <- -99
         est.delta <- FALSE
     }
@@ -252,7 +252,7 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU
 
     # fit statistics
     res0 <- btm_fit_statistics( probs=probs, dat0=dat0, ind1=ind1, ind2=ind2,
-                TP=TP, judge=judge, wgt.ties=wgt.ties )
+                    TP=TP, judge=judge, wgt.ties=wgt.ties )
     effects$outfit <- res0$outfit
     effects$infit <- res0$infit
     multiple_judges <- res0$multiple_judges
@@ -267,10 +267,10 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU
     #--- output list
     effects <- effects[ order(effects$propscore, decreasing=TRUE), ]
     res <- list( effects=effects, pars=pars, summary.effects=summary.effects,
-                mle.rel=mle.rel, sepG=sep.rel, probs=probs, data=dat0,
-                multiple_judges=multiple_judges, fit_judges=fit_judges,
-                residuals=residuals, eps=eps, ignore.ties=ignore.ties,
-                wgt.ties=wgt.ties, time_alg=time_alg, ll=ll, dat=dat)
+                    mle.rel=mle.rel, sepG=sep.rel, probs=probs, data=dat0,
+                    multiple_judges=multiple_judges, fit_judges=fit_judges,
+                    residuals=residuals, eps=eps, ignore.ties=ignore.ties,
+                    wgt.ties=wgt.ties, time_alg=time_alg, ll=ll, dat=dat)
     res$CALL <- CALL
     res$iter <- iter
     ic <- list( n=length(teams), D=nrow(dat0) )

diff --git a/R/ccov.np.R b/R/ccov.np.R
@@ -1,5 +1,5 @@
 ## File Name: ccov.np.R
-## File Version: 1.221
+## File Version: 1.225
 
 
 #---- nonparametric estimation of conditional covariance
@@ -36,12 +36,12 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200),
     }
     icc_items <- matrix( 0, length(thetagrid), I )
     if ( I >=20 ){
-        display <- seq( 1, I, floor( I/20 ) )[ 2:20 ]
+        display <- seq( 1, I, floor( I/20 ) )[ 2L:20 ]
     } else {
         display <- 20
     }
     i <- 1
-    for ( ii in 1:I ){
+    for ( ii in 1L:I ){
         obs_ii <- ! is.na( data[,ii] )
         x <- score[ obs_ii ]
         y <- data[ obs_ii, ii ]
@@ -62,10 +62,11 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200),
         utils::flush.console()
     }
     # calculation of conditional covariance
-    ccov.table <- data.frame( 'item1ID'=rep( 1:I, I), 'item2ID'=rep( 1:I, each=I ) )
+    ccov.table <- data.frame( item1ID=rep(1L:I, I),
+                                item2ID=rep( 1L:I, each=I) )
     ccov.table <- ccov.table[ ccov.table$item1ID < ccov.table$item2ID, ]
     ccov.table$N <- apply( ccov.table, 1, FUN=function(ll){
-                    sum( rowSums( is.na( data[, c( ll[1], ll[2] ) ] ) )==0 ) } )
+                            sum( rowSums( is.na( data[, c( ll[1], ll[2] ) ] ) )==0 ) } )
     ccov.table <- ccov.table[ ccov.table$N > 0, ]
     ccov.table$item1 <- colnames(data)[ ccov.table$item1ID ]
     ccov.table$item2 <- colnames(data)[ ccov.table$item2ID ]
@@ -76,9 +77,9 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200),
     ccov.matrix <- prod.matrix <- matrix( 0, nrow=length(thetagrid), ncol=FF )
     ii <- 1
     ccov_sum_score <- rep(NA, FF)
-    for (ff in 1:FF){
+    for (ff in 1L:FF){
         if (FF>20){
-            display <- seq( 1, FF, floor( FF/20 ) )[ 2:20 ]
+            display <- seq( 1, FF, floor( FF/20 ) )[ 2L:20 ]
         } else {
             display <- seq(1,FF)
         }

diff --git a/R/ccov_np_compute_ccov_sum_score.R b/R/ccov_np_compute_ccov_sum_score.R
@@ -1,5 +1,5 @@
 ## File Name: ccov_np_compute_ccov_sum_score.R
-## File Version: 0.161
+## File Version: 0.162
 
 ccov_np_compute_ccov_sum_score <- function(score, data, use_rcpp=TRUE)
 {
@@ -8,7 +8,7 @@ ccov_np_compute_ccov_sum_score <- function(score, data, use_rcpp=TRUE)
     NS <- length(scores)
     ccov_ff <- rep(NA,NS)
     if (!use_rcpp){
-        for (ss in 1:NS){
+        for (ss in 1L:NS){
             i1 <- which(score==scores[ss])
             s1 <- stats::cov.wt(x=data[i1,], method='ML')
             ccov_ff[ss] <- s1$cov[1,2]

diff --git a/R/decategorize.R b/R/decategorize.R
@@ -1,5 +1,5 @@
 ## File Name: decategorize.R
-## File Version: 0.122
+## File Version: 0.124
 
 #* decategorize
 decategorize <- function( dat, categ_design=NULL )
@@ -10,9 +10,9 @@ decategorize <- function( dat, categ_design=NULL )
 
     #** handle categories
     if ( ! is.null( dfr ) ){
-        vars <- sort( unique( paste( dfr$variable )))
+        vars <- sort(unique(paste(dfr$variable)))
         VV <- length(vars)
-        for (vv in 1:VV){
+        for (vv in 1L:VV){
             dfr.vv <- dfr[ paste(dfr$variable)==vars[vv], ]
             dat4[, vars[vv] ] <- dfr.vv[ match( dat3[,vars[vv]], dfr.vv$recode ), 'orig']
         }

diff --git a/R/dif.strata.variance.R b/R/dif.strata.variance.R
@@ -1,5 +1,5 @@
 ## File Name: dif.strata.variance.R
-## File Version: 0.151
+## File Version: 0.153
 
 
 
@@ -15,7 +15,7 @@ dif.strata.variance <- function( dif, se.dif, itemcluster )
     stratadif$M.DIF <- stats::aggregate( dif, list(itemcluster), mean, na.rm=TRUE )[,2]
     # DIF in strata
     SS <- nrow(stratadif)
-    for (ss in 1:SS){
+    for (ss in 1L:SS){
         items.ss <- which( itemcluster==stratadif[ss,'strata']  )
         dif.ss <- dif[ items.ss ]
         difv.ss <- dif.variance( dif=dif.ss, se.dif=se.dif[ items.ss ] )
@@ -25,8 +25,9 @@ dif.strata.variance <- function( dif, se.dif, itemcluster )
     stratadif[ is.na(stratadif ) ] <- 0
 
     sd_ni1 <- stratadif$N.Items-1
-    weighted.DIFSD <- sum(stratadif$N.Items/sum(stratadif$N.Items)*stratadif$weighted.tau)
-    unweighted.DIFSD <- sum( sd_ni1/ (sum(stratadif$N.Items)-1)*stratadif$unweighted.tau)
+    NI <- sum(stratadif$N.Items)
+    weighted.DIFSD <- sum(stratadif$N.Items/NI*stratadif$weighted.tau)
+    unweighted.DIFSD <- sum( sd_ni1/(NI-1)*stratadif$unweighted.tau)
 
     #-- output
     res <- list( stratadif=stratadif, weighted.DIFSD=weighted.DIFSD,

diff --git a/R/mcmc.2pno.R b/R/mcmc.2pno.R
@@ -1,5 +1,5 @@
 ## File Name: mcmc.2pno.R
-## File Version: 1.26
+## File Version: 1.271
 ##############################################
 # MCMC estimation 2PNO model
 mcmc.2pno <- function(dat, weights=NULL, burnin=500, iter=1000, N.sampvalues=1000,
@@ -24,7 +24,7 @@ mcmc.2pno <- function(dat, weights=NULL, burnin=500, iter=1000, N.sampvalues=100
     # item parameters in matrix form
     aM <- matrix( a, nrow=N, ncol=I, byrow=TRUE)
     bM <- matrix( b, nrow=N, ncol=I, byrow=TRUE)
-    theta <- qnorm( ( rowMeans( dat0,na.rm=TRUE ) + .01 ) / 1.02 )
+    theta <- stats::qnorm( ( rowMeans( dat0,na.rm=TRUE ) + .01 ) / 1.02 )
     # define lower and upper thresholds
     ZZ <- 1000
     threshlow <- -ZZ + ZZ*dat

diff --git a/R/mgsem.R b/R/mgsem.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem.R
-## File Version: 0.553
+## File Version: 0.554
 
 mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
         estimator="ML", p_me=2, p_pen=1, pen_type="scad",
@@ -35,7 +35,7 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL,
         groups <- names(suffstat)
         G <- length(groups)
         if (is.null(groups)){
-            groups <- paste0('Group',1:G)
+            groups <- paste0('Group',1L:G)
         }
         data_proc <- NULL
         is_data <- FALSE

diff --git a/R/mgsem_bdiag.R b/R/mgsem_bdiag.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_bdiag.R
-## File Version: 0.05
+## File Version: 0.06
 
 mgsem_bdiag <- function(x1, x2)
 {
@@ -8,7 +8,7 @@ mgsem_bdiag <- function(x1, x2)
     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
+    mat[1L:n1,1L:n1] <- x1
+    mat[n1+1L:n2,n1+1L:n2] <- x2
     return(mat)
 }
diff --git a/R/mgsem_cd_opt.R b/R/mgsem_cd_opt.R
@@ -1,5 +1,5 @@
 ## File Name: mgsem_cd_opt.R
-## File Version: 0.177
+## File Version: 0.178
 
 mgsem_cd_opt <- function(x, opt_fun_args, tol=1e-4, eps_approx=1e-20,
                 maxiter=100, h=1e-4, verbose=TRUE, interval_length=0.025,
@@ -29,7 +29,7 @@ mgsem_cd_opt <- function(x, opt_fun_args, tol=1e-4, eps_approx=1e-20,
 
         x0 <- x
 
-        for (pp in 1:NP){
+        for (pp in 1L:NP){
             partable <- cd_fun_args$partable
             x1 <- x[pp]
             interval <- x1 + interval_length*c(-1,1)

diff --git a/R/mgsem_create_index.R b/R/mgsem_create_index.R
@@ -1,12 +1,12 @@
 ## File Name: mgsem_create_index.R
-## File Version: 0.03
+## File Version: 0.04
 
 
 mgsem_create_index <- function(x, all=TRUE, start=0, symm=FALSE, onlydiag=FALSE)
 {
     ND <- prod(dim(x))
     if (start>0){
-        v <- matrix( start + ( 1:ND ) - 1, nrow=dim(x)[1], ncol=dim(x)[2])
+        v <- matrix( start + ( 1L:ND ) - 1, nrow=dim(x)[1], ncol=dim(x)[2])
         if (symm){
             v <- v + t(v)
             v <- as.vector(v)
@@ -23,7 +23,7 @@ mgsem_create_index <- function(x, all=TRUE, start=0, symm=FALSE, onlydiag=FALSE)
     if (onlydiag){
         ND <- dim(x)[1]
         if (start>0){
-            v <- diag( start - 1 + 1:ND )
+            v <- diag( start - 1 + 1L:ND )
         } else {
             v <- diag( rep(1, dim(x)[1]))
         }