From 097cd32b5b3bd1e0c030c0c6653ef533628e9e17 Mon Sep 17 00:00:00 2001 From: Robitzsch Date: Fri, 1 Mar 2024 16:11:00 +0100 Subject: [PATCH] 4.2-13 --- DESCRIPTION | 4 +- R/RcppExports.R | 2 +- R/lsem.MGM.stepfunctions.R | 4 +- R/lsem.bootstrap.R | 4 +- R/lsem.estimate.R | 53 ++++--- R/lsem.permutationTest.R | 8 +- R/lsem.test.R | 8 +- R/lsem_aggregate_statistics.R | 14 ++ R/lsem_estimate_proc_args.R | 12 +- R/lsem_fit_initial_model.R | 7 +- ..._fit_initial_model_sufficient_statistics.R | 36 +++-- R/lsem_fitsem.R | 32 +++-- R/lsem_fitsem_compute_sufficient_statistics.R | 59 +++++--- ...fitsem_joint_estimation_prepare_partable.R | 10 +- R/lsem_group_moderator.R | 32 +++-- R/lsem_local_weights.R | 96 +++++++++---- R/lsem_parameter_summary.R | 4 +- R/lsem_residualize.R | 134 ++++++++++++------ R/m_est.R | 10 +- R/rasch.pairwise.itemcluster.R | 11 +- R/summary.lsem.R | 5 +- R/xxirt.R | 14 +- R/xxirt_em_algorithm.R | 6 +- R/xxirt_newton_raphson.R | 8 +- R/xxirt_nr_grad_fun_Rcpp.R | 4 +- R/xxirt_nr_optim_fun.R | 14 +- R/xxirt_proc_ParTable.R | 4 +- README.md | 4 +- inst/NEWS | 8 +- man/lsem.estimate.Rd | 9 +- src/RcppExports.cpp | 2 +- 31 files changed, 409 insertions(+), 209 deletions(-) create mode 100644 R/lsem_aggregate_statistics.R diff --git a/DESCRIPTION b/DESCRIPTION index e203c3bc..63cdea6a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sirt Type: Package Title: Supplementary Item Response Theory Models -Version: 4.2-2 -Date: 2024-02-15 17:09:03 +Version: 4.2-13 +Date: 2024-03-01 15:48:53 Author: Alexander Robitzsch [aut,cre] () Maintainer: Alexander Robitzsch Description: diff --git a/R/RcppExports.R b/R/RcppExports.R index b55e09b2..5570cca6 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,5 +1,5 @@ ## File Name: RcppExports.R -## File Version: 4.002002 +## File Version: 4.002013 # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 diff --git a/R/lsem.MGM.stepfunctions.R b/R/lsem.MGM.stepfunctions.R index b53a31d5..ae4ce01a 100644 --- a/R/lsem.MGM.stepfunctions.R +++ b/R/lsem.MGM.stepfunctions.R @@ -1,5 +1,5 @@ ## File Name: lsem.MGM.stepfunctions.R -## File Version: 0.131 +## File Version: 0.132 #**** define object with LSEM step functions lsem.MGM.stepfunctions <- function( object, moderator.grid ) @@ -12,7 +12,7 @@ lsem.MGM.stepfunctions <- function( object, moderator.grid ) dfr <- NULL G <- length(moderator.grid) - for (gg in 1:G){ + for (gg in 1L:G){ mod.gg <- moderator.grid[gg] ind.gg <- which( ( moderator.grouped$min <=mod.gg ) & ( moderator.grouped$max > mod.gg ) ) diff --git a/R/lsem.bootstrap.R b/R/lsem.bootstrap.R index 9a47d7f1..baf37d24 100644 --- a/R/lsem.bootstrap.R +++ b/R/lsem.bootstrap.R @@ -1,5 +1,5 @@ ## File Name: lsem.bootstrap.R -## File Version: 0.431 +## File Version: 0.432 lsem.bootstrap <- function(object, R=100, verbose=TRUE, cluster=NULL, @@ -135,7 +135,7 @@ lsem.bootstrap <- function(object, R=100, verbose=TRUE, cluster=NULL, verbose=verbose, arglist=arglist) parallel::stopCluster(cl) - for (rr in 1:R){ + for (rr in 1L:R){ res_out_rr <- res_all[[rr]] parameters_boot[,rr] <- res_out_rr$parameters_boot fitstats_joint_boot[,rr] <- res_out_rr$fitstats_joint_boot diff --git a/R/lsem.estimate.R b/R/lsem.estimate.R index 943ae677..ffe114cc 100644 --- a/R/lsem.estimate.R +++ b/R/lsem.estimate.R @@ -1,5 +1,5 @@ ## File Name: lsem.estimate.R -## File Version: 1.074 +## File Version: 1.108 # estimate LSEM model lsem.estimate <- function( data, moderator, moderator.grid, @@ -23,10 +23,19 @@ lsem.estimate <- function( data, moderator, moderator.grid, } use_lavaan_survey <- FALSE + #- check if list of imputed datasets is available + is_imputed <- ! ( is.list(data) & is.data.frame(data) ) + #- data cleaning - data <- as.data.frame(data) - data <- data[ ! is.na(data[,moderator]), ] - moderator_variable <- data[,moderator] + if (!is_imputed){ + data <- as.data.frame(data) + data <- data[ ! is.na(data[,moderator]), ] + moderator_variable <- data[,moderator] + Nimp <- 0 + } else { + moderator_variable <- NULL + Nimp <- length(data) + } #- process arguments res <- lsem_estimate_proc_args( lavaan.args=lavaan.args, @@ -35,7 +44,8 @@ lsem.estimate <- function( data, moderator, moderator.grid, use_lavaan_survey=use_lavaan_survey, est_joint=est_joint, par_invariant=par_invariant, par_linear=par_linear, par_quadratic=par_quadratic, partable_joint=partable_joint, - moderator.grid=moderator.grid, se=se, verbose=verbose ) + moderator.grid=moderator.grid, se=se, verbose=verbose, + is_imputed=is_imputed) sufficient_statistics <- res$sufficient_statistics use_lavaan_survey <- res$use_lavaan_survey variables_model <- res$variables_model @@ -51,7 +61,8 @@ lsem.estimate <- function( data, moderator, moderator.grid, # group moderator if type='MGM' out <- lsem_group_moderator( data=data, type=type, moderator.grid=moderator.grid, - moderator=moderator, residualize=residualize, h=h ) + moderator=moderator, residualize=residualize, h=h, + is_imputed=is_imputed, Nimp=Nimp ) data <- out$data moderator.grouped <- out$moderator.grouped h <- out$h @@ -63,7 +74,8 @@ lsem.estimate <- function( data, moderator, moderator.grid, 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) + variables_model=variables_model, is_imputed=is_imputed, + Nimp=Nimp) G <- out$G data <- out$data weights <- out$weights @@ -91,7 +103,8 @@ lsem.estimate <- function( data, moderator, moderator.grid, variables_model=variables_model, sampling_weights=sampling_weights, has_meanstructure=has_meanstructure, sufficient_statistics=sufficient_statistics, est_joint=est_joint, - se=se, use_lavaan_survey=use_lavaan_survey, ... ) + se=se, use_lavaan_survey=use_lavaan_survey, + is_imputed=is_imputed, Nimp=Nimp, ... ) nobs <- unlist(lavfit@Data@nobs) # extract variables which are in model and data frame @@ -127,7 +140,8 @@ lsem.estimate <- function( data, moderator, moderator.grid, loc_linear_smooth=loc_linear_smooth, pd=pd, residualized_intercepts=residualized_intercepts, has_meanstructure=has_meanstructure, est_DIF=est_DIF, - residualize=residualize, ... ) + residualize=residualize, is_imputed=is_imputed, + Nimp=Nimp, moderator=moderator, ... ) dif_effects <- out2$dif_effects parameters <- out2$parameters is_meanstructure <- out2$is_meanstructure @@ -139,8 +153,11 @@ lsem.estimate <- function( data, moderator, moderator.grid, parameters_summary <- lsem_parameter_summary( parameters=parameters, moderator.density=out$moderator.density, verbose=verbose ) - out$moderator.density$Neff <- colSums(weights) - + weights0 <- weights + if (is_imputed){ + weights0 <- lsem_aggregate_statistics(x=weights) + } + out$moderator.density$Neff <- colSums(weights0) obji0 <- obji <- out$moderator.density obji$moderator <- obji$moderator obji$wgt <- obji$wgt @@ -148,10 +165,14 @@ lsem.estimate <- function( data, moderator, moderator.grid, Y <- obji0[,-1] dfr <- data.frame( M=colMeans(Y), SD=apply( Y, 2, stats::sd ), min=apply( Y, 2, min ), max=apply( Y, 2, max ) ) - x <- data[,moderator] - dfr0 <- data.frame(M=mean( x, na.rm=TRUE ), SD=out$sd.moderator, - min=min( x, na.rm=TRUE ), max=max( x, na.rm=TRUE ) ) - obji <- rbind( dfr0, dfr ) + if (is_imputed){ + x <- (data[[1]])[, moderator] + } else { + x <- data[,moderator] + } + dfr0 <- data.frame(M=mean(x, na.rm=TRUE ), SD=out$sd.moderator, + min=min(x, na.rm=TRUE ), max=max(x, na.rm=TRUE ) ) + obji <- rbind( dfr0, dfr) rownames(obji) <- NULL moderator.stat <- data.frame(variable=c('moderator','wgt', 'Neff'), obji ) @@ -188,7 +209,7 @@ lsem.estimate <- function( data, moderator, moderator.grid, partable_joint=partable_joint, dif_effects=dif_effects, sample_stats=sample_stats, loc_linear_smooth=loc_linear_smooth, - se=se, compute_se=compute_se, + se=se, compute_se=compute_se, is_imputed=is_imputed, Nimp=Nimp, class_boot=FALSE, type=type, CALL=CALL ) class(res) <- 'lsem' return(res) diff --git a/R/lsem.permutationTest.R b/R/lsem.permutationTest.R index 6013e9c9..69ac324d 100644 --- a/R/lsem.permutationTest.R +++ b/R/lsem.permutationTest.R @@ -1,5 +1,5 @@ ## File Name: lsem.permutationTest.R -## File Version: 0.592 +## File Version: 0.593 #*** permutation test for LSEM model @@ -18,13 +18,13 @@ lsem.permutationTest <- function( lsem.object, B=1000, residualize=TRUE, object <- lsem.object arglist <- list() EE <- length(entr) - for (ee in 1:EE){ + for (ee in 1L:EE){ arglist[[ entr[ee] ]] <- object[[ entr[ee] ]] } arglist2 <- lsem.object$lavaan.args NL <- length(arglist2) if (NL > 0){ - for (ll in 1:NL){ + for (ll in 1L:NL){ arglist[[ names(arglist2)[ll] ]] <- arglist2[[ names(arglist2)[ll] ]] } } @@ -140,7 +140,7 @@ lsem.permutationTest <- function( lsem.object, B=1000, residualize=TRUE, parallel::stopCluster(cl) - for (bb in 1:B){ + for (bb in 1L:B){ parameters_permutation[, bb] <- res_all[[bb]]$est parameters_summary_M[,bb] <- res_all[[bb]]$M parameters_summary_SD[,bb] <- res_all[[bb]]$SD diff --git a/R/lsem.test.R b/R/lsem.test.R index a7bbfc85..2a3d8ee6 100644 --- a/R/lsem.test.R +++ b/R/lsem.test.R @@ -1,5 +1,5 @@ ## File Name: lsem.test.R -## File Version: 0.133 +## File Version: 0.134 #**** test LSEM model based on bootstrap lsem.test <- function( mod, bmod, models=NULL ) @@ -23,7 +23,7 @@ lsem.test <- function( mod, bmod, models=NULL ) #* define design matrix A <- matrix(0, nrow=NG-1, ncol=NG) - for (gg in 1:(NG-1)){ + for (gg in 1L:(NG-1)){ A[gg,gg] <- -1 A[gg,gg+1] <- 1 } @@ -73,7 +73,7 @@ lsem.test <- function( mod, bmod, models=NULL ) if (! bmod_missing){ est_boot <- matrix(NA, nrow=NC, ncol=R) rr <- 1 - for (rr in 1:R){ + for (rr in 1L:R){ dat$y <- parameters_boot[ind_pp,rr] dat$w <- bmod$moderator_density_boot[,rr] mod12 <- stats::lm(formula=model_mm, data=dat) @@ -89,7 +89,7 @@ lsem.test <- function( mod, bmod, models=NULL ) # global Wald test for all parameters without intercept V <- stats::cov(t(est_boot)) A <- matrix(0, nrow=NC-1, ncol=NC) - for (cc in 1:(NC-1)){ + for (cc in 1L:(NC-1)){ A[cc,cc+1] <- 1 } res <- lsem_wald_test(theta=coef11, V=V, A=A) diff --git a/R/lsem_aggregate_statistics.R b/R/lsem_aggregate_statistics.R new file mode 100644 index 00000000..f973b302 --- /dev/null +++ b/R/lsem_aggregate_statistics.R @@ -0,0 +1,14 @@ +## File Name: lsem_aggregate_statistics.R +## File Version: 0.02 + + +lsem_aggregate_statistics <- function(x) +{ + Nimp <- length(x) + w <- 0 + for (ii in 1L:Nimp){ + w <- w + x[[ii]] + } + w <- w / Nimp + return(w) +} diff --git a/R/lsem_estimate_proc_args.R b/R/lsem_estimate_proc_args.R index eeb2c766..aebf0ae9 100644 --- a/R/lsem_estimate_proc_args.R +++ b/R/lsem_estimate_proc_args.R @@ -1,11 +1,15 @@ ## File Name: lsem_estimate_proc_args.R -## File Version: 0.406 +## File Version: 0.416 lsem_estimate_proc_args <- function(lavaan.args, sufficient_statistics, pseudo_weights, lavmodel, data, use_lavaan_survey, est_joint=FALSE, par_invariant=NULL, par_linear=NULL, par_quadratic=NULL, - partable_joint=NULL, se=NULL, G=NULL, moderator.grid=NULL, verbose=TRUE) + partable_joint=NULL, se=NULL, G=NULL, moderator.grid=NULL, verbose=TRUE, + is_imputed=FALSE) { + if (is_imputed){ + data <- data[[1]] + } use_pseudo_weights <- pseudo_weights > 0 if ( sufficient_statistics | use_pseudo_weights ){ @@ -54,7 +58,7 @@ lsem_estimate_proc_args <- function(lavaan.args, sufficient_statistics, data_ordered <- rep(FALSE, ncol(data1)) names(data_ordered) <- colnames(data1) NV <- ncol(data1) - for (vv in 1:NV){ + for (vv in 1L:NV){ data_ordered[vv] <- is.factor(data1[,vv]) } some_ordinal <- FALSE @@ -90,6 +94,6 @@ lsem_estimate_proc_args <- function(lavaan.args, sufficient_statistics, variables_ordered=variables_ordered, est_joint=est_joint, partable=partable, has_meanstructure=has_meanstructure, se=se, compute_se=compute_se, pseudo_weights=pseudo_weights, - some_ordinal=some_ordinal) + some_ordinal=some_ordinal, is_imputed=is_imputed) return(res) } diff --git a/R/lsem_fit_initial_model.R b/R/lsem_fit_initial_model.R index 33c53b3d..94b4a432 100644 --- a/R/lsem_fit_initial_model.R +++ b/R/lsem_fit_initial_model.R @@ -1,9 +1,9 @@ ## File Name: lsem_fit_initial_model.R -## File Version: 0.205 +## File Version: 0.212 lsem_fit_initial_model <- function(lavmodel__, lavaan_est_fun, dat, variables_model, sampling_weights, has_meanstructure, sufficient_statistics, est_joint=FALSE, - se="standard", use_lavaan_survey=FALSE, ...) + se="standard", use_lavaan_survey=FALSE, is_imputed=FALSE, Nimp=0, ...) { if (est_joint){ has_meanstructure <- TRUE @@ -12,7 +12,8 @@ lsem_fit_initial_model <- function(lavmodel__, lavaan_est_fun, dat, variables_mo #- compute sufficient statistics res <- lsem_fit_initial_model_sufficient_statistics(dat=dat, variables_model=variables_model, sampling_weights=sampling_weights, - has_meanstructure=has_meanstructure) + has_meanstructure=has_meanstructure, + is_imputed=is_imputed, Nimp=Nimp) wmean <- res$wmean wcov <- res$wcov Nobs <- res$Nobs diff --git a/R/lsem_fit_initial_model_sufficient_statistics.R b/R/lsem_fit_initial_model_sufficient_statistics.R index 7c24ed4e..990faf81 100644 --- a/R/lsem_fit_initial_model_sufficient_statistics.R +++ b/R/lsem_fit_initial_model_sufficient_statistics.R @@ -1,19 +1,37 @@ ## File Name: lsem_fit_initial_model_sufficient_statistics.R -## File Version: 0.064 +## File Version: 0.075 lsem_fit_initial_model_sufficient_statistics <- function(dat, variables_model, - sampling_weights, has_meanstructure) + sampling_weights, has_meanstructure, is_imputed=FALSE, Nimp=0) { - data_suff <- dat[, variables_model] - dat_resp <- 1 - is.na(data_suff) - wmean <- lsem_weighted_mean( x=data_suff, weights=sampling_weights, - x_resp=dat_resp)$mean - res <- lsem_weighted_cov( x=data_suff, weights=sampling_weights, x_resp=dat_resp) - wcov <- res$cov - Nobs <- round(res$Nobs) + + dat0 <- dat + if (!is_imputed){ + dat0 <- list(dat0) + } + Nimp <- max(1, Nimp) + wmean <- list() + wcov <- list() + Nobs <- list() + for (ii in 1L:Nimp){ + dat <- dat0[[ii]] + data_suff <- dat[, variables_model] + dat_resp <- 1 - is.na(data_suff) + wmean[[ii]] <- lsem_weighted_mean( x=data_suff, weights=sampling_weights, + x_resp=dat_resp)$mean + res <- lsem_weighted_cov( x=data_suff, weights=sampling_weights, x_resp=dat_resp) + wcov[[ii]] <- res$cov + Nobs[[ii]] <- round(res$Nobs) + } + + Nobs <- lsem_aggregate_statistics(x=Nobs) + wcov <- lsem_aggregate_statistics(x=wcov) + wmean <- lsem_aggregate_statistics(x=wmean) + if (! has_meanstructure){ wmean <- NULL } + #--- output res <- list(wmean=wmean, wcov=wcov, Nobs=Nobs ) return(res) diff --git a/R/lsem_fitsem.R b/R/lsem_fitsem.R index 6408d9eb..dce2806f 100644 --- a/R/lsem_fitsem.R +++ b/R/lsem_fitsem.R @@ -1,5 +1,5 @@ ## File Name: lsem_fitsem.R -## File Version: 0.643 +## File Version: 0.658 lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.grid, verbose, pars, standardized, variables_model, sufficient_statistics, @@ -8,7 +8,8 @@ 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, residualized_intercepts=NULL, residualize=TRUE, ... ) + est_DIF=FALSE, residualized_intercepts=NULL, residualize=TRUE, + is_imputed=FALSE, Nimp=0, moderator=NULL, ... ) { parameters <- NULL fits <- NULL @@ -29,13 +30,17 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr #- sufficient statistics if (sufficient_statistics){ + sample_stats <- lsem_fitsem_compute_sufficient_statistics(G=G, dat=dat, - variables_model=variables_model, weights=weights, - moderator_variable=moderator_variable, - loc_linear_smooth=loc_linear_smooth, moderator.grid=moderator.grid, - pd=pd, residualized_intercepts=residualized_intercepts, - has_meanstructure=has_meanstructure, - residualize=residualize) + variables_model=variables_model, weights=weights, + moderator_variable=moderator_variable, + loc_linear_smooth=loc_linear_smooth, + moderator.grid=moderator.grid, + pd=pd, residualized_intercepts=residualized_intercepts, + has_meanstructure=has_meanstructure, + residualize=residualize, is_imputed=is_imputed, Nimp=Nimp, + moderator=moderator) + } if (est_joint & (! sufficient_statistics)){ N <- nrow(dat) @@ -100,8 +105,13 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr } #-- separate estimation: loop over groups - for (gg in 1:G){ - dat$weight <- weights[,gg] + for (gg in 1L:G){ + + if (! is_imputed){ + dat$weight <- weights[,gg] + } else { + dat$weight <- weights[[1]][,gg] + } #***** fit the model using weighted data if (( ! sufficient_statistics) & ( est_separate)){ @@ -144,7 +154,7 @@ lsem_fitsem <- function( dat, weights, lavfit, fit_measures, NF, G, moderator.gr } else { ind <- NULL par_gg <- NULL - for (pp in 1:NP){ + for (pp in 1L:NP){ ind_pp <- which(pars==pars0[pp]) npp <- length(ind_pp) ind <- c( ind, ind_pp) diff --git a/R/lsem_fitsem_compute_sufficient_statistics.R b/R/lsem_fitsem_compute_sufficient_statistics.R index bc73b196..f4728c01 100644 --- a/R/lsem_fitsem_compute_sufficient_statistics.R +++ b/R/lsem_fitsem_compute_sufficient_statistics.R @@ -1,31 +1,56 @@ ## File Name: lsem_fitsem_compute_sufficient_statistics.R -## File Version: 0.109 +## File Version: 0.124 lsem_fitsem_compute_sufficient_statistics <- function(G, dat, variables_model, weights, moderator_variable=NULL, loc_linear_smooth=NULL, moderator.grid=NULL, pd=FALSE, residualized_intercepts=NULL, has_meanstructure=FALSE, - residualize=TRUE) + residualize=TRUE, is_imputed=FALSE, Nimp=0, moderator=NULL) { wmean <- wcov <- Nobs <- as.list(1:G) - data_suff <- dat[, variables_model] - dat_resp <- 1 - is.na(data_suff) - for (gg in 1:G){ - weights_gg <- weights[,gg] - # res <- lsem_weighted_mean( x=data_suff, weights=weights_gg, x_resp=dat_resp) - 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, - residualized_intercepts=residualized_intercepts, - has_meanstructure=has_meanstructure, residualize=residualize) - wmean[[gg]] <- res$mean - wcov[[gg]] <- res$cov - Nobs[[gg]] <- round(res$Nobs) + + dat0 <- dat + if (!is_imputed){ + dat0 <- list(dat0) } + Nimp <- max(1, Nimp) + Nobs_tt <- wcov_tt <- wmean_tt <- list() + + for (gg in 1L:G){ + + for (ii in 1L:Nimp){ + if (!is_imputed){ + weights_gg <- weights[,gg] + } else { + weights_gg <- weights[[ii]][,gg, drop=TRUE] + } + dat <- dat0[[ii]] + data_suff <- dat[, variables_model] + dat_resp <- 1 - is.na(data_suff) + if (is_imputed){ + moderator_variable <- dat[, moderator] + } + + # res <- lsem_weighted_mean( x=data_suff, weights=weights_gg, x_resp=dat_resp) + 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, + residualized_intercepts=residualized_intercepts, + has_meanstructure=has_meanstructure, residualize=residualize) + wmean_tt[[ii]] <- res$mean + wcov_tt[[ii]] <- res$cov + Nobs_tt[[ii]] <- res$Nobs + } # end ii + + wmean[[gg]] <- lsem_aggregate_statistics(x=wmean_tt) + wcov[[gg]] <- lsem_aggregate_statistics(x=wcov_tt) + Nobs[[gg]] <- round( lsem_aggregate_statistics(x=Nobs_tt) ) + + } # end gg #** adapt if mean structure is requested if ( has_meanstructure & residualize ){ - for (gg in 1:G){ + for (gg in 1L:G){ wmean[[gg]] <- residualized_intercepts[gg,] } } diff --git a/R/lsem_fitsem_joint_estimation_prepare_partable.R b/R/lsem_fitsem_joint_estimation_prepare_partable.R index 4192d8d4..266aa8e1 100644 --- a/R/lsem_fitsem_joint_estimation_prepare_partable.R +++ b/R/lsem_fitsem_joint_estimation_prepare_partable.R @@ -1,5 +1,5 @@ ## File Name: lsem_fitsem_joint_estimation_prepare_partable.R -## File Version: 0.318 +## File Version: 0.319 lsem_fitsem_joint_estimation_prepare_partable <- function(partable, G, par_invariant=NULL, par_linear=NULL, par_quadratic=NULL, pw_linear=1, @@ -12,7 +12,7 @@ lsem_fitsem_joint_estimation_prepare_partable <- function(partable, G, partable1 <- lsem_fitsem_joint_estimation_prepare_partable_include_group_label( partable=partable, gg=1, label_list=label_list) partable_mg <- partable1 - for (gg in 2:G){ + for (gg in 2L:G){ partable_gg <- partable partable_gg <- lsem_fitsem_joint_estimation_prepare_partable_include_group_label( partable=partable_gg, gg=gg, label_list=label_list) @@ -45,7 +45,7 @@ lsem_fitsem_joint_estimation_prepare_partable <- function(partable, G, if ( NI > 0 ){ partable1 <- partable_mg[1,] NV <- ncol(partable_mg) - for (vv in 1:NV){ + for (vv in 1L:NV){ if (is.numeric(partable1[1,vv])){ partable1[1,vv] <- 0 } else { @@ -58,7 +58,7 @@ lsem_fitsem_joint_estimation_prepare_partable <- function(partable, G, partable1$op <- '==' partable1c <- partable1 - for (vv in 1:NI){ + for (vv in 1L:NI){ par_vec_vv <- par_vec[vv] par_vec_names_vv <- par_vec_names[vv] ind_vv <- which( paste(partable_mg$par)==par_vec_names[vv] ) @@ -99,7 +99,7 @@ lsem_fitsem_joint_estimation_prepare_partable <- function(partable, G, par_segments[indices2, 'con_include'] <- 0 } - for (ll in 2:LV2){ + for (ll in 2L:LV2){ partable1c$con <- vv if (par_vec_vv=='inv'){ partable1c$lhs <- plabels[1] diff --git a/R/lsem_group_moderator.R b/R/lsem_group_moderator.R index eed4fadf..0e265c91 100644 --- a/R/lsem_group_moderator.R +++ b/R/lsem_group_moderator.R @@ -1,24 +1,40 @@ ## File Name: lsem_group_moderator.R -## File Version: 0.172 +## File Version: 0.178 #***** grouping a moderator variable lsem_group_moderator <- function( data, type, moderator.grid, - moderator, residualize, h ) + moderator, residualize, h, is_imputed=FALSE, Nimp=0 ) { + data0 <- data moderator.grouped <- NULL if (type=='MGM'){ G1 <- length(moderator.grid) moderator.grouped <- data.frame( min=moderator.grid[-G1], max=moderator.grid[-1] ) moderator.grouped$mid <- rowMeans( moderator.grouped) - v1 <- data[, moderator ] - v2 <- moderator.grouped$mid[1] - for (gg in 2:G1){ - v2 <- ifelse( v1 > moderator.grouped$max[gg-1], - moderator.grouped$mid[gg], v2 ) + if (! is_imputed){ + v1 <- data[, moderator ] + v2 <- moderator.grouped$mid[1] + for (gg in 2L:G1){ + v2 <- ifelse( v1 > moderator.grouped$max[gg-1], + moderator.grouped$mid[gg], v2 ) + } + data[,moderator] <- v2 + } + if (is_imputed){ + for (ii in 1L:Nimp){ + data <- data0[[ii]] + v1 <- data[, moderator ] + v2 <- moderator.grouped$mid[1] + for (gg in 2L:G1){ + v2 <- ifelse( v1 > moderator.grouped$max[gg-1], + moderator.grouped$mid[gg], v2 ) + } + data[,moderator] <- v2 + data0[[ii]] <- data + } } - data[,moderator] <- v2 # residualize <- FALSE h <- 1E-5 moderator.grid <- moderator.grouped$mid diff --git a/R/lsem_local_weights.R b/R/lsem_local_weights.R index 7694eb6f..bfbc676b 100644 --- a/R/lsem_local_weights.R +++ b/R/lsem_local_weights.R @@ -1,8 +1,9 @@ ## File Name: lsem_local_weights.R -## File Version: 0.202 +## File Version: 0.219 lsem_local_weights <- function(data.mod, moderator.grid, h, - sampling_weights=NULL, bw=NULL, kernel="gaussian") + sampling_weights=NULL, bw=NULL, kernel="gaussian", + is_imputed=FALSE, Nimp=0, data=NULL, moderator=NULL) { eps <- 1E-8 N <- length(data.mod) @@ -11,35 +12,68 @@ lsem_local_weights <- function(data.mod, moderator.grid, h, sampling_weights <- rep(1,N) no_sampling_weights <- TRUE } - sampling_weights <- sampling_weights / sum(sampling_weights) * N - # select nearest neighbor in moderator group for calculating residuals - G <- length(moderator.grid) - modgrid_index <- rep(1,N) - for (gg in 2:G){ - modgrid_index <- ifelse( abs( data.mod - moderator.grid[ modgrid_index ] ) < - abs( data.mod - moderator.grid[ gg ] ), - modgrid_index, gg ) - } - # compute weights for every grid point gg - weights <- matrix( NA, nrow=N, ncol=G ) - sd.moderator <- TAM::weighted_sd(x=data.mod, w=sampling_weights) - m.moderator <- TAM::weighted_mean(x=data.mod, w=sampling_weights) - bw_basis <- sd.moderator * N^(-1/5) - if (is.null(bw)){ - bw <- h * bw_basis - } else { - h <- bw / bw_basis - } - weights1 <- sampling_weights / sum(sampling_weights) - - moderator.density <- stats::density( data.mod, weights=weights1, - from=min(moderator.grid), to=max(moderator.grid ), n=G )$y - moderator.density <- data.frame( moderator=moderator.grid, - wgt=moderator.density / sum(moderator.density) ) - for (gg in 1:G){ - xgg <- moderator.grid[gg] - wgt <- lsem_kernel_weights(x=data.mod, x0=xgg, bw=bw, kernel=kernel) - weights[,gg] <- ifelse( wgt < eps, eps, wgt ) + bw0 <- bw + + # define list for imputed datasets + Nimp <- max(1, Nimp) + res0 <- list() + + for (ii in 1L:Nimp){ + if (is_imputed){ + data.mod <- ( data[[ii]] )[, moderator ] + } + bw <- bw0 + + sampling_weights <- sampling_weights / sum(sampling_weights) * N + # select nearest neighbor in moderator group for calculating residuals + G <- length(moderator.grid) + modgrid_index <- rep(1,N) + for (gg in 2L:G){ + modgrid_index <- ifelse( abs( data.mod - moderator.grid[ modgrid_index ] ) < + abs( data.mod - moderator.grid[ gg ] ), + modgrid_index, gg ) + } + res0$modgrid_index[[ii]] <- modgrid_index + + # compute weights for every grid point gg + weights <- matrix( NA, nrow=N, ncol=G ) + sd.moderator <- TAM::weighted_sd(x=data.mod, w=sampling_weights) + m.moderator <- TAM::weighted_mean(x=data.mod, w=sampling_weights) + bw_basis <- sd.moderator * N^(-1/5) + if (is.null(bw)){ + bw <- h * bw_basis + } else { + h <- bw / bw_basis + } + weights1 <- sampling_weights / sum(sampling_weights) + + moderator.density <- stats::density( data.mod, weights=weights1, + from=min(moderator.grid), to=max(moderator.grid ), n=G )$y + moderator.density <- data.frame( moderator=moderator.grid, + wgt=moderator.density / sum(moderator.density) ) + for (gg in 1L:G){ + xgg <- moderator.grid[gg] + wgt <- lsem_kernel_weights(x=data.mod, x0=xgg, bw=bw, kernel=kernel) + weights[,gg] <- ifelse( wgt < eps, eps, wgt ) + } + + sampling_weights <- sampling_weights / sum(sampling_weights) * N + weights <- weights * sampling_weights + + res0$weights[[ii]] <- weights + res0$moderator.density[[ii]] <- moderator.density + res0$m.moderator[[ii]] <- m.moderator + res0$sd.moderator[[ii]] <- sd.moderator + res0$bw[[ii]] <- bw + + } # end ii + if (is_imputed){ + m.moderator <- lsem_aggregate_statistics(x=res0$m.moderator) + sd.moderator <- lsem_aggregate_statistics(x=res0$sd.moderator) + moderator.density <- lsem_aggregate_statistics(x=res0$moderator.density) + bw <- lsem_aggregate_statistics(x=res0$bw) + weights <- res0$weights + modgrid_index <- res0$modgrid_index } #--- output diff --git a/R/lsem_parameter_summary.R b/R/lsem_parameter_summary.R index 3c6d48e8..0f100600 100644 --- a/R/lsem_parameter_summary.R +++ b/R/lsem_parameter_summary.R @@ -1,5 +1,5 @@ ## File Name: lsem_parameter_summary.R -## File Version: 0.235 +## File Version: 0.236 ## lsem parameter summary @@ -12,7 +12,7 @@ lsem_parameter_summary <- function( parameters, moderator.density, verbose ) utils::flush.console() } parameters_summary <- NULL - for (pp in 1:NP){ + for (pp in 1L:NP){ par.pp <- parameters[ parameters$parindex==pp, ] pars1 <- data.frame( par=paste(par.pp$par[1]), parindex=pp) x <- par.pp[,'est'] diff --git a/R/lsem_residualize.R b/R/lsem_residualize.R index e473a694..2c388247 100644 --- a/R/lsem_residualize.R +++ b/R/lsem_residualize.R @@ -1,20 +1,30 @@ ## File Name: lsem_residualize.R -## File Version: 0.433 +## File Version: 0.477 #**** residualize data lsem_residualize <- function( data, moderator, moderator.grid, lavmodel, h=1.1, bw=NULL, residualize=TRUE, sampling_weights=NULL, - eps=1E-10, verbose=TRUE, kernel="gaussian", variables_model ) + eps=1E-10, verbose=TRUE, kernel="gaussian", variables_model, + is_imputed=FALSE, Nimp=0 ) { # lavaanify model lavaanstr <- sirt_import_lavaan_lavaanify(model=lavmodel) vars <- variables_model - data.mod <- data[, moderator, drop=TRUE ] + + # values of moderator variable + if (! is_imputed){ + data.mod <- data[, moderator, drop=TRUE ] + Nimp <- 1 + } else { # use first dataset for imputed data + data.mod <- ( data[[1]] )[, moderator, drop=TRUE ] + } # compute local weights res <- lsem_local_weights(data.mod=data.mod, moderator.grid=moderator.grid, - h=h, sampling_weights=sampling_weights, bw=bw, kernel=kernel) + h=h, sampling_weights=sampling_weights, bw=bw, kernel=kernel, + is_imputed=is_imputed, Nimp=Nimp, data=data, + moderator=moderator) weights <- res$weights modgrid_index <- res$modgrid_index N <- res$N @@ -27,50 +37,90 @@ lsem_residualize <- function( data, moderator, moderator.grid, sampling_weights <- res$sampling_weights no_sampling_weights <- res$no_sampling_weights - sampling_weights <- sampling_weights / sum(sampling_weights) * N - weights <- weights * sampling_weights + res0 <- as.list(1:Nimp) + data0 <- data # residualize - dat2 <- data - V <- length(vars) - 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') - utils::flush.console() + for (ii in 1L:Nimp){ + dat2 <- data + if (is_imputed){ + dat2 <- data <- data0[[ii]] } - N <- nrow(data) - for (vv in 1:V){ - var.vv <- vars[vv] - ind_vv <- which( ! is.na( data[,var.vv] ) ) - y0 <- rep(NA,N) - for (gg in 1:G){ - x <- dat2[,moderator] - data1 <- data - data1$x <- x - res_formula <- paste0( var.vv, ' ~ x + I(x^2)' ) - weights_gg <- weights[,gg] - data1$weights_gg <- weights_gg - mod <- stats::lm( formula=res_formula, data=data1, - weights=weights_gg ) - dfr_pred <- data.frame( x=moderator.grid[gg] ) - m1 <- stats::predict( mod, dfr_pred ) - residualized_intercepts[gg,vv] <- m1 - y <- stats::resid(mod) - y0[ ind_vv ] <- y - dat2[, var.vv ] <- ifelse( modgrid_index==gg, y0, dat2[, var.vv ] ) + + V <- length(vars) + residualized_intercepts <- matrix( 0, nrow=G, ncol=V) + colnames(residualized_intercepts) <- vars + rownames(residualized_intercepts) <- round( moderator.grid, 3 ) + if (residualize){ + if (verbose & ii==1){ + cat('** Residualize Data\n') + utils::flush.console() + } + N <- nrow(data) + for (vv in 1L:V){ + var.vv <- vars[vv] + ind_vv <- which( ! is.na( data[,var.vv] ) ) + y0 <- rep(NA,N) + for (gg in 1L:G){ + x <- dat2[,moderator] + data1 <- data + data1$x <- x + res_formula <- paste0( var.vv, ' ~ x + I(x^2)' ) + if (!is_imputed){ + weights_gg <- weights[,gg] + } else { + weights_gg <- (weights[[ii]])[,gg, drop=TRUE] + } + data1$weights_gg <- weights_gg + mod <- stats::lm( formula=res_formula, data=data1, + weights=weights_gg ) + dfr_pred <- data.frame( x=moderator.grid[gg] ) + m1 <- stats::predict( mod, dfr_pred ) + residualized_intercepts[gg,vv] <- m1 + y <- stats::resid(mod) + y0[ ind_vv ] <- y + if (!is_imputed){ + modgrid_index1 <- modgrid_index + } else { + modgrid_index1 <- modgrid_index[[ii]] + } + + dat2[, var.vv] <- ifelse( modgrid_index1==gg, y0, dat2[, var.vv] ) + } + } # end vv + } # end residualize=TRUE + res <- list( resid_vars=vars, data=dat2, weights_grid=weights, bw=bw, h=h, + moderator.density=moderator.density, sd.moderator=sd.moderator, G=G, N=N, + residualized_intercepts=residualized_intercepts, + sampling_weights=sampling_weights, + no_sampling_weights=no_sampling_weights, + m.moderator=m.moderator, residualize=residualize, + is_imputed=is_imputed, Nimp=Nimp) + res0[[ii]] <- res + + } # end imputation loop ii + + #--- process output + if (! is_imputed){ + res <- res0[[1]] + } else { + res <- res0[[1]] + entries <- c('data','residualized_intercepts') + for (ee in entries){ + v1 <- list() + for (ii in 1L:Nimp){ + v1[[ii]] <- res0[[ii]][[ee]] } + res[[ee]] <- v1 } - } - #--- OUTPUT - res <- list( resid_vars=vars, data=dat2, weights_grid=weights, bw=bw, h=h, - moderator.density=moderator.density, sd.moderator=sd.moderator, G=G, N=N, - residualized_intercepts=residualized_intercepts, - sampling_weights=sampling_weights, no_sampling_weights=no_sampling_weights, - m.moderator=m.moderator, residualize=residualize ) + #- aggregate residualized intercepts + ee <- 'residualized_intercepts' + res[[ee]] <- lsem_aggregate_statistics(x=res[[ee]]) + + } # end process output + + #-- out return(res) } diff --git a/R/m_est.R b/R/m_est.R index 76456498..55fa28a1 100644 --- a/R/m_est.R +++ b/R/m_est.R @@ -1,5 +1,5 @@ ## File Name: m_est.R -## File Version: 0.150 +## File Version: 0.151 m_est <- function(data, par, optfun_case=NULL, gradfun_case=NULL, bread=TRUE, optimizer="optim", @@ -25,7 +25,7 @@ m_est <- function(data, par, optfun_case=NULL, gradfun_case=NULL, grad_case <- matrix(NA, nrow=N, ncol=NP) colnames(grad_case) <- names(par) args1 <- args - for (pp in 1:NP){ + for (pp in 1L:NP){ par1 <- m_est_add_increment(x=par, pos=pp, h=h) args1$par <- par1 val2 <- do.call(what=what_optfun, args=args1) @@ -84,14 +84,14 @@ m_est <- function(data, par, optfun_case=NULL, gradfun_case=NULL, hess <- array(0, dim=c(N,NP,NP)) h <- 1e-4 pp <- 1 - for (pp in 1:NP){ + for (pp in 1L:NP){ par_pp <- m_est_add_increment(x=coef, pos=pp, h=h) scores2 <- -gradfun_case(data=data, par=par_pp) hess[,pp,] <- (scores2-scores)/h } B <- matrix(0, nrow=NP, ncol=NP) - for (pp in 1:NP){ - for (hh in 1:NP){ + for (pp in 1L:NP){ + for (hh in 1L:NP){ B[pp,hh] <- mean( hess[,pp,hh] ) } } diff --git a/R/rasch.pairwise.itemcluster.R b/R/rasch.pairwise.itemcluster.R index ed674b4d..8b625cb5 100644 --- a/R/rasch.pairwise.itemcluster.R +++ b/R/rasch.pairwise.itemcluster.R @@ -1,5 +1,5 @@ ## File Name: rasch.pairwise.itemcluster.R -## File Version: 0.32 +## File Version: 0.35 #***** Pairwise estimation with itemclusters @@ -27,7 +27,7 @@ rasch.pairwise.itemcluster <- function( dat, itemcluster=NULL, zerosum <- FALSE } # create count tables - weights <- weights / sum(weights) * nrow(dat) +# weights <- weights / sum(weights) * nrow(dat) sw <- sqrt(weights) dat0 <- sw*(dat==0) dat1 <- sw*(dat==1) @@ -45,15 +45,16 @@ rasch.pairwise.itemcluster <- function( dat, itemcluster=NULL, eps0 <- eps <- exp(b) max.change <- 10 iter <- 1 - + tol <- 1e-10 #**** start algorithm while( max.change > conv ){ b0 <- b eps0 <- eps m1 <- sirt_matrix2( eps0, nrow=I) + matrix( eps0, nrow=I, ncol=I ) + m1 <- m1+tol g1 <- rowSums(nij/m1) - eps <- Aij_rowsums/g1 - b <- log(eps) + eps <- Aij_rowsums/(g1+tol) + b <- log(eps+tol) # include item parameter constraints if ( ! is.null(b.fixed) ){ eps[ b.fixed[,1] ] <- b.fixed[,3] diff --git a/R/summary.lsem.R b/R/summary.lsem.R index f1a4e615..7ddc83de 100644 --- a/R/summary.lsem.R +++ b/R/summary.lsem.R @@ -1,5 +1,5 @@ ## File Name: summary.lsem.R -## File Version: 0.412 +## File Version: 0.415 #-- summary lsem @@ -38,6 +38,9 @@ summary.lsem <- function( object, file=NULL, digits=3, ... ) cat( paste0( 'Used observations in analysis', sp_eq, round(object$nobs, digits) ), '\n') cat('Used sampling weights:', ! object$no_sampling_weights, '\n') + if (object$is_imputed){ + cat( paste0( 'Number of imputed datasets: ', object$Nimp, '\n') ) + } if ( object$type=='LSEM'){ cat( paste0( 'Bandwidth factor', sp_eq, round(object$h,digits) ), '\n') cat( paste0( 'Bandwidth', sp_eq, round(object$bw,digits) ), '\n') diff --git a/R/xxirt.R b/R/xxirt.R index 75b68399..ac8f9b2b 100644 --- a/R/xxirt.R +++ b/R/xxirt.R @@ -1,5 +1,5 @@ ## File Name: xxirt.R -## File Version: 1.111 +## File Version: 1.112 #--- user specified item response model @@ -51,7 +51,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, partable <- xxirt_createParTable( dat=dat, itemtype=itemtype, customItems=customItems ) } - + # process partable and itemtype res <- xxirt_proc_ParTable( itemtype=itemtype, partable=partable, items=items ) itemtype <- res$itemtype @@ -66,7 +66,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, # create item list item_list <- xxirt_createItemList( customItems=customItems, itemtype=itemtype, items=items, partable=partable ) - + # shortcut for calculating expected counts dat1_resp <- xxirt_prepare_response_data(G=G, group_index=group_index, weights=weights, dat1=dat1, dat_resp=dat_resp, maxK=maxK ) @@ -74,12 +74,12 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, #*** starting values item parameters par0 <- xxirt_partable_extract_freeParameters( partable=partable ) par1 <- xxirt_ThetaDistribution_extract_freeParameters( customTheta=customTheta ) - + #*** if (is.null(customTheta$some_bound)){ customTheta$some_bound <- FALSE } - + #*** verbose verbose1 <- verbose==1 verbose2 <- verbose==2 @@ -96,7 +96,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, } do_cv <- cv_kfold>0 em_count <- 1 - + while(em_iterate){ #--- create list with arguments for EM algorithm @@ -118,7 +118,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, em_args$verbose2 <- FALSE em_args$verbose3 <- FALSE } - + #--- run EM algorithm em_out <- res <- do.call(what=xxirt_em_algorithm, args=em_args) diff --git a/R/xxirt_em_algorithm.R b/R/xxirt_em_algorithm.R index 6bdfe217..7fdea80e 100644 --- a/R/xxirt_em_algorithm.R +++ b/R/xxirt_em_algorithm.R @@ -1,5 +1,5 @@ ## File Name: xxirt_em_algorithm.R -## File Version: 0.087 +## File Version: 0.088 xxirt_em_algorithm <- function(maxit, verbose1, verbose2, verbose3, disp, item_list, items, Theta, ncat, partable, partable_index, dat, resp_index, @@ -25,7 +25,7 @@ xxirt_em_algorithm <- function(maxit, verbose1, verbose2, verbose3, disp, item_l probs_items <- xxirt_compute_itemprobs( item_list=item_list, items=items, Theta=Theta, ncat=ncat, partable=partable, partable_index=partable_index ) - + #*** compute individual likelihood p.xi.aj <- xxirt_compute_likelihood( probs_items=probs_items, dat=dat, resp_index=resp_index, dat_resp_bool=dat_resp_bool ) @@ -33,7 +33,7 @@ xxirt_em_algorithm <- function(maxit, verbose1, verbose2, verbose3, disp, item_l #*** compute prior distribution 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, diff --git a/R/xxirt_newton_raphson.R b/R/xxirt_newton_raphson.R index f5b5ba1a..e839b5c4 100644 --- a/R/xxirt_newton_raphson.R +++ b/R/xxirt_newton_raphson.R @@ -1,5 +1,5 @@ ## File Name: xxirt_newton_raphson.R -## File Version: 0.215 +## File Version: 0.216 xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, @@ -27,7 +27,7 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, NPI <- length(par_items) NPT <- length(par_Theta) NP <- NPI+NPT - + em_args$parindex_items <- 1:NPI if (NPI==0){ em_args$parindex_items <- NULL @@ -89,11 +89,11 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, upper <- c( partable_free$upper, customTheta$upper ) } else { lower <- c( customTheta$lower ) - upper <- c( customTheta$upper ) + upper <- c( customTheta$upper ) } if (!is.null(lower)){ names(lower) <- names(x) } if (!is.null(upper)){ names(upper) <- names(x) } - + if (verbose){ cat( paste0('****** Newton-Raphson Optimization ********\n')) utils::flush.console() diff --git a/R/xxirt_nr_grad_fun_Rcpp.R b/R/xxirt_nr_grad_fun_Rcpp.R index cde5c2c8..37b3f54d 100644 --- a/R/xxirt_nr_grad_fun_Rcpp.R +++ b/R/xxirt_nr_grad_fun_Rcpp.R @@ -1,5 +1,5 @@ ## File Name: xxirt_nr_grad_fun_Rcpp.R -## File Version: 0.173 +## File Version: 0.174 xxirt_nr_grad_fun_Rcpp <- function(x, em_args, eps=1e-100) { @@ -11,7 +11,7 @@ xxirt_nr_grad_fun_Rcpp <- function(x, em_args, eps=1e-100) eps2 <- 1e-300 grad <- 0*x - + #*** compute prior distribution prior_Theta0 <- xxirt_compute_prior_Theta_from_x(x=x, em_args=em_args) diff --git a/R/xxirt_nr_optim_fun.R b/R/xxirt_nr_optim_fun.R index 6b4ab6c6..75e02d2d 100644 --- a/R/xxirt_nr_optim_fun.R +++ b/R/xxirt_nr_optim_fun.R @@ -1,31 +1,31 @@ ## File Name: xxirt_nr_optim_fun.R -## File Version: 0.123 +## File Version: 0.124 xxirt_nr_optim_fun <- function(x, em_args, output_all=FALSE) { pen_val <- 0 - + #*** compute prior distribution prior_Theta <- xxirt_compute_prior_Theta_from_x(x=x, em_args=em_args) #* compute item response probabilities probs_items <- xxirt_compute_prob_item_from_x(x=x, em_args=em_args) - + #*** individual likelihood p.xi.aj <- xxirt_compute_likelihood( probs_items=probs_items, dat=em_args$dat, dat_resp_bool=em_args$dat_resp_bool ) - + #*** compute likelihood function ll_case <- xxirt_compute_casewise_likelihood(prior_Theta=prior_Theta, group=em_args$group, p.xi.aj=p.xi.aj) ll0 <- ll <- -sum( em_args$weights*log(ll_case) ) - + #*** add prior distributions x1 <- x[ em_args$parindex_items ] partable1 <- xxirt_partable_include_freeParameters( partable=em_args$partable, x=x1) par_prior <- sum( xxirt_mstep_itemParameters_evalPrior(partable1, h=0) ) ll <- ll+par_prior - + #*** add penalty function if (!is.null(em_args$penalty_fun_item)){ pen_val <- em_args$penalty_fun_item(x=x1) @@ -34,7 +34,7 @@ xxirt_nr_optim_fun <- function(x, em_args, output_all=FALSE) res <- ll if (output_all){ res <- list(opt_fun=ll, pen_val=pen_val, par_prior=par_prior, ll=ll0) - } + } return(res) } diff --git a/R/xxirt_proc_ParTable.R b/R/xxirt_proc_ParTable.R index 15946c52..dd3a0e28 100644 --- a/R/xxirt_proc_ParTable.R +++ b/R/xxirt_proc_ParTable.R @@ -1,5 +1,5 @@ ## File Name: xxirt_proc_ParTable.R -## File Version: 0.463 +## File Version: 0.464 #-- process parameter table xxirt_proc_ParTable <- function( itemtype, partable, items ) @@ -22,7 +22,7 @@ xxirt_proc_ParTable <- function( itemtype, partable, items ) partable[ ! partable$est, 'prior' ] <- NA #*** list with parameter table indices partable$parfree <- 1*partable$est - partable_index <- as.list( 1:I ) + partable_index <- as.list( 1:I ) for (ii in 1:I){ partable_index[[ii]] <- which( partable$itemnr==ii ) } diff --git a/README.md b/README.md index 5600b0aa..3a9beaf7 100644 --- a/README.md +++ b/README.md @@ -22,9 +22,9 @@ The CRAN version can be installed from within R using: utils::install.packages("sirt") ``` -#### GitHub version `sirt` 4.2-2 (2024-02-15) +#### GitHub version `sirt` 4.2-13 (2024-03-01) -[![](https://img.shields.io/badge/github%20version-4.2--2-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   +[![](https://img.shields.io/badge/github%20version-4.2--13-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/inst/NEWS b/inst/NEWS index 01f13183..fd775dc4 100644 --- a/inst/NEWS +++ b/inst/NEWS @@ -37,13 +37,15 @@ CHANGELOG sirt ------------------------------------------------------------------------ -VERSIONS sirt 4.2 | 2024-02-15 | Last: sirt 4.2-2 +VERSIONS sirt 4.2 | 2024-03-01 | Last: sirt 4.2-13 ------------------------------------------------------------------------ FIXED * fixed a bug in handling sufficient statistics in lsem.permutationTest() (thanks to @martinabader, #22) - - +FIXED * fixed numerical instabilities in rasch.pairwise.itemcluster() + in small samples +ADDED * added option for using a list of imputed datasets in + lsem.estimate() DATA * included/modified datasets: --- EXAMP * included/modified examples: --- diff --git a/man/lsem.estimate.Rd b/man/lsem.estimate.Rd index ef2693aa..7a09bbe8 100644 --- a/man/lsem.estimate.Rd +++ b/man/lsem.estimate.Rd @@ -1,5 +1,5 @@ %% File Name: lsem.estimate.Rd -%% File Version: 0.769 +%% File Version: 0.773 \name{lsem.estimate} \alias{lsem.estimate} @@ -71,7 +71,7 @@ lsem.bootstrap(object, R=100, verbose=TRUE, cluster=NULL, %- maybe also 'usage' for other objects documented here. \arguments{ \item{data}{ -Data frame +Data frame or a list of imputed datasets } \item{moderator}{ Variable name of the moderator @@ -242,8 +242,9 @@ competing factor analytic approaches for the investigation of measurement invari Alexander Robitzsch, Oliver Luedtke, Andrea Hildebrandt } -%% \note{ -%% } +% \note{ +% } + %% Note that the specified model in \code{lavaan} syntax should be explicitly saved as %% the object \code{lavmodel} in the global envirnoment as it is done in the examples. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d68223e3..7110c44d 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,5 +1,5 @@ //// File Name: RcppExports.cpp -//// File Version: 4.002002 +//// File Version: 4.002013 // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393