3.13-24

alexanderrobitzsch · Nov 28, 2022 · 83f9aed · 83f9aed
1 parent d0afec2
commit 83f9aed
Show file tree

Hide file tree

Showing 25 changed files with 299 additions and 76 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.13-11
-Date: 2022-10-11 22:44:04
+Version: 3.13-24
+Date: 2022-11-28 16:03:25
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

diff --git a/R/IRT.likelihood_sirt.R b/R/IRT.likelihood_sirt.R
@@ -1,5 +1,5 @@
 ## File Name: IRT.likelihood_sirt.R
-## File Version: 0.17
+## File Version: 0.18
 
 
 ################################################
@@ -20,6 +20,7 @@
 
 ########################################################
 # likelihood rasch.copula2
+# @ checked ARb 2022-10-13
 IRT.likelihood.rasch.copula2 <- function( object, ... )
 {
     ll <- object$f.yi.qk

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.013011
+## File Version: 3.013024
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
@@ -399,3 +399,7 @@ sirt_rcpp_xxirt_compute_likelihood <- function(dat, dat_resp_bool, probs, TP, ma
     .Call('_sirt_sirt_rcpp_xxirt_compute_likelihood', PACKAGE='sirt', dat, dat_resp_bool, probs, TP, maxK)
 }
 
+sirt_rcpp_xxirt_hessian_reduced_probs <- function(dat, dat_resp_bool, probs_ratio, TP, maxK, itemnr, itemnr2, use_itemnr2, p_xi_aj) {
+    .Call('_sirt_sirt_rcpp_xxirt_hessian_reduced_probs', PACKAGE='sirt', dat, dat_resp_bool, probs_ratio, TP, maxK, itemnr, itemnr2, use_itemnr2, p_xi_aj)
+}
+
diff --git a/R/mcmcirt_create_partable_Theta.R b/R/mcmcirt_create_partable_Theta.R
@@ -1,8 +1,8 @@
 ## File Name: mcmcirt_create_partable_Theta.R
-## File Version: 0.12
+## File Version: 0.13
 
 
-mcmcirt_create_partable_Theta <- function(par, est, prior=NULL, 
+mcmcirt_create_partable_Theta <- function(par, est, prior=NULL,
         prior_par1=NULL, prior_par2=NULL, sd_proposal=NULL)
 {
     NP <- length(par)

diff --git a/R/rasch.copula2.R b/R/rasch.copula2.R
@@ -1,5 +1,5 @@
 ## File Name: rasch.copula2.R
-## File Version: 6.317
+## File Version: 6.319
 
 
 
@@ -290,11 +290,14 @@ rasch.copula2 <- function( dat, itemcluster, weights=NULL,
                 rescop <- .ll.rasch.copula20( theta.k, b0, alpha1, alpha2, a, dat2.li, itemcluster0,
                             CC, dp.ld, dat2.ld, dat3.ld, dat2.ld.resp, dat2.li.resp, delta, wgt.theta, I,
                             bdat2.li, bdat2.li.resp, pattern, GG, copula.type,
-                            Ncat.ld    )
+                            Ncat.ld)
                 res.posterior <- rescop
-                            }
+            }
             # is this really necessary?
 #            wgt.theta <- rescop$pik
+Revalprstr("res.posterior")
+
+stop()
 
             rest1 <- .update.ll.rasch.copula21( theta.k, b0 + h*est.bb, alpha1, alpha2,
                 a, dat2.li, itemcluster0,
@@ -694,8 +697,7 @@ rasch.copula2 <- function( dat, itemcluster, weights=NULL,
 
     #---- results item parameters
     N_item <- colSums((!is.na(dat00))*w)
-    item <- data.frame( item=colnames(dat),
-                N=N_item,
+    item <- data.frame( item=colnames(dat), N=N_item,
                 p=colSums( w*dat00, na.rm=TRUE )/N_item,
                 b=b, est.b=est.b, a=a, est.a=est.a )
     item$thresh <- item$a * item$b

diff --git a/R/rasch.mml2.R b/R/rasch.mml2.R
@@ -1,5 +1,5 @@
 ## File Name: rasch.mml2.R
-## File Version: 7.688
+## File Version: 7.699
 
 
 # Semiparametric Maximum Likelihood Estimation in the Rasch type Model
@@ -404,6 +404,7 @@ rasch.mml2 <- function( dat, theta.k=seq(-6,6,len=21), group=NULL, weights=NULL,
     #-- indicators of estimated parameters
     est_parameters <- list( a=sum(est.a)>0, c=sum(est.c)>0, d=sum(est.d)>0)
 
+    aa0 <- Sys.time()
 
     #******************************************************
     #*************** MML Iteration Algorithm **************
@@ -523,6 +524,7 @@ rasch.mml2 <- function( dat, theta.k=seq(-6,6,len=21), group=NULL, weights=NULL,
         }
 # cat("m step") ; zz1 <- Sys.time(); print(zz1-zz0) ; zz0 <- zz1
 
+
         #***************************************
         # update mean and covariance in multidimensional models
         if ( D > 1){
@@ -669,8 +671,7 @@ rasch.mml2 <- function( dat, theta.k=seq(-6,6,len=21), group=NULL, weights=NULL,
                                     }
                             }
                          }  # end non-normal distribution
-  #cat("trait distribution estimation") ; zz1 <- Sys.time(); print(zz1-zz0) ; zz0 <- zz1
-
+#  cat("trait distribution estimation") ; zz1 <- Sys.time(); print(zz1-zz0) ; zz0 <- zz1
 
         #---- estimation of alpha, c and d parameters
         alpha.change <- 0
@@ -893,6 +894,8 @@ rasch.mml2 <- function( dat, theta.k=seq(-6,6,len=21), group=NULL, weights=NULL,
         ############################################################################
         ##**************************************************************************
 
+# cat(" ++++ total estimation time") ; aa1 <- Sys.time(); print(aa1-aa0) ; aa0 <- aa1
+
             if ( irtmodel %in% irtmodel_missing){
                 m1$center <- FALSE
                 G <- 1

diff --git a/R/rasch_mml2_estep_missing1.R b/R/rasch_mml2_estep_missing1.R
@@ -1,5 +1,5 @@
 ## File Name: rasch_mml2_estep_missing1.R
-## File Version: 1.111
+## File Version: 1.119
 
 
 
@@ -11,19 +11,24 @@ rasch_mml2_estep_missing1 <- function( dat2, dat2.resp, theta.k, b, beta, delta.
     if (is.null(est.a)){
         est.a <- rep(0,I)
     }
-a0 <- Sys.time()
+# a0 <- Sys.time()
     # probability correct response
     pjk <- rasch_mml2_calcprob_missing1( theta.k=theta.k, b=b, beta=beta,
                 delta.miss=delta.miss, pjk=pjk, fixed.a=fixed.a, irtmodel=irtmodel )
+# cat( " *** calc probs") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
+
 
     #** calculate likelihood
     probs_ <- as.matrix( array( pjk, dim=c(I,CC*TP) ) )
     f.yi.qk <- sirt_rcpp_probs_pcm_groups_C( dat=dat2, dat_resp=dat2.resp, group=group_,
                     probs=probs_, CC=CC, TP=TP )
+#  cat( " *** calc like") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
 
     #*** calculate expected counts
     e1 <- sirt_rcpp_calccounts_pcm_groups_C( dat=dat2, dat_resp=dat2.resp, group=group_,
                 fyiqk=f.yi.qk, pik=pi.k, CC=CC, weights=weights )
+# cat( " *** calc counts") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
+
     e1$f.yi.qk <- f.yi.qk
     v1 <- array( e1$nik, dim=c(I,CC,TP) )
     e1$pjk <- pjk
@@ -37,6 +42,6 @@ a0 <- Sys.time()
 }
 
 
-# cat( " * calc counts") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
+
 
 .e.step.missing1 <- rasch_mml2_estep_missing1
diff --git a/R/rasch_mml2_mstep_missing1.R b/R/rasch_mml2_mstep_missing1.R
@@ -1,5 +1,5 @@
 ## File Name: rasch_mml2_mstep_missing1.R
-## File Version: 0.156
+## File Version: 0.159
 
 
 #*** M-step for missing data model
@@ -65,7 +65,7 @@ rasch_mml2_mstep_missing1 <- function( theta.k, n.ik, mitermax, conv1,
                         n.ik=n.ik, diffindex=diffindex, max.increment=max.increment,
                         numdiff.parm=numdiff.parm)
         args0 <- res$args0
-        beta <- args0[[entry]]
+        beta <- sirt_squeeze(args0[[entry]], lower=min.beta)
         max_incr_beta <- res$max_increment
         se.beta <- res$se
 
@@ -84,7 +84,7 @@ rasch_mml2_mstep_missing1 <- function( theta.k, n.ik, mitermax, conv1,
             max_incr_delta <- res$max_increment
             se.delta <- res$se
         }
-        #--- update delta
+        #--- update a
         if (est_a){
             a0 <- fixed.a
             max.increment <- max_incr_a

diff --git a/R/rasch_mml2_mstep_one_step.R b/R/rasch_mml2_mstep_one_step.R
@@ -1,10 +1,13 @@
 ## File Name: rasch_mml2_mstep_one_step.R
-## File Version: 1.071
+## File Version: 1.076
 
 
 rasch_mml2_mstep_one_step <- function(args0, prob_fun, entry, n.ik,
         diffindex, max.increment, numdiff.parm)
 {
+
+    # gg0 <- Sys.time()
+
     h <- numdiff.parm
     val0 <- args0[[entry]]
     pjk <- do.call(what=prob_fun, args=args0)
@@ -17,9 +20,12 @@ rasch_mml2_mstep_one_step <- function(args0, prob_fun, entry, n.ik,
     res <- rasch_mml2_numdiff_index( pjk=pjk, pjk1=pjk1, pjk2=pjk2, n.ik=n.ik,
                     diffindex=diffindex, max.increment=max.increment,
                     numdiff.parm=numdiff.parm )
+
     # update
     args0[[entry]] <- args0[[entry]] + res$increment
     res$args0 <- args0
     #- output
     return(res)
 }
+
+# cat( " @@@@@ calc probs") ; gg1 <- Sys.time(); print(gg1-gg0) ; gg0 <- gg1
diff --git a/R/rasch_mml2_numdiff_index.R b/R/rasch_mml2_numdiff_index.R
@@ -1,26 +1,52 @@
 ## File Name: rasch_mml2_numdiff_index.R
-## File Version: 1.112
+## File Version: 1.127
 
 
 
 #** general function for numerical differentiation
 #** diffindex aggregates across super items
 rasch_mml2_numdiff_index <- function( pjk, pjk1, pjk2, n.ik, diffindex,
-        max.increment, numdiff.parm, eps=1e-16 )
+        max.increment, numdiff.parm, eps=1e-16, shortcut=TRUE )
 {
     eps2 <- 1e-10
     h <- numdiff.parm
     an.ik <- n.ik
+
     ll0 <- rowSums( an.ik * log(pjk+eps) )
-    ll1 <- rowSums( an.ik * log(pjk1+eps) )
-    ll2 <- rowSums( an.ik * log(pjk2+eps) )
-    ll0 <- stats::aggregate( ll0, list(diffindex), sum )[,2]
-    ll1 <- stats::aggregate( ll1, list(diffindex), sum )[,2]
-    ll2 <- stats::aggregate( ll2, list(diffindex), sum )[,2]
-    d1 <- ( ll1 - ll2  ) / ( 2 * h )    # negative sign?
-    # second order derivative
-    # f(x+h)+f(x-h)=2*f(x) + f''(x)*h^2
-    d2 <- ( ll1 + ll2 - 2*ll0 ) / h^2
+    ll0 <- rowsum(ll0, diffindex)[,1]
+
+    if (! shortcut){
+
+        ll1 <- rowSums( an.ik * log(pjk1+eps) )
+        ll2 <- rowSums( an.ik * log(pjk2+eps) )
+        # ll0 <- stats::aggregate( ll0, list(diffindex), sum )[,2]
+        # ll1 <- stats::aggregate( ll1, list(diffindex), sum )[,2]
+        # ll2 <- stats::aggregate( ll2, list(diffindex), sum )[,2]
+        ll1 <- rowsum(ll1, diffindex)[,1]
+        ll2 <- rowsum(ll2, diffindex)[,1]
+
+        d1 <- ( ll1 - ll2  ) / ( 2 * h )    # negative sign?
+
+        # second order derivative
+        # f(x+h)+f(x-h)=2*f(x) + f''(x)*h^2
+        d2 <- ( ll1 + ll2 - 2*ll0 ) / h^2
+
+    }
+
+    if (shortcut){
+
+        p1 <- ( pjk1 - pjk2 ) / (2*h)
+        ll1a <- rowSums( an.ik * p1 / pjk )
+        ll1a <- rowsum(ll1a, diffindex)[,1]
+        d1 <- ll1a
+
+        p2 <- ( pjk1 + pjk2 - 2*pjk) / h^2
+        ll1a <- rowSums( an.ik * (p2*pjk-p1^2)/pjk^2 )
+        ll1a <- rowsum(ll1a, diffindex)[,1]
+        d2 <- ll1a
+
+    }
+
     # change in item difficulty
     d2[ abs(d2) < eps2 ] <- eps2
     increment <- - d1 / d2

diff --git a/R/rm.sdt.R b/R/rm.sdt.R
@@ -1,5 +1,5 @@
 ## File Name: rm.sdt.R
-## File Version: 8.886
+## File Version: 8.888
 
 #################################################################
 # Hierarchical rater model

diff --git a/R/xxirt.R b/R/xxirt.R
@@ -1,5 +1,5 @@
 ## File Name: xxirt.R
-## File Version: 0.947
+## File Version: 0.950
 
 
 #--- user specified item response model