3.13-105

alexanderrobitzsch · Mar 19, 2023 · 9dd6b9f · 9dd6b9f
1 parent 6d3bde5
commit 9dd6b9f
Show file tree

Hide file tree

Showing 28 changed files with 118 additions and 91 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.13-96
-Date: 2023-03-17 10:17:36
+Version: 3.13-105
+Date: 2023-03-19 12:31:21
 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: 3.013096
+## File Version: 3.013105
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 

diff --git a/R/linking_haberman_itempars_prepare.R b/R/linking_haberman_itempars_prepare.R
@@ -1,5 +1,5 @@
 ## File Name: linking_haberman_itempars_prepare.R
-## File Version: 0.12
+## File Version: 0.131
 
 
 linking_haberman_itempars_prepare <- function(b, a=NULL, wgt=NULL)
@@ -8,7 +8,7 @@ linking_haberman_itempars_prepare <- function(b, a=NULL, wgt=NULL)
     NS <- ncol(b)
     if ( is.null(rownames(b) )){
         i0 <- ceiling( log10(I) + 1 )
-        rownames(b) <- paste0("I", 10^i0 + 1:I)
+        rownames(b) <- paste0('I', 10^i0 + 1:I)
     }
     if ( is.null(colnames(b) )){
         colnames(b) <- 1:NS

diff --git a/R/linking_haberman_summary_estimation_information.R b/R/linking_haberman_summary_estimation_information.R
@@ -1,11 +1,11 @@
 ## File Name: linking_haberman_summary_estimation_information.R
-## File Version: 0.04
+## File Version: 0.051
 
 linking_haberman_summary_estimation_information <- function(res_opt)
 {
-    cat("Estimation type", "=", res_opt$estimation,"\n")
-    cat("Number of iterations", "=", res_opt$iter,"\n")
-    cat("Used trimming factor ('BSQ','HUB')", "=", res_opt$cutoff,"\n")
-    cat("Trimming factor estimated ('BSQ','HUB')", "=", res_opt$k_estimate,"\n")
-    cat("Proportion retained observation ('LTS')", "=", res_opt$lts_prop,"\n")
+    cat('Estimation type', '=', res_opt$estimation,'\n')
+    cat('Number of iterations', '=', res_opt$iter,'\n')
+    cat('Used trimming factor (\'BSQ\',\'HUB\')', '=', res_opt$cutoff,'\n')
+    cat('Trimming factor estimated (\'BSQ\',\'HUB\')', '=', res_opt$k_estimate,'\n')
+    cat('Proportion retained observation (\'LTS\')', '=', res_opt$lts_prop,'\n')
 }
diff --git a/R/linking_haebara_gradient_function_R.R b/R/linking_haebara_gradient_function_R.R
@@ -1,5 +1,5 @@
 ## File Name: linking_haebara_gradient_function_R.R
-## File Version: 0.292
+## File Version: 0.293
 
 
 linking_haebara_gradient_function_R <- function(NI, NS, dist, aM, bM, theta,
@@ -19,15 +19,15 @@ linking_haebara_gradient_function_R <- function(NI, NS, dist, aM, bM, theta,
                 b_exp <- ( b[ii] - mu[ss] ) / sigma[ss]
                 p_exp <- stats::plogis( a_exp * (theta - b_exp ) )
                 der <- p_exp*(1-p_exp)
-                if (dist=="L2"){
+                if (dist=='L2'){
                     der_basis <- -2*(p_obs - p_exp)*prob_theta*der
                 }
-                if (dist=="L1"){
+                if (dist=='L1'){
                     diff2 <- p_obs - p_exp
                     dist2 <- diff2^2
                     der_basis <- -(dist2+eps)^(-.5)*diff2*prob_theta*der
                 }
-                if (dist=="Lp"){
+                if (dist=='Lp'){
                     diff2 <- p_obs - p_exp
                     dist2 <- diff2^2
                     der_basis <- -pow*(dist2+eps)^(pow/2-1)*diff2*prob_theta*der

diff --git a/R/linking_haebara_optim_function_R.R b/R/linking_haebara_optim_function_R.R
@@ -1,5 +1,5 @@
 ## File Name: linking_haebara_optim_function_R.R
-## File Version: 0.170
+## File Version: 0.171
 
 
 linking_haebara_optim_function_R <- function(NI, NS, dist, aM, bM, theta,
@@ -19,13 +19,13 @@ linking_haebara_optim_function_R <- function(NI, NS, dist, aM, bM, theta,
                 ## a_exp[ii]=a[ii]*sigma[ss]
                 ## a_exp[ii] * b_exp[ii]=a[ii]*(b[ii]-mu[ss])
                 dist2 <- (p_obs - p_exp)^2
-                if (dist=="L2"){
+                if (dist=='L2'){
                     dist1 <- sum( dist2*prob_theta )
                 }
-                if (dist=="L1"){
+                if (dist=='L1'){
                     dist1 <- sum( sqrt( dist2 + eps )*prob_theta )
                 }
-                if (dist=="Lp"){
+                if (dist=='Lp'){
                     dist1 <- sum( ( dist2 + eps )^(pow/2)*prob_theta )
                 }
                 val <- val + wgtM[ii,ss]*dist1

diff --git a/R/linking_haebara_summary_optimization.R b/R/linking_haebara_summary_optimization.R
@@ -1,15 +1,15 @@
 ## File Name: linking_haebara_summary_optimization.R
-## File Version: 0.08
+## File Version: 0.091
 
 linking_haebara_summary_optimization <- function(object, digits)
 {
-    cat("Distance function type", "=", object$dist, "\n" )
-    if (object$dist=="L1"){
-        cat("Epsilon Value", "=", object$eps, "\n" )
+    cat('Distance function type', '=', object$dist, '\n')
+    if (object$dist=='L1'){
+        cat('Epsilon Value', '=', object$eps, '\n')
     }
-    cat("Optimization Function Value", "=", round(object$res_optim$value, digits), "\n" )
-    cat("Optimizer", "=", object$res_optim$optimizer, "\n" )
-    cat("use_rcpp", "=", object$use_rcpp, "\n" )
-    cat("Number of iterations", "=", object$res_optim$iter, "\n" )
-    cat("Converged", "=", object$res_optim$converged, "\n" )
+    cat('Optimization Function Value', '=', round(object$res_optim$value, digits), '\n')
+    cat('Optimizer', '=', object$res_optim$optimizer, '\n')
+    cat('use_rcpp', '=', object$use_rcpp, '\n')
+    cat('Number of iterations', '=', object$res_optim$iter, '\n')
+    cat('Converged', '=', object$res_optim$converged, '\n')
 }
diff --git a/R/logLik_sirt.R b/R/logLik_sirt.R
@@ -1,45 +1,42 @@
 ## File Name: logLik_sirt.R
-## File Version: 0.14
+## File Version: 0.152
 
 
-###############################################################
-# log-likelihood function rasch.copula2
+#--- log-likelihood function rasch.copula2
 logLik.rasch.copula2 <- function (object, ...)
 {
     # extract log-likelihood
     out <- - object$ic$deviance / 2
     # number of parameters
-    attr(out, "df") <- object$ic$np
+    attr(out, 'df') <- object$ic$np
     # extract number of observations
-    attr(out, "nobs") <- object$ic$n
-    class(out) <- "logLik"
+    attr(out, 'nobs') <- object$ic$n
+    class(out) <- 'logLik'
     return(out)
 }
 logLik.rasch.copula3 <- logLik.rasch.copula2
-################################################################
 
-#####################################################
-# logLik.rasch.mml
+#-- logLik.rasch.mml
 logLik.rasch.mml <- logLik.rasch.copula2
-#####################################################
-# smirt
+
+#-- smirt
 logLik.smirt <- logLik.rasch.copula2
-#####################################################
-# rasch.mirtlc
+
+#-- rasch.mirtlc
 logLik.rasch.mirtlc <- logLik.rasch.copula2
-#####################################################
-# gom
+
+#-- gom
 logLik.gom <- logLik.rasch.copula2
-#####################################################
-# rm.facets
+
+#-- rm.facets
 logLik.rm.facets <- logLik.rasch.copula2
-#####################################################
-# rm.sdt
+
+#-- rm.sdt
 logLik.rm.sdt <- logLik.rasch.copula2
-#####################################################
-# prob.guttman
+
+#-- prob.guttman
 logLik.prob.guttman <- logLik.rasch.copula2
-#####################################################
-# xxirt
+
+#-- xxirt
 logLik.xxirt <- logLik.rasch.copula2
-#####################################################
+
diff --git a/R/lsem.bootstrap.R b/R/lsem.bootstrap.R
@@ -1,5 +1,5 @@
 ## File Name: lsem.bootstrap.R
-## File Version: 0.347
+## File Version: 0.351
 
 
 lsem.bootstrap <- function(object, R=100, verbose=TRUE, cluster=NULL,
@@ -104,6 +104,7 @@ lsem.bootstrap <- function(object, R=100, verbose=TRUE, cluster=NULL,
     object$class_boot <- TRUE
     object$fitstats_joint <- fitstats_joint
     object$repl_design <- repl_design
+    object$repl_factor <- repl_factor
     object$repl_design_used <- repl_design_used
     s2 <- Sys.time()
     object$s1 <- s1

diff --git a/R/lsem.estimate.R b/R/lsem.estimate.R
@@ -1,5 +1,5 @@
 ## File Name: lsem.estimate.R
-## File Version: 1.057
+## File Version: 1.067
 
 # estimate LSEM model
 lsem.estimate <- function( data, moderator, moderator.grid,
@@ -59,9 +59,9 @@ lsem.estimate <- function( data, moderator, moderator.grid,
 
     # residualize input data
     out <- lsem_residualize( data=data, moderator=moderator,
-                    moderator.grid=moderator.grid,
-                    lavmodel=lavmodel, h=h, bw=bw, residualize=residualize, eps=eps,
-                    verbose=verbose, sampling_weights=sampling_weights, kernel=kernel,
+                    moderator.grid=moderator.grid, lavmodel=lavmodel, h=h, bw=bw,
+                    residualize=residualize, eps=eps, verbose=verbose,
+                    sampling_weights=sampling_weights, kernel=kernel,
                     variables_model=variables_model)
     G <- out$G
     data <- out$data
@@ -123,7 +123,9 @@ lsem.estimate <- function( data, moderator, moderator.grid,
                     pw_linear=pw_linear, pw_quadratic=pw_quadratic,
                     se=se, moderator_variable=moderator_variable,
                     loc_linear_smooth=loc_linear_smooth, pd=pd,
-                    has_meanstructure=has_meanstructure, est_DIF=est_DIF, ... )
+                    residualized_intercepts=residualized_intercepts,
+                    has_meanstructure=has_meanstructure, est_DIF=est_DIF,
+                    residualize=residualize, ... )
     dif_effects <- out2$dif_effects
     parameters <- out2$parameters
     is_meanstructure <- out2$is_meanstructure

diff --git a/R/lsem.test.R b/R/lsem.test.R
@@ -1,10 +1,15 @@
 ## File Name: lsem.test.R
-## File Version: 0.118
+## File Version: 0.133
 
 #**** test LSEM model based on bootstrap
 lsem.test <- function( mod, bmod, models=NULL )
 {
     parameters <- mod$parameters
+    repl_factor <- bmod$repl_factor
+    R <- bmod$R
+    if (is.null(repl_factor)){
+        repl_factor <- 1/(R-1)
+    }
 
     parnames <- unique(paste(parameters$par))
     w <- mod$moderator.density$wgt
@@ -31,7 +36,8 @@ lsem.test <- function( mod, bmod, models=NULL )
             ind_pp <- which( parameters$par==pp )
             parameters_pp <- parameters[ind_pp, ]
             theta <- parameters_pp$est
-            V <- stats::cov( t( parameters_boot[ind_pp, ] ) )
+            par_boot_pp <- t( parameters_boot[ind_pp, ] )
+            V <- stats::cov( par_boot_pp )*(R-1)*repl_factor
             M <- TAM::weighted_mean(x=theta, w=w)
             SD <- TAM::weighted_sd(x=theta, w=w)
             dfr1 <- data.frame(par=pp, M=M, SD=SD, chisq=NA, df=NA, p=NA)

diff --git a/R/lsem_estimate_proc_args.R b/R/lsem_estimate_proc_args.R
@@ -1,5 +1,5 @@
 ## File Name: lsem_estimate_proc_args.R
-## File Version: 0.403
+## File Version: 0.406
 
 lsem_estimate_proc_args <- function(lavaan.args, sufficient_statistics,
     pseudo_weights, lavmodel, data, use_lavaan_survey, est_joint=FALSE,
@@ -18,7 +18,7 @@ lsem_estimate_proc_args <- function(lavaan.args, sufficient_statistics,
     #- variables in model
     partable <- sirt_import_lavaan_lavaanify(model=lavmodel)
     variables_model <- intersect( union( partable$lhs, partable$rhs ), colnames(data) )
-    has_meanstructure <- '~' %in% partable$op
+    has_meanstructure <- '~1' %in% paste(partable$op)
 
     #- joint estimation
     par_vec <- union(union(par_invariant, par_linear), par_quadratic)

diff --git a/R/lsem_fitsem.R b/R/lsem_fitsem.R
@@ -1,5 +1,5 @@
 ## File Name: lsem_fitsem.R
-## File Version: 0.623
+## File Version: 0.639
 
 lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.grid,
                 verbose, pars, standardized, variables_model, sufficient_statistics,
@@ -8,7 +8,7 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr
                 partable_joint=NULL, pw_linear=1, pw_quadratic=1,
                 se="standard", moderator_variable=NULL,
                 loc_linear_smooth=NULL, pd=FALSE, has_meanstructure=FALSE,
-                est_DIF=FALSE, ... )
+                est_DIF=FALSE, residualized_intercepts=NULL, residualize=TRUE, ... )
 {
     parameters <- NULL
     fits <- NULL
@@ -33,7 +33,9 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr
                     variables_model=variables_model, weights=weights,
                     moderator_variable=moderator_variable,
                     loc_linear_smooth=loc_linear_smooth, moderator.grid=moderator.grid,
-                    pd=pd)
+                    pd=pd, residualized_intercepts=residualized_intercepts,
+                    has_meanstructure=has_meanstructure,
+                    residualize=residualize)
     }
     if (est_joint & (! sufficient_statistics)){
         N <- nrow(dat)
@@ -53,6 +55,7 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr
                                 par_quadratic=par_quadratic, pw_linear=pw_linear,
                                 pw_quadratic=pw_quadratic)
         }
+
         #- fit model
         fit_args <- list( partable_joint=partable_joint,
                             is_meanstructure=is_meanstructure, data_joint=data_joint,

diff --git a/R/lsem_fitsem_compute_sufficient_statistics.R b/R/lsem_fitsem_compute_sufficient_statistics.R
@@ -1,9 +1,10 @@
 ## File Name: lsem_fitsem_compute_sufficient_statistics.R
-## File Version: 0.094
+## File Version: 0.109
 
 lsem_fitsem_compute_sufficient_statistics <- function(G, dat, variables_model,
     weights, moderator_variable=NULL, loc_linear_smooth=NULL, moderator.grid=NULL,
-    pd=FALSE)
+    pd=FALSE, residualized_intercepts=NULL,    has_meanstructure=FALSE,
+    residualize=TRUE)
 {
     wmean <- wcov <- Nobs <- as.list(1:G)
     data_suff <- dat[, variables_model]
@@ -14,12 +15,21 @@ lsem_fitsem_compute_sufficient_statistics <- function(G, dat, variables_model,
         res <- lsem_weighted_cov( x=data_suff, weights=weights_gg, x_resp=dat_resp,
                     moderator_variable=moderator_variable,
                     loc_linear_smooth=loc_linear_smooth,
-                    moderator_value=moderator.grid[gg], pd=pd)
+                    moderator_value=moderator.grid[gg], pd=pd,
+                    residualized_intercepts=residualized_intercepts,
+                    has_meanstructure=has_meanstructure, residualize=residualize)
         wmean[[gg]] <- res$mean
         wcov[[gg]] <- res$cov
         Nobs[[gg]] <- round(res$Nobs)
     }
 
+    #** adapt if mean structure is requested
+    if ( has_meanstructure & residualize ){
+        for (gg in 1:G){
+            wmean[[gg]] <- residualized_intercepts[gg,]
+        }
+    }
+
     #- output
     res <- list(wmean=wmean, wcov=wcov, Nobs=Nobs)
     return(res)

diff --git a/R/lsem_fitsem_joint_estimation.R b/R/lsem_fitsem_joint_estimation.R
@@ -1,5 +1,5 @@
 ## File Name: lsem_fitsem_joint_estimation.R
-## File Version: 0.181
+## File Version: 0.184
 
 lsem_fitsem_joint_estimation <- function(partable_joint,
     is_meanstructure, sample_stats, lavaan_est_fun, se,
@@ -9,6 +9,7 @@ lsem_fitsem_joint_estimation <- function(partable_joint,
     wmean <- sample_stats$wmean
     wcov <- sample_stats$wcov
     Nobs <- sample_stats$Nobs
+
     if (is_meanstructure){
         sample_mean <- wmean
     } else {

diff --git a/R/lsem_residualize.R b/R/lsem_residualize.R
@@ -1,5 +1,5 @@
 ## File Name: lsem_residualize.R
-## File Version: 0.428
+## File Version: 0.433
 
 
 #**** residualize data
@@ -36,7 +36,6 @@ lsem_residualize <- function( data, moderator, moderator.grid,
     residualized_intercepts <- matrix( 0, nrow=G, ncol=V)
     colnames(residualized_intercepts) <- vars
     rownames(residualized_intercepts) <- round( moderator.grid, 3 )
-
     if (residualize){
         if (verbose){
             cat('** Residualize Data\n')