From 83f9aed9020f6b61095f572ebc44fbdb2a2c5a12 Mon Sep 17 00:00:00 2001 From: Robitzsch Date: Mon, 28 Nov 2022 16:45:24 +0100 Subject: [PATCH] 3.13-24 --- DESCRIPTION | 4 +- R/IRT.likelihood_sirt.R | 3 +- R/RcppExports.R | 6 +- R/mcmcirt_create_partable_Theta.R | 4 +- R/rasch.copula2.R | 12 +- R/rasch.mml2.R | 9 +- R/rasch_mml2_estep_missing1.R | 11 +- R/rasch_mml2_mstep_missing1.R | 6 +- R/rasch_mml2_mstep_one_step.R | 8 +- R/rasch_mml2_numdiff_index.R | 48 ++++++-- R/rm.sdt.R | 2 +- R/xxirt.R | 2 +- R/xxirt_hessian.R | 140 ++++++++++++++++++++---- R/xxirt_hessian_compute_loglike.R | 11 ++ R/xxirt_mstep_ThetaParameters.R | 4 +- README.md | 4 +- docs/404.html | 2 +- docs/authors.html | 8 +- docs/index.html | 2 +- docs/pkgdown.yml | 2 +- inst/NEWS | 4 +- src/RcppExports.cpp | 22 +++- src/probs_multcat_items_counts_rcpp.cpp | 4 +- src/sirt_rcpp_rm_sdt.cpp | 7 +- src/sirt_rcpp_xxirt.cpp | 50 ++++++++- 25 files changed, 299 insertions(+), 76 deletions(-) create mode 100644 R/xxirt_hessian_compute_loglike.R diff --git a/DESCRIPTION b/DESCRIPTION index 3283b3c2..1293fa54 100644 --- 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] () Maintainer: Alexander Robitzsch Description: diff --git a/R/IRT.likelihood_sirt.R b/R/IRT.likelihood_sirt.R index e89b4664..4b8499d2 100644 --- 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 index abdfa473..b44bb3a2 100644 --- 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 index ac72622e..7a846714 100644 --- 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 index 70c94944..2898303b 100644 --- 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 index 1f5e4dcc..2c2f3f56 100644 --- 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 index adc3e270..36da5392 100644 --- 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 index ddb2b6eb..d2109622 100644 --- 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 index a0bc8012..27bed9ce 100644 --- 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 index 24d08e95..1448caf3 100644 --- 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 index c83e0013..7b8213b6 100644 --- 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 index 710d20c4..b53502cc 100644 --- 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 diff --git a/R/xxirt_hessian.R b/R/xxirt_hessian.R index bd12c0f0..bfdeb37a 100644 --- a/R/xxirt_hessian.R +++ b/R/xxirt_hessian.R @@ -1,9 +1,10 @@ ## File Name: xxirt_hessian.R -## File Version: 0.363 +## File Version: 0.568 #--- computation of hessian matrix xxirt_hessian <- function( object ) { + h <- 1e-5 item_list <- object$item_list items <- object$items Theta <- object$Theta @@ -27,42 +28,133 @@ xxirt_hessian <- function( object ) NPI <- length(par_items) NPT <- length(par_Theta) + #** preliminary computations + probs_items <- xxirt_compute_itemprobs( item_list=item_list, + items=items, Theta=Theta, ncat=ncat, + partable=partable, partable_index=partable_index ) + + par1 <- xxirt_partable_extract_freeParameters( partable=partable ) + + p.xi.aj <- xxirt_compute_likelihood( probs_items=probs_items, dat=dat, + resp_index=resp_index, dat_resp_bool=dat_resp_bool ) + prior_Theta <- xxirt_compute_priorDistribution( Theta=Theta, + customTheta=customTheta, G=G ) + prior1 <- t( prior_Theta[, group ] ) + probs_items0 <- probs_items + p_xi_aj <- p.xi.aj + + #--- compute Hessian matrix + par1 <- xxirt_partable_extract_freeParameters( partable=partable ) + par2 <- xxirt_parTheta_extract_freeParameters( customTheta=customTheta ) + par0 <- par <- c(par1, par2) + + #** detect whether there are item-wise parameters + a1 <- aggregate( partable$itemnr, list(partable$parindex), min ) + a2 <- aggregate( partable$itemnr, list(partable$parindex), max ) + item_wise <- sum(a2[,2]> a1[,2])==0 + #*********************** - fct_irt <- function(x){ + fct_irt <- function(x, prior1, par0, item_wise=FALSE, probs_items0, + p_xi_aj) + { + eps <- 1e-10 + I <- dim(probs_items0)[1] #*** include free paramaters in partable partable <- xxirt_partable_include_freeParameters( partable, x[ 1:NPI ] ) #**** include parameter in customTheta customTheta$par[ customTheta$est ] <- x[ (NPI+1):(NPI+NPT) ] - #*** item probabilities - probs_items <- xxirt_compute_itemprobs( item_list=item_list, - items=items, Theta=Theta, ncat=ncat, - partable=partable, partable_index=partable_index ) + #*** different parameters + ind <- which(abs(x-par0)>eps) + max_ind <- 0 + min_ind <- 0 + IA <- 100 + if (length(ind)>0){ + max_ind <- max(ind) + min_ind <- min(ind) + items_active <- unique(partable[ partable$parindex %in% ind, "itemnr" ]) + if (length(items_active)>0){ + IA <- length(items_active) + } + } + + + #*** compute individual likelihood - p.xi.aj <- xxirt_compute_likelihood( probs_items=probs_items, dat=dat, + # if ( (! item_wise) | (min_ind==0) | (max_ind>=NPI+1 ) ){ + + use_shortcut <- TRUE + use_shortcut <- use_shortcut & (IA <=2) + + # use_shortcut <- FALSE + + if (! use_shortcut){ + if (min_ind <=NPI){ + probs_items <- xxirt_compute_itemprobs( item_list=item_list, + items=items, Theta=Theta, ncat=ncat, + partable=partable, partable_index=partable_index ) + } + p.xi.aj <- xxirt_compute_likelihood( probs_items=probs_items, dat=dat, resp_index=resp_index, dat_resp_bool=dat_resp_bool ) + } else { + itemnr <- partable[ partable$parindex==min_ind, "itemnr" ] + itemnr2 <- partable[ partable$parindex==max_ind, "itemnr" ] + maxK <- dim(probs_items)[2] + TP <- dim(probs_items)[3] + eps <- 1e-16 + use_itemnr2 <- FALSE + if (length(itemnr2)>0){ + if (itemnr2!=itemnr){ + use_itemnr2 <- TRUE + } + } else { + if (itemnr>1){ + itemnr2 <- 1 + } else { + itemnr2 <- 2 + } + } + + item_index <- c(itemnr, itemnr2) + if (item_index[2]==item_index[1]){ + if (item_index[1]==1){ + item_index[2] <- 2 + } else { + item_index[2] <- 1 + } + } + probs_items_temp <- xxirt_compute_itemprobs( item_list=item_list, + items=items, Theta=Theta, ncat=ncat, + partable=partable, partable_index=partable_index, + item_index=item_index) + probs_ratio1 <- probs_items0 + + probs_ratio1[item_index,,] <- probs_items_temp/( probs_items0[item_index,,]+eps) + probs_ratio1 <- matrix( probs_ratio1, nrow=I, ncol=maxK*TP ) + + p.xi.aj <- sirt_rcpp_xxirt_hessian_reduced_probs(dat=dat, + dat_resp_bool=dat_resp_bool, probs_ratio=probs_ratio1, + TP=TP, maxK=maxK, itemnr=itemnr-1, itemnr2=itemnr2-1, + use_itemnr2=use_itemnr2, p_xi_aj=p_xi_aj) + + } + + #*** compute prior distribution - prior_Theta <- xxirt_compute_priorDistribution( Theta=Theta, + if (max_ind >=NPI+1){ + prior_Theta <- xxirt_compute_priorDistribution( Theta=Theta, customTheta=customTheta, G=G ) - #*** compute posterior distribution and expected counts - res <- xxirt_compute_posterior( prior_Theta=prior_Theta, p.xi.aj=p.xi.aj, - group=group,G=G, weights=weights, dat1=dat1, - dat_resp=dat_resp, maxK=maxK,group_index=group_index, - dat1_resp=dat1_resp) - n.ik <- res$n.ik - p.aj.xi <- res$p.aj.xi - N.ik <- res$N.ik - N.k <- res$N.k - post_unnorm <- res$post_unnorm - dev <- sum( weights * log( rowSums( post_unnorm ) ) ) + prior1 <- t( prior_Theta[, group ] ) + } + dev <- xxirt_hessian_compute_loglike(p.xi.aj=p.xi.aj, prior1=prior1, + weights=weights) return(dev) } #******************************* - #--- compute Hessian matrix - par1 <- xxirt_partable_extract_freeParameters( partable=partable ) - par2 <- xxirt_parTheta_extract_freeParameters( customTheta=customTheta ) - par <- c(par1, par2) - hess <- CDM::numerical_Hessian( par=par, FUN=fct_irt ) + + hess <- CDM::numerical_Hessian( par=par, FUN=fct_irt, h=h, prior1=prior1, + par0=par, item_wise=item_wise, probs_items0=probs_items0, + p_xi_aj=p_xi_aj) rownames(hess) <- colnames(hess) <- names(par) return(hess) } diff --git a/R/xxirt_hessian_compute_loglike.R b/R/xxirt_hessian_compute_loglike.R new file mode 100644 index 00000000..2f50e578 --- /dev/null +++ b/R/xxirt_hessian_compute_loglike.R @@ -0,0 +1,11 @@ +## File Name: xxirt_hessian_compute_loglike.R +## File Version: 0.03 + + +xxirt_hessian_compute_loglike <- function(p.xi.aj, prior1, weights) +{ + post_unnorm <- prior1 * p.xi.aj + # compute log-likelihood + dev <- sum( weights * log( rowSums( post_unnorm ) ) ) + return(dev) +} diff --git a/R/xxirt_mstep_ThetaParameters.R b/R/xxirt_mstep_ThetaParameters.R index 14568662..1c5d50f6 100644 --- a/R/xxirt_mstep_ThetaParameters.R +++ b/R/xxirt_mstep_ThetaParameters.R @@ -1,5 +1,5 @@ ## File Name: xxirt_mstep_ThetaParameters.R -## File Version: 0.182 +## File Version: 0.183 ####################################################################### xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, @@ -14,7 +14,7 @@ xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, NP <- length(customTheta$prior) pen <- 0 if ( NP > 0 ){ - for (pp in 1:NP){ + for (pp in 1:NP){ if ( ! is.na( customTheta$prior[pp] ) ) { prior_pp <- customTheta$prior[pp] # val <- par1[ names(customTheta$prior) ] diff --git a/README.md b/README.md index 171cd3a5..cf077bbd 100644 --- a/README.md +++ b/README.md @@ -25,9 +25,9 @@ The CRAN version can be installed from within R using: utils::install.packages("sirt") ``` -#### GitHub version `sirt` 3.13-11 (2022-10-11) +#### GitHub version `sirt` 3.13-24 (2022-11-28) -[![](https://img.shields.io/badge/github%20version-3.13--11-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   +[![](https://img.shields.io/badge/github%20version-3.13--24-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   The version hosted [here](https://github.com/alexanderrobitzsch/sirt) is the development version of `sirt`. The GitHub version can be installed using `devtools` as: diff --git a/docs/404.html b/docs/404.html index a3832cbc..a47d889b 100644 --- a/docs/404.html +++ b/docs/404.html @@ -71,7 +71,7 @@ sirt - 3.13-11 + 3.13-24 diff --git a/docs/authors.html b/docs/authors.html index fb3cb7d4..aa93f364 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -71,7 +71,7 @@ sirt - 3.13-11 + 3.13-24 @@ -111,12 +111,12 @@

Citation

Source: inst/CITATION -

Robitzsch, A. (2022). sirt: Supplementary Item Response Theory Models. R package version 3.13-11. https://CRAN.R-project.org/package=sirt

- 
@Manual{sirt_3.13-11,
+    
Robitzsch, A. (2022). sirt: Supplementary Item Response Theory Models. R package version 3.13-24. https://CRAN.R-project.org/package=sirt

+    
@Manual{sirt_3.13-24,
   title = {sirt: Supplementary Item Response Theory Models},
   author = {Alexander Robitzsch},
   year = {2022},
-  note = {R package version 3.13-11},
+  note = {R package version 3.13-24},
   url = {https://CRAN.R-project.org/package=sirt},
 }

 
diff --git a/docs/index.html b/docs/index.html
index 768d729d..040dd30f 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -47,7 +47,7 @@
       
       
         sirt
-        3.13-11
+        3.13-24
       
     
 
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index ea84c57c..7af75ca5 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -2,5 +2,5 @@ pandoc: 1.13.1
 pkgdown: 1.5.1
 pkgdown_sha: ~
 articles: []
-last_built: 2022-10-11T21:09Z
+last_built: 2022-11-28T15:27Z
 
diff --git a/inst/NEWS b/inst/NEWS
index 04f99c4d..89b2dbc8 100644
--- a/inst/NEWS
+++ b/inst/NEWS
@@ -35,7 +35,7 @@ CHANGELOG sirt
 
 
 ------------------------------------------------------------------------
-VERSIONS sirt 3.13  |  2022-10-11  |  Last: sirt 3.13-11
+VERSIONS sirt 3.13  |  2022-11-28  |  Last: sirt 3.13-24
 ------------------------------------------------------------------------
 
 NOTE  * included arguments 'pw_linear', 'pw_quadratic' that allow 
@@ -53,6 +53,8 @@ NOTE  * included argument 'np_fun_item' in xxirt() that counts the
         number of item parameters
 FIXED * fixed a bug in cfa_meas_inv() (thanks to @mbluemke,
         https://github.com/alexanderrobitzsch/sirt/issues/19)
+FIXED * fixed a bug in rm.sdt() caused by numeric instabilities in
+        deviance computation (thanks to Yusuke Shono)
 
 
 DATA  * included/modified datasets: ---
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index e0cd792a..ce05ffa2 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,5 +1,5 @@
 //// File Name: RcppExports.cpp
-//// 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
 
@@ -1531,6 +1531,25 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// sirt_rcpp_xxirt_hessian_reduced_probs
+Rcpp::NumericMatrix sirt_rcpp_xxirt_hessian_reduced_probs(Rcpp::IntegerMatrix dat, Rcpp::LogicalMatrix dat_resp_bool, Rcpp::NumericMatrix probs_ratio, int TP, int maxK, int itemnr, int itemnr2, bool use_itemnr2, Rcpp::NumericMatrix p_xi_aj);
+RcppExport SEXP _sirt_sirt_rcpp_xxirt_hessian_reduced_probs(SEXP datSEXP, SEXP dat_resp_boolSEXP, SEXP probs_ratioSEXP, SEXP TPSEXP, SEXP maxKSEXP, SEXP itemnrSEXP, SEXP itemnr2SEXP, SEXP use_itemnr2SEXP, SEXP p_xi_ajSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::IntegerMatrix >::type dat(datSEXP);
+    Rcpp::traits::input_parameter< Rcpp::LogicalMatrix >::type dat_resp_bool(dat_resp_boolSEXP);
+    Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type probs_ratio(probs_ratioSEXP);
+    Rcpp::traits::input_parameter< int >::type TP(TPSEXP);
+    Rcpp::traits::input_parameter< int >::type maxK(maxKSEXP);
+    Rcpp::traits::input_parameter< int >::type itemnr(itemnrSEXP);
+    Rcpp::traits::input_parameter< int >::type itemnr2(itemnr2SEXP);
+    Rcpp::traits::input_parameter< bool >::type use_itemnr2(use_itemnr2SEXP);
+    Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type p_xi_aj(p_xi_ajSEXP);
+    rcpp_result_gen = Rcpp::wrap(sirt_rcpp_xxirt_hessian_reduced_probs(dat, dat_resp_bool, probs_ratio, TP, maxK, itemnr, itemnr2, use_itemnr2, p_xi_aj));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_sirt_gooijer_csn_table", (DL_FUNC) &_sirt_gooijer_csn_table, 7},
@@ -1632,6 +1651,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_sirt_sirt_rcpp_rm_sdt_calc_probs_grm_item", (DL_FUNC) &_sirt_sirt_rcpp_rm_sdt_calc_probs_grm_item, 8},
     {"_sirt_sirt_rcpp_xxirt_compute_posterior_expected_counts", (DL_FUNC) &_sirt_sirt_rcpp_xxirt_compute_posterior_expected_counts, 3},
     {"_sirt_sirt_rcpp_xxirt_compute_likelihood", (DL_FUNC) &_sirt_sirt_rcpp_xxirt_compute_likelihood, 5},
+    {"_sirt_sirt_rcpp_xxirt_hessian_reduced_probs", (DL_FUNC) &_sirt_sirt_rcpp_xxirt_hessian_reduced_probs, 9},
     {NULL, NULL, 0}
 };
 
diff --git a/src/probs_multcat_items_counts_rcpp.cpp b/src/probs_multcat_items_counts_rcpp.cpp
index 3c5c5ea8..ca890316 100644
--- a/src/probs_multcat_items_counts_rcpp.cpp
+++ b/src/probs_multcat_items_counts_rcpp.cpp
@@ -1,5 +1,5 @@
 //// File Name: probs_multcat_items_counts_rcpp.cpp
-//// File Version: 3.175
+//// File Version: 3.178
 
 
 #include 
@@ -258,3 +258,5 @@ Rcpp::List sirt_rcpp_calccounts_pcm_groups_C( Rcpp::IntegerMatrix dat,
         );
 }
 ///********************************************************************
+
+
diff --git a/src/sirt_rcpp_rm_sdt.cpp b/src/sirt_rcpp_rm_sdt.cpp
index fcd2d694..5d5e9b13 100644
--- a/src/sirt_rcpp_rm_sdt.cpp
+++ b/src/sirt_rcpp_rm_sdt.cpp
@@ -1,5 +1,5 @@
 //// File Name: sirt_rcpp_rm_sdt.cpp
-//// File Version: 0.415
+//// File Version: 0.416
 
 
 // [[Rcpp::depends(RcppArmadillo)]]
@@ -54,7 +54,8 @@ Rcpp::List sirt_rcpp_rm_sdt_posterior( Rcpp::NumericVector prob_item,
     double p1=0;
     double temp=0;
     double ll=0;
-
+    double eps_like=1e-300;
+    
     for (int nn=0; nn