diff --git a/DESCRIPTION b/DESCRIPTION index d2f10b4..913bdc9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sirt Type: Package Title: Supplementary Item Response Theory Models -Version: 4.2-40 -Date: 2024-03-19 13:56:57 +Version: 4.2-51 +Date: 2024-04-12 16:00:09 Author: Alexander Robitzsch [aut,cre] () Maintainer: Alexander Robitzsch Description: diff --git a/NAMESPACE b/NAMESPACE index 1c48b64..70a3ed7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -398,6 +398,7 @@ export(xxirt_createParTable) export(xxirt_createThetaDistribution) export(xxirt_hessian) export(xxirt_modifyParTable) +export(xxirt_sandwich_pml) export(yen.q3) # defunct diff --git a/R/IRT.anova.sirt.R b/R/IRT.anova.sirt.R new file mode 100644 index 0000000..db10da5 --- /dev/null +++ b/R/IRT.anova.sirt.R @@ -0,0 +1,41 @@ +## File Name: IRT.anova.sirt.R +## File Version: 0.02 + + + +IRT.anova.sirt <- function (object, ...) +{ + cl <- match.call() + cl2 <- paste(cl)[-1] + if (length(list(object, ...)) !=2) { + stop('anova method can only be applied for comparison of two models.\n') + } + objects <- list(object, ...) + model1a <- objects[[1]] + model2a <- objects[[2]] + model1 <- IRT.IC(model1a) + model2 <- IRT.IC(model2a) + dfr1 <- data.frame(Model=cl2[1], loglike=model1['loglike'], + Deviance=-2 * model1['loglike']) + dfr1$Npars <- model1['Npars'] + dfr1$AIC <- model1['AIC'] + dfr1$BIC <- model1['BIC'] + dfr2 <- data.frame(Model=cl2[2], loglike=model2['loglike'], + Deviance=-2 * model2['loglike']) + dfr2$Npars <- model2['Npars'] + dfr2$AIC <- model2['AIC'] + dfr2$BIC <- model2['BIC'] + dfr <- rbind(dfr1, dfr2) + dfr <- dfr[ order(dfr$Npars), ] + dfr$Chisq <- NA + dfr$df <- NA + dfr$p <- NA + dfr[1, 'Chisq'] <- dfr[1, 'Deviance'] - dfr[2, 'Deviance'] + dfr[1, 'df'] <- abs(dfr[1, 'Npars'] - dfr[2, 'Npars']) + dfr[1, 'p'] <- round(1 - stats::pchisq(dfr[1, 'Chisq'], df=dfr[1,'df']), 5) + for (vv in 2L:(ncol(dfr))){ + dfr[, vv] <- round(dfr[, vv], 5) + } + rownames(dfr) <- NULL + invisible(dfr) +} diff --git a/R/RcppExports.R b/R/RcppExports.R index 1767ae4..4edca54 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,5 +1,5 @@ ## File Name: RcppExports.R -## File Version: 4.002040 +## File Version: 4.002051 # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 @@ -407,3 +407,15 @@ sirt_rcpp_xxirt_newton_raphson_derivative_par <- function(dat, dat_resp_bool, ra .Call('_sirt_sirt_rcpp_xxirt_newton_raphson_derivative_par', PACKAGE='sirt', dat, dat_resp_bool, ratio, p_xi_aj, item, prior_Theta, group0, weights, ll_case0, eps) } +sirt_rcpp_xxirt_nr_pml_opt_fun <- function(prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps) { + .Call('_sirt_sirt_rcpp_xxirt_nr_pml_opt_fun', PACKAGE='sirt', prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps) +} + +sirt_rcpp_xxirt_nr_pml_grad_fun_eval <- function(prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps, NP, der_prior_Theta, val1, val2, pp_Theta, der_probs_items, index_freq1, index_freq2) { + .Call('_sirt_sirt_rcpp_xxirt_nr_pml_grad_fun_eval', PACKAGE='sirt', prior_Theta, probs_items, freq1, freq2, W1, W2_long, G, K, I, TP, NI2, eps, NP, der_prior_Theta, val1, val2, pp_Theta, der_probs_items, index_freq1, index_freq2) +} + +sirt_rcpp_xxirt_nr_pml_casewise_opt_fun <- function(prior_Theta, probs_items, W1, W2_long, G, K, I, TP, NI2, eps, group0, weights, dat1, dat_resp) { + .Call('_sirt_sirt_rcpp_xxirt_nr_pml_casewise_opt_fun', PACKAGE='sirt', prior_Theta, probs_items, W1, W2_long, G, K, I, TP, NI2, eps, group0, weights, dat1, dat_resp) +} + diff --git a/R/anova_sirt.R b/R/anova_sirt.R index 2dd392e..419db49 100644 --- a/R/anova_sirt.R +++ b/R/anova_sirt.R @@ -1,5 +1,5 @@ ## File Name: anova_sirt.R -## File Version: 0.242 +## File Version: 0.244 ############################################################## # anova rasch.mml @@ -39,44 +39,6 @@ anova.xxirt <- anova.rasch.mml ############################################################## -################################################################## -IRT.anova.sirt <- function (object, ...) -{ - cl <- match.call() - cl2 <- paste(cl)[-1] - if (length(list(object, ...)) !=2) { - stop('anova method can only be applied for comparison of two models.\n') - } - objects <- list(object, ...) - model1a <- objects[[1]] - model2a <- objects[[2]] - model1 <- IRT.IC(model1a) - model2 <- IRT.IC(model2a) - dfr1 <- data.frame(Model=cl2[1], loglike=model1['loglike'], - Deviance=-2 * model1['loglike']) - dfr1$Npars <- model1['Npars'] - dfr1$AIC <- model1['AIC'] - dfr1$BIC <- model1['BIC'] - dfr2 <- data.frame(Model=cl2[2], loglike=model2['loglike'], - Deviance=-2 * model2['loglike']) - dfr2$Npars <- model2['Npars'] - dfr2$AIC <- model2['AIC'] - dfr2$BIC <- model2['BIC'] - dfr <- rbind(dfr1, dfr2) - dfr <- dfr[ order(dfr$Npars), ] - dfr$Chisq <- NA - dfr$df <- NA - dfr$p <- NA - dfr[1, 'Chisq'] <- dfr[1, 'Deviance'] - dfr[2, 'Deviance'] - dfr[1, 'df'] <- abs(dfr[1, 'Npars'] - dfr[2, 'Npars']) - dfr[1, 'p'] <- round(1 - stats::pchisq(dfr[1, 'Chisq'], df=dfr[1,'df']), 5) - for (vv in 2:(ncol(dfr))){ - dfr[, vv] <- round(dfr[, vv], 5) - } - rownames(dfr) <- NULL - invisible(dfr) -} - ############################################################## # Likelihood ratio test for rasch.copula2 objects anova.rasch.copula2 <- function( object, ... ) @@ -116,7 +78,7 @@ anova.rasch.copula2 <- function( object, ... ) dfr[1,'Chisq'] <- dfr[1,'Deviance'] - dfr[2,'Deviance'] dfr[1,'df'] <- abs( dfr[1,'Npars'] - dfr[2,'Npars'] ) dfr[1, 'p' ] <- round( 1 - stats::pchisq( dfr[1,'Chisq'], df=dfr[1,'df'] ), 5 ) - for ( vv in 2:( ncol(dfr))){ + for ( vv in 2L:( ncol(dfr))){ dfr[,vv] <- round( dfr[,vv], 5 ) } print(dfr) diff --git a/R/btm.R b/R/btm.R index 2e838fe..8f19f49 100644 --- a/R/btm.R +++ b/R/btm.R @@ -1,5 +1,5 @@ ## File Name: btm.R -## File Version: 1.537 +## File Version: 1.539 #--- Bradley-Terry model in sirt @@ -20,7 +20,7 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU } if ( ignore.ties ){ - admiss <- admiss[1:2] + admiss <- admiss[c(1,2)] delta <- -99 est.delta <- FALSE } @@ -252,7 +252,7 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU # fit statistics res0 <- btm_fit_statistics( probs=probs, dat0=dat0, ind1=ind1, ind2=ind2, - TP=TP, judge=judge, wgt.ties=wgt.ties ) + TP=TP, judge=judge, wgt.ties=wgt.ties ) effects$outfit <- res0$outfit effects$infit <- res0$infit multiple_judges <- res0$multiple_judges @@ -267,10 +267,10 @@ btm <- function( data, judge=NULL, ignore.ties=FALSE, fix.eta=NULL, fix.delta=NU #--- output list effects <- effects[ order(effects$propscore, decreasing=TRUE), ] res <- list( effects=effects, pars=pars, summary.effects=summary.effects, - mle.rel=mle.rel, sepG=sep.rel, probs=probs, data=dat0, - multiple_judges=multiple_judges, fit_judges=fit_judges, - residuals=residuals, eps=eps, ignore.ties=ignore.ties, - wgt.ties=wgt.ties, time_alg=time_alg, ll=ll, dat=dat) + mle.rel=mle.rel, sepG=sep.rel, probs=probs, data=dat0, + multiple_judges=multiple_judges, fit_judges=fit_judges, + residuals=residuals, eps=eps, ignore.ties=ignore.ties, + wgt.ties=wgt.ties, time_alg=time_alg, ll=ll, dat=dat) res$CALL <- CALL res$iter <- iter ic <- list( n=length(teams), D=nrow(dat0) ) diff --git a/R/ccov.np.R b/R/ccov.np.R index c5a916b..f069329 100644 --- a/R/ccov.np.R +++ b/R/ccov.np.R @@ -1,5 +1,5 @@ ## File Name: ccov.np.R -## File Version: 1.221 +## File Version: 1.225 #---- nonparametric estimation of conditional covariance @@ -36,12 +36,12 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200), } icc_items <- matrix( 0, length(thetagrid), I ) if ( I >=20 ){ - display <- seq( 1, I, floor( I/20 ) )[ 2:20 ] + display <- seq( 1, I, floor( I/20 ) )[ 2L:20 ] } else { display <- 20 } i <- 1 - for ( ii in 1:I ){ + for ( ii in 1L:I ){ obs_ii <- ! is.na( data[,ii] ) x <- score[ obs_ii ] y <- data[ obs_ii, ii ] @@ -62,10 +62,11 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200), utils::flush.console() } # calculation of conditional covariance - ccov.table <- data.frame( 'item1ID'=rep( 1:I, I), 'item2ID'=rep( 1:I, each=I ) ) + ccov.table <- data.frame( item1ID=rep(1L:I, I), + item2ID=rep( 1L:I, each=I) ) ccov.table <- ccov.table[ ccov.table$item1ID < ccov.table$item2ID, ] ccov.table$N <- apply( ccov.table, 1, FUN=function(ll){ - sum( rowSums( is.na( data[, c( ll[1], ll[2] ) ] ) )==0 ) } ) + sum( rowSums( is.na( data[, c( ll[1], ll[2] ) ] ) )==0 ) } ) ccov.table <- ccov.table[ ccov.table$N > 0, ] ccov.table$item1 <- colnames(data)[ ccov.table$item1ID ] ccov.table$item2 <- colnames(data)[ ccov.table$item2ID ] @@ -76,9 +77,9 @@ ccov.np <- function( data, score, bwscale=1.1, thetagrid=seq( -3,3,len=200), ccov.matrix <- prod.matrix <- matrix( 0, nrow=length(thetagrid), ncol=FF ) ii <- 1 ccov_sum_score <- rep(NA, FF) - for (ff in 1:FF){ + for (ff in 1L:FF){ if (FF>20){ - display <- seq( 1, FF, floor( FF/20 ) )[ 2:20 ] + display <- seq( 1, FF, floor( FF/20 ) )[ 2L:20 ] } else { display <- seq(1,FF) } diff --git a/R/ccov_np_compute_ccov_sum_score.R b/R/ccov_np_compute_ccov_sum_score.R index 4b9e6a2..1e874c2 100644 --- a/R/ccov_np_compute_ccov_sum_score.R +++ b/R/ccov_np_compute_ccov_sum_score.R @@ -1,5 +1,5 @@ ## File Name: ccov_np_compute_ccov_sum_score.R -## File Version: 0.161 +## File Version: 0.162 ccov_np_compute_ccov_sum_score <- function(score, data, use_rcpp=TRUE) { @@ -8,7 +8,7 @@ ccov_np_compute_ccov_sum_score <- function(score, data, use_rcpp=TRUE) NS <- length(scores) ccov_ff <- rep(NA,NS) if (!use_rcpp){ - for (ss in 1:NS){ + for (ss in 1L:NS){ i1 <- which(score==scores[ss]) s1 <- stats::cov.wt(x=data[i1,], method='ML') ccov_ff[ss] <- s1$cov[1,2] diff --git a/R/decategorize.R b/R/decategorize.R index 477a8c1..827eaa7 100644 --- a/R/decategorize.R +++ b/R/decategorize.R @@ -1,5 +1,5 @@ ## File Name: decategorize.R -## File Version: 0.122 +## File Version: 0.124 #* decategorize decategorize <- function( dat, categ_design=NULL ) @@ -10,9 +10,9 @@ decategorize <- function( dat, categ_design=NULL ) #** handle categories if ( ! is.null( dfr ) ){ - vars <- sort( unique( paste( dfr$variable ))) + vars <- sort(unique(paste(dfr$variable))) VV <- length(vars) - for (vv in 1:VV){ + for (vv in 1L:VV){ dfr.vv <- dfr[ paste(dfr$variable)==vars[vv], ] dat4[, vars[vv] ] <- dfr.vv[ match( dat3[,vars[vv]], dfr.vv$recode ), 'orig'] } diff --git a/R/dif.strata.variance.R b/R/dif.strata.variance.R index 8630f04..90420db 100644 --- a/R/dif.strata.variance.R +++ b/R/dif.strata.variance.R @@ -1,5 +1,5 @@ ## File Name: dif.strata.variance.R -## File Version: 0.151 +## File Version: 0.153 @@ -15,7 +15,7 @@ dif.strata.variance <- function( dif, se.dif, itemcluster ) stratadif$M.DIF <- stats::aggregate( dif, list(itemcluster), mean, na.rm=TRUE )[,2] # DIF in strata SS <- nrow(stratadif) - for (ss in 1:SS){ + for (ss in 1L:SS){ items.ss <- which( itemcluster==stratadif[ss,'strata'] ) dif.ss <- dif[ items.ss ] difv.ss <- dif.variance( dif=dif.ss, se.dif=se.dif[ items.ss ] ) @@ -25,8 +25,9 @@ dif.strata.variance <- function( dif, se.dif, itemcluster ) stratadif[ is.na(stratadif ) ] <- 0 sd_ni1 <- stratadif$N.Items-1 - weighted.DIFSD <- sum(stratadif$N.Items/sum(stratadif$N.Items)*stratadif$weighted.tau) - unweighted.DIFSD <- sum( sd_ni1/ (sum(stratadif$N.Items)-1)*stratadif$unweighted.tau) + NI <- sum(stratadif$N.Items) + weighted.DIFSD <- sum(stratadif$N.Items/NI*stratadif$weighted.tau) + unweighted.DIFSD <- sum( sd_ni1/(NI-1)*stratadif$unweighted.tau) #-- output res <- list( stratadif=stratadif, weighted.DIFSD=weighted.DIFSD, diff --git a/R/mcmc.2pno.R b/R/mcmc.2pno.R index fbd88cf..74c298d 100644 --- a/R/mcmc.2pno.R +++ b/R/mcmc.2pno.R @@ -1,5 +1,5 @@ ## File Name: mcmc.2pno.R -## File Version: 1.26 +## File Version: 1.271 ############################################## # MCMC estimation 2PNO model mcmc.2pno <- function(dat, weights=NULL, burnin=500, iter=1000, N.sampvalues=1000, @@ -24,7 +24,7 @@ mcmc.2pno <- function(dat, weights=NULL, burnin=500, iter=1000, N.sampvalues=100 # item parameters in matrix form aM <- matrix( a, nrow=N, ncol=I, byrow=TRUE) bM <- matrix( b, nrow=N, ncol=I, byrow=TRUE) - theta <- qnorm( ( rowMeans( dat0,na.rm=TRUE ) + .01 ) / 1.02 ) + theta <- stats::qnorm( ( rowMeans( dat0,na.rm=TRUE ) + .01 ) / 1.02 ) # define lower and upper thresholds ZZ <- 1000 threshlow <- -ZZ + ZZ*dat diff --git a/R/mgsem.R b/R/mgsem.R index 031a0e6..d3358c7 100644 --- a/R/mgsem.R +++ b/R/mgsem.R @@ -1,5 +1,5 @@ ## File Name: mgsem.R -## File Version: 0.553 +## File Version: 0.554 mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL, estimator="ML", p_me=2, p_pen=1, pen_type="scad", @@ -35,7 +35,7 @@ mgsem <- function(suffstat, model, data=NULL, group=NULL, weights=NULL, groups <- names(suffstat) G <- length(groups) if (is.null(groups)){ - groups <- paste0('Group',1:G) + groups <- paste0('Group',1L:G) } data_proc <- NULL is_data <- FALSE diff --git a/R/mgsem_bdiag.R b/R/mgsem_bdiag.R index a5621a8..0f6a5b5 100644 --- a/R/mgsem_bdiag.R +++ b/R/mgsem_bdiag.R @@ -1,5 +1,5 @@ ## File Name: mgsem_bdiag.R -## File Version: 0.05 +## File Version: 0.06 mgsem_bdiag <- function(x1, x2) { @@ -8,7 +8,7 @@ mgsem_bdiag <- function(x1, x2) n2 <- ncol(x2) mat <- matrix(0,nrow=n1+n2,ncol=n1+n2) rownames(mat) <- colnames(mat) - mat[1:n1,1:n1] <- x1 - mat[n1+1:n2,n1+1:n2] <- x2 + mat[1L:n1,1L:n1] <- x1 + mat[n1+1L:n2,n1+1L:n2] <- x2 return(mat) } diff --git a/R/mgsem_cd_opt.R b/R/mgsem_cd_opt.R index 2092b89..d2bb796 100644 --- a/R/mgsem_cd_opt.R +++ b/R/mgsem_cd_opt.R @@ -1,5 +1,5 @@ ## File Name: mgsem_cd_opt.R -## File Version: 0.177 +## File Version: 0.178 mgsem_cd_opt <- function(x, opt_fun_args, tol=1e-4, eps_approx=1e-20, maxiter=100, h=1e-4, verbose=TRUE, interval_length=0.025, @@ -29,7 +29,7 @@ mgsem_cd_opt <- function(x, opt_fun_args, tol=1e-4, eps_approx=1e-20, x0 <- x - for (pp in 1:NP){ + for (pp in 1L:NP){ partable <- cd_fun_args$partable x1 <- x[pp] interval <- x1 + interval_length*c(-1,1) diff --git a/R/mgsem_create_index.R b/R/mgsem_create_index.R index eeec533..6fd2065 100644 --- a/R/mgsem_create_index.R +++ b/R/mgsem_create_index.R @@ -1,12 +1,12 @@ ## File Name: mgsem_create_index.R -## File Version: 0.03 +## File Version: 0.04 mgsem_create_index <- function(x, all=TRUE, start=0, symm=FALSE, onlydiag=FALSE) { ND <- prod(dim(x)) if (start>0){ - v <- matrix( start + ( 1:ND ) - 1, nrow=dim(x)[1], ncol=dim(x)[2]) + v <- matrix( start + ( 1L:ND ) - 1, nrow=dim(x)[1], ncol=dim(x)[2]) if (symm){ v <- v + t(v) v <- as.vector(v) @@ -23,7 +23,7 @@ mgsem_create_index <- function(x, all=TRUE, start=0, symm=FALSE, onlydiag=FALSE) if (onlydiag){ ND <- dim(x)[1] if (start>0){ - v <- diag( start - 1 + 1:ND ) + v <- diag( start - 1 + 1L:ND ) } else { v <- diag( rep(1, dim(x)[1])) } diff --git a/R/mgsem_duplication_matrix.R b/R/mgsem_duplication_matrix.R index 4180bbf..07a0734 100644 --- a/R/mgsem_duplication_matrix.R +++ b/R/mgsem_duplication_matrix.R @@ -1,16 +1,16 @@ ## File Name: mgsem_duplication_matrix.R -## File Version: 0.11 +## File Version: 0.121 mgsem_duplication_matrix <- function(x) { NV <- ncol(x) - dfr1 <- data.frame( index1=rep(1:NV, NV), index2=rep(1:NV, each=NV)) + dfr1 <- data.frame( index1=rep(1L:NV, NV), index2=rep(1L:NV, each=NV)) N1 <- nrow(dfr1) dfr2 <- dfr1[ dfr1[,1] >=dfr1[,2], ] ND <- nrow(dfr2) dupl <- matrix(0, nrow=ND, ncol=N1) dd <- 1 - for (dd in 1:ND){ + for (dd in 1L:ND){ h1 <- dfr2$index1[dd] h2 <- dfr2$index2[dd] i1 <- which( ( dfr1$index1==h1 ) & ( dfr1$index2==h2 ) ) diff --git a/R/mgsem_evaluate_penalties.R b/R/mgsem_evaluate_penalties.R index c49a230..6269ebc 100644 --- a/R/mgsem_evaluate_penalties.R +++ b/R/mgsem_evaluate_penalties.R @@ -1,5 +1,5 @@ ## File Name: mgsem_evaluate_penalties.R -## File Version: 0.345 +## File Version: 0.347 mgsem_evaluate_penalties <- function(x, partable, prior_list, technical, @@ -8,7 +8,7 @@ mgsem_evaluate_penalties <- function(x, partable, prior_list, technical, { ND <- nrow(partable) if (is.null(loop_parms)){ - loop_parms <- (1:ND)[ partable$unique==1] + loop_parms <- (1L:ND)[ partable$unique==1] } NP <- max(partable$index) if (!deriv){ @@ -137,7 +137,7 @@ mgsem_evaluate_penalties <- function(x, partable, prior_list, technical, der_z <- val NP <- length(x) val <- rep(0,NP) - for (hh in 1:NDP){ + for (hh in 1L:NDP){ i1 <- diffpar_pen_list_entries$index1[hh] val[i1] <- val[i1] + der_z[hh] i2 <- diffpar_pen_list_entries$index2[hh] diff --git a/R/mgsem_grad_fun.R b/R/mgsem_grad_fun.R index 6f97e86..4a1d3f3 100644 --- a/R/mgsem_grad_fun.R +++ b/R/mgsem_grad_fun.R @@ -1,5 +1,5 @@ ## File Name: mgsem_grad_fun.R -## File Version: 0.172 +## File Version: 0.173 mgsem_grad_fun <- function(x, opt_fun_args, output_all=FALSE) @@ -22,7 +22,7 @@ mgsem_grad_fun <- function(x, opt_fun_args, output_all=FALSE) est_total0 <- res$est_total0 dermoments0 <- list() - for (gg in 1:G){ + for (gg in 1L:G){ grad_suffstat_fun_args <- list(suffstat=opt_fun_args$suffstat[[gg]], Mu=implied0[[gg]]) if (estimator=='ME'){ @@ -36,7 +36,7 @@ mgsem_grad_fun <- function(x, opt_fun_args, output_all=FALSE) args=grad_suffstat_fun_args) } - for (dd in 1:ND){ + for (dd in 1L:ND){ group_dd <- partable$group[dd] i1_dd <- partable$i1[dd] diff --git a/R/mgsem_list_elements_est_total_implied.R b/R/mgsem_list_elements_est_total_implied.R index 6bb4e87..539f628 100644 --- a/R/mgsem_list_elements_est_total_implied.R +++ b/R/mgsem_list_elements_est_total_implied.R @@ -1,12 +1,12 @@ ## File Name: mgsem_list_elements_est_total_implied.R -## File Version: 0.02 +## File Version: 0.03 mgsem_list_elements_est_total_implied <- function(model, is_B) { G <- length(model)-1 est_total0 <- list() implied0 <- list() - for (gg in 1:G){ + for (gg in 1L:G){ est0 <- model[[1]]$est est_gg <- model[[gg+1]]$est est_total0[[gg]] <- mgsem_add_list_entries(list1=est0, add_list=est_gg, diff --git a/R/mgsem_loglike_suffstat_derivative.R b/R/mgsem_loglike_suffstat_derivative.R index 30c6946..63b2972 100644 --- a/R/mgsem_loglike_suffstat_derivative.R +++ b/R/mgsem_loglike_suffstat_derivative.R @@ -1,5 +1,5 @@ ## File Name: mgsem_loglike_suffstat_derivative.R -## File Version: 0.171 +## File Version: 0.172 mgsem_loglike_suffstat_derivative <- function(suffstat, Mu, Sigma ) @@ -20,22 +20,10 @@ mgsem_loglike_suffstat_derivative <- function(suffstat, Mu, Sigma ) dermean <- as.vector( N*( crossprod(m1, S1 ))) #*** covariance matrix - # dercov <- dercov0 <- matrix(0, nrow=p, ncol=p) y <- S1 %*% m1 S2 <- S %*% S1 S3 <- S1 %*% S2 - # for (ii in 1:p){ - # for (jj in ii:p){ - # # Sigma_der <- add_increment(dercov0, 1, ii,jj, symm=TRUE) - # # t2 <- -( t(y) %*% Sigma_der %*% y )[1,1] - # t2 <- -ifelse(ii==jj, y[ii]^2, y[ii]*y[jj]+y[jj]*y[ii]) - # # t3 <- - sum( diag( S3 %*% Sigma_der ) ) - # t3 <- -ifelse(ii==jj, S3[ii,jj], 2*S3[ii,jj] ) - # # t4 <- sum(diag(S1 %*% Sigma_der)) - # t4 <- ifelse(ii==jj, S1[ii,jj], 2*S1[ii,jj] ) - # dercov[ii,jj] <- dercov[jj,ii] <- -N/2*(t2+t3+t4) - # } - # } + #-- output res <- list(dermean=dermean, y=y, S1=S1, S2=S2, S3=S3) return(res) diff --git a/R/mgsem_opt_fun.R b/R/mgsem_opt_fun.R index b088b91..bb937da 100644 --- a/R/mgsem_opt_fun.R +++ b/R/mgsem_opt_fun.R @@ -1,5 +1,5 @@ ## File Name: mgsem_opt_fun.R -## File Version: 0.274 +## File Version: 0.275 mgsem_opt_fun <- function(x, opt_fun_args, output_all=FALSE) @@ -86,7 +86,7 @@ mgsem_opt_fun <- function(x, opt_fun_args, output_all=FALSE) if (output_all){ # chi square statistic and RMSEA p_mu <- 0 - for (gg in 1:G){ + for (gg in 1L:G){ mu1 <- suffstat[[gg]]$M p_mu <- p_mu + sum(abs(mu1)>1e-14) } diff --git a/R/mgsem_output_proc_casewise_likelihood.R b/R/mgsem_output_proc_casewise_likelihood.R index 2e2ea77..f6fa4cf 100644 --- a/R/mgsem_output_proc_casewise_likelihood.R +++ b/R/mgsem_output_proc_casewise_likelihood.R @@ -1,5 +1,5 @@ ## File Name: mgsem_output_proc_casewise_likelihood.R -## File Version: 0.07 +## File Version: 0.08 mgsem_output_proc_casewise_likelihood <- function(data_proc, implied, estimator="ML") @@ -12,7 +12,7 @@ mgsem_output_proc_casewise_likelihood <- function(data_proc, implied, estimator= data <- data_proc$data idgroup <- data_proc$idgroup case_ll <- rep(NA,N) - for (gg in 1:G){ + for (gg in 1L:G){ implied_gg <- implied[[gg]] ind_gg <- which(idgroup==gg) Mu <- as.vector(implied_gg$Mu[,1]) diff --git a/R/mgsem_output_proc_residuals.R b/R/mgsem_output_proc_residuals.R index c29bb70..e75d09e 100644 --- a/R/mgsem_output_proc_residuals.R +++ b/R/mgsem_output_proc_residuals.R @@ -1,5 +1,5 @@ ## File Name: mgsem_output_proc_residuals.R -## File Version: 0.05 +## File Version: 0.06 mgsem_output_proc_residuals <- function(implied, suffstat) @@ -7,7 +7,7 @@ mgsem_output_proc_residuals <- function(implied, suffstat) G <- length(implied) residuals_groupwise <- list() fit_total <- list(fit_Mu=0, fit_Sigma=0) - for (gg in 1:G){ + for (gg in 1L:G){ N <- suffstat[[gg]]$N res1 <- list() res1$Mu <- suffstat[[gg]]$M-implied[[gg]]$Mu diff --git a/R/mgsem_partable2model.R b/R/mgsem_partable2model.R index fc37656..922604c 100644 --- a/R/mgsem_partable2model.R +++ b/R/mgsem_partable2model.R @@ -1,5 +1,5 @@ ## File Name: mgsem_partable2model.R -## File Version: 0.151 +## File Version: 0.152 mgsem_partable2model <- function(partable, model, index=FALSE) @@ -10,7 +10,7 @@ mgsem_partable2model <- function(partable, model, index=FALSE) entries <- c('est','index') } for (entry in entries){ - for (dd in 1:ND){ + for (dd in 1L:ND){ hh <- partable[dd,'group']+1 type <- paste(partable[dd,'type']) mat <- model[[hh]][[entry]][[type]] diff --git a/R/mgsem_proc_data.R b/R/mgsem_proc_data.R index c80914c..81696bc 100644 --- a/R/mgsem_proc_data.R +++ b/R/mgsem_proc_data.R @@ -1,5 +1,5 @@ ## File Name: mgsem_proc_data.R -## File Version: 0.055 +## File Version: 0.057 mgsem_proc_data <- function(data, group, weights) { @@ -14,7 +14,7 @@ mgsem_proc_data <- function(data, group, weights) } idgroup <- match(group, groups) suffstat <- list() - for (gg in 1:G){ + for (gg in 1L:G){ ind_gg <- which(group==groups[gg]) dat_gg <- data[ ind_gg, ] w <- weights[ind_gg] diff --git a/R/mgsem_proc_model.R b/R/mgsem_proc_model.R index bfe520d..80761b7 100644 --- a/R/mgsem_proc_model.R +++ b/R/mgsem_proc_model.R @@ -1,5 +1,5 @@ ## File Name: mgsem_proc_model.R -## File Version: 0.307 +## File Version: 0.308 mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group, W, prior_list=NULL, pen_type="lasso", fixed_parms=FALSE, @@ -27,7 +27,7 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group, W, model <- mgsem_proc_model_single_group(model=model) #--- loop over groups - for (gg in 0:G){ + for (gg in 0L:G){ group <- gg hh <- gg+1 @@ -54,7 +54,7 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group, W, M2 <- index[[type]] n1 <- nrow(M1) n2 <- ncol(M2) - for (ii in 1:n1){ + for (ii in 1L:n1){ if (symm){ hh <- ii } else { @@ -135,8 +135,8 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group, W, diffpar_pen$coef_indices <- coef_indices dp1 <- NULL NW <- ncol(W) - for (ww in 1:NW){ - for (uu in 1:NW){ + for (ww in 1L:NW){ + for (uu in 1L:NW){ val <- W[ww,uu] if (abs(val) > 1e-14){ dp2 <- data.frame(index1=ww, index2=uu, W=val) @@ -172,7 +172,7 @@ mgsem_proc_model <- function(model, G=G, random_sd=1e-1, technical, N_group, W, model <- mgsem_partable2model(partable=dfr, model=model, index=TRUE) #*** unique parameters - loop_parms <- (1:ND)[ dfr$unique==1] + loop_parms <- (1L:ND)[ dfr$unique==1] #- rewrite penalty parameters into model matrices entries <- c('pen_l2', 'pen_lp', 'pen_difflp') diff --git a/R/mgsem_proc_model_is_B.R b/R/mgsem_proc_model_is_B.R index 7de5723..872206f 100644 --- a/R/mgsem_proc_model_is_B.R +++ b/R/mgsem_proc_model_is_B.R @@ -1,12 +1,12 @@ ## File Name: mgsem_proc_model_is_B.R -## File Version: 0.03 +## File Version: 0.04 mgsem_proc_model_is_B <- function(model) { H <- length(model) is_B <- 0 - for (hh in 1:H){ + for (hh in 1L:H){ is_B <- is_B + ( ! is.null( model[[hh]][['est']][['B']] ) ) } is_B <- ( is_B > 0 ) diff --git a/R/mgsem_proc_model_partable_define_index.R b/R/mgsem_proc_model_partable_define_index.R index d9f5350..9d5d631 100644 --- a/R/mgsem_proc_model_partable_define_index.R +++ b/R/mgsem_proc_model_partable_define_index.R @@ -1,5 +1,5 @@ ## File Name: mgsem_proc_model_partable_define_index.R -## File Version: 0.121 +## File Version: 0.123 mgsem_proc_model_partable_define_index <- function(partable) @@ -9,13 +9,13 @@ mgsem_proc_model_partable_define_index <- function(partable) N1 <- length(ind1) M1 <- max(dfr$index)+10 if (length(ind1)>0){ - dfr$index[ind1] <- M1+1:N1 + dfr$index[ind1] <- M1+1L:N1 } dfr$index <- match(dfr$index, unique(dfr$index)) dfr$unique <- 1 * (! duplicated(dfr$index) ) #** find duplicated parameters ND <- nrow(dfr) - for (dd in 1:ND){ + for (dd in 1L:ND){ if (dfr[dd,'unique']==0){ i1 <- which( dfr[,'index']==dfr[dd,'index']) i2 <- which( dfr[,'unique']==1) diff --git a/R/mgsem_proc_model_update_penalties_matrix.R b/R/mgsem_proc_model_update_penalties_matrix.R index cc31e14..8a84264 100644 --- a/R/mgsem_proc_model_update_penalties_matrix.R +++ b/R/mgsem_proc_model_update_penalties_matrix.R @@ -1,12 +1,12 @@ ## File Name: mgsem_proc_model_update_penalties_matrix.R -## File Version: 0.072 +## File Version: 0.073 mgsem_proc_model_update_penalties_matrix <- function(partable, entries, model) { ND <- nrow(partable) for (entry in entries){ - for (dd in 1:ND){ + for (dd in 1L:ND){ if (partable[dd,'unique']==1){ group <- partable$group[dd]+1 type_dd <- paste(partable$type[dd]) diff --git a/R/mgsem_proc_suffstat.R b/R/mgsem_proc_suffstat.R index bd1163f..8df48e7 100644 --- a/R/mgsem_proc_suffstat.R +++ b/R/mgsem_proc_suffstat.R @@ -1,12 +1,12 @@ ## File Name: mgsem_proc_suffstat.R -## File Version: 0.098 +## File Version: 0.099 mgsem_proc_suffstat <- function(suffstat) { G <- length(suffstat) N_group <- rep(0,G) - for (gg in 1:G){ + for (gg in 1L:G){ suffstat_gg <- suffstat[[gg]] N_group[gg] <- suffstat_gg$N # rename entries diff --git a/R/mgsem_suffstat_covariance_matrix.R b/R/mgsem_suffstat_covariance_matrix.R index b2bef27..6ff607d 100644 --- a/R/mgsem_suffstat_covariance_matrix.R +++ b/R/mgsem_suffstat_covariance_matrix.R @@ -1,13 +1,12 @@ ## File Name: mgsem_suffstat_covariance_matrix.R -## File Version: 0.099 +## File Version: 0.101 mgsem_suffstat_covariance_matrix <- function(suffstat) { G <- length(suffstat) dfr <- NULL - # gg <- 1 - for (gg in 1:G){ + for (gg in 1L:G){ suffstat_gg <- suffstat[[gg]] names_gg <- names(suffstat_gg) mu <- suffstat_gg[[ intersect( c('mu','M'), names_gg) ]] @@ -16,23 +15,23 @@ mgsem_suffstat_covariance_matrix <- function(suffstat) vars <- names(mu) I1 <- length(mu) if (is.null(vars)){ - vars <- paste0('V',1:I1) + vars <- paste0('V',1L:I1) } #** mu label <- paste0('mu_G', gg, '_', vars) se <- sqrt( diag(Sigma) / N ) - dfr1 <- data.frame(type='mu', group=gg, index1=1:I1, index2=NA, label=label, + dfr1 <- data.frame(type='mu', group=gg, index1=1L:I1, index2=NA, label=label, par=mu, se=se) V1 <- Sigma / N rownames(V1) <- colnames(V1) <- label #** sigma - c1 <- rbind( 1:I1, 1:I1) + c1 <- rbind( 1L:I1, 1L:I1) c2 <- utils::combn(x=I1, m=2) c3 <- rbind( t(c1), t(c2)) c3 <- c3[ order(c3[,1]), ] label <- paste0('Sigma_G', gg, '_', vars[c3[,1]], '_', vars[c3[,2]] ) - par <- Sigma[ c3[,1:2] ] + par <- Sigma[ c3[,1L:2] ] #* compute duplication matrix K <- mgsem_duplication_matrix(x=Sigma) diff --git a/R/mgsem_vcov_me.R b/R/mgsem_vcov_me.R index 4320bc6..97a71d5 100644 --- a/R/mgsem_vcov_me.R +++ b/R/mgsem_vcov_me.R @@ -1,5 +1,5 @@ ## File Name: mgsem_vcov_me.R -## File Version: 0.169 +## File Version: 0.171 mgsem_vcov_me <- function(coef, opt_fun_args, suffstat_vcov, comp_se, se_delta_formula=FALSE) @@ -18,7 +18,7 @@ mgsem_vcov_me <- function(coef, opt_fun_args, suffstat_vcov, comp_se, grad_der_coef <- matrix(0, nrow=NP, ncol=NP) rownames(grad_der_coef) <- colnames(grad_der_coef) <- parnames h <- opt_fun_args$technical$h - for (pp in 1:NP){ + for (pp in 1L:NP){ coef1 <- mgsem_add_increment(x=coef, h=h, i1=pp) res1 <- mgsem_grad_fun(x=coef1, opt_fun_args=opt_fun_args, output_all=FALSE) coef2 <- mgsem_add_increment(x=coef, h=-h, i1=pp) @@ -38,11 +38,11 @@ mgsem_vcov_me <- function(coef, opt_fun_args, suffstat_vcov, comp_se, colnames(grad_der_suffstat) <- suffstat_pars$label opt_fun_args1 <- opt_fun_args suffstat <- opt_fun_args$suffstat - for (pp in 1:SP){ + for (pp in 1L:SP){ suffstat_pars_pp <- suffstat_pars[pp,] group_pp <- suffstat_pars_pp$group val_pp <- list(NA,2) - for (oo in 1:2){ + for (oo in 1L:2){ suffstat1 <- suffstat u <- h if (oo==2){ diff --git a/R/noharm_sirt_optim_function_R.R b/R/noharm_sirt_optim_function_R.R index a12a0d8..291a5bd 100644 --- a/R/noharm_sirt_optim_function_R.R +++ b/R/noharm_sirt_optim_function_R.R @@ -1,16 +1,17 @@ ## File Name: noharm_sirt_optim_function_R.R -## File Version: 0.02 +## File Version: 0.05 noharm_sirt_optim_function_R <- function(gamma_val, delta, I, wgtm, pm, b0.jk, b1.jk, b2.jk, b3.jk) { val <- 0 - for (ii in 1:(I-1)){ + for (ii in 1L:(I-1)){ for (jj in (ii+1):I){ - if (wgtm[ii,jj] >0 ){ + if (wgtm[ii,jj]>0 ){ x_ij <- gamma_val[ii,jj] / sqrt( delta[ii] * delta[jj] ) - pm_exp <- b0.jk[ii,jj] + b1.jk[ii,jj]*x_ij + b2.jk[ii,jj]*x_ij^2 + b3.jk[ii,jj]*x_ij^3 + pm_exp <- b0.jk[ii,jj] + b1.jk[ii,jj]*x_ij + b2.jk[ii,jj]*x_ij^2 + + b3.jk[ii,jj]*x_ij^3 val <- val + wgtm[ii,jj]*(pm[ii,jj] - pm_exp)^2 } } diff --git a/R/rasch.copula3.R b/R/rasch.copula3.R index 582c59f..bd0109a 100644 --- a/R/rasch.copula3.R +++ b/R/rasch.copula3.R @@ -1,5 +1,5 @@ ## File Name: rasch.copula3.R -## File Version: 6.53 +## File Version: 6.542 @@ -464,7 +464,7 @@ rasch.copula3 <- function( dat, itemcluster, dims=NULL, ll1[cc] <- sum( rjkCC[[cc]] * log( rest1$pjk.theta.kCC[[cc]] ) ) ll2[cc] <- sum( rjkCC[[cc]] * log( rest2$pjk.theta.kCC[[cc]] ) ) } - a1 <- aggregate( cbind( ll0, ll1, ll2 ), list(est.delta), sum, na.rm=T) + a1 <- stats::aggregate( cbind( ll0, ll1, ll2 ), list(est.delta), sum, na.rm=T) ll0 <- a1[,2] ll1 <- a1[,3] ll2 <- a1[,4] @@ -573,7 +573,7 @@ rasch.copula3 <- function( dat, itemcluster, dims=NULL, ll0 <- rescop$ll # mu + h w1 <- wgt.theta - w2 <- dnorm( theta.k, mean=mu[gg] + h, sd=sigma[gg] ) + w2 <- stats::dnorm( theta.k, mean=mu[gg] + h, sd=sigma[gg] ) w1[,gg] <- w2 / sum(w2) rescop <- .ll.rasch.copula20( theta.k, b, alpha1, alpha2, a, dat2.li, itemcluster0, CC, dp.ld, dat2.ld, dat3.ld, dat2.ld.resp, dat2.li.resp, delta, @@ -621,7 +621,7 @@ rasch.copula3 <- function( dat, itemcluster, dims=NULL, prbar <- floor( prbar ) prbar <- c( prbar[1], diff(prbar) ) if (progress){ cat(" Estimation of sigma: |") } - for (gg in 2:G){ + for (gg in 2L:G){ # est.aa <- est.a * (est.a==aa ) rescop <- .ll.rasch.copula20( theta.k, b, alpha1, alpha2, a, dat2.li, itemcluster0, CC, dp.ld, dat2.ld, dat3.ld, dat2.ld.resp, dat2.li.resp, delta, wgt.theta, I, @@ -653,7 +653,7 @@ rasch.copula3 <- function( dat, itemcluster, dims=NULL, sigma.change <- ifelse( abs( sigma.change ) > .3, .3*sign(sigma.change), sigma.change ) sigma.change <- sigma.change * ( ( 1:G )==gg ) sigma <- sigma + sigma.change - w2 <- dnorm( theta.k, mean=mu[gg], sd=sigma[gg] ) + w2 <- stats::dnorm( theta.k, mean=mu[gg], sd=sigma[gg] ) wgt.theta[,gg] <- w2 / sum(w2) # cat( aa, " ") ; if (progress){ diff --git a/R/rmvn.R b/R/rmvn.R index 0499899..9845309 100644 --- a/R/rmvn.R +++ b/R/rmvn.R @@ -1,5 +1,5 @@ ## File Name: rmvn.R -## File Version: 0.04 +## File Version: 0.05 rmvn <- function(N, mu, Sigma, exact=TRUE) { @@ -10,7 +10,7 @@ rmvn <- function(N, mu, Sigma, exact=TRUE) #-- compute data with exact zero means and identity covariance matrix if (exact){ - c1 <- stats::cov.wt(dat0, method="ML") + c1 <- stats::cov.wt(dat0, method='ML') dat0 <- dat0 - matrix( c1$center, nrow=N, ncol=D, byrow=TRUE) COV0 <- c1$cov c0 <- svd(COV0) diff --git a/R/sirt_EAP.R b/R/sirt_EAP.R index 8fea90e..c074626 100644 --- a/R/sirt_EAP.R +++ b/R/sirt_EAP.R @@ -1,5 +1,5 @@ ## File Name: sirt_EAP.R -## File Version: 0.01 +## File Version: 0.02 sirt_EAP <- function(post, theta) { @@ -7,7 +7,7 @@ sirt_EAP <- function(post, theta) N <- nrow(post) D <- ncol(theta) EAP <- matrix(NA, nrow=N, ncol=D) - for (tt in 1:D){ + for (tt in 1L:D){ thetaM <- sirt_matrix2(theta[,tt], nrow=N) EAP[,tt] <- rowSums(thetaM * post) } diff --git a/R/sirt_MAP.R b/R/sirt_MAP.R index 842ba1d..283b243 100644 --- a/R/sirt_MAP.R +++ b/R/sirt_MAP.R @@ -1,12 +1,12 @@ ## File Name: sirt_MAP.R -## File Version: 0.01 +## File Version: 0.02 sirt_MAP <- function(post, theta) { TP <- ncol(post) maxval <- post[,1] indval <- 1 - for (tt in 2:TP){ + for (tt in 2L:TP){ m0 <- maxval maxval <- ifelse( post[,tt] > m0, post[,tt], m0 ) indval <- ifelse( post[,tt] > m0, tt, indval ) diff --git a/R/sirt_add_list_elements.R b/R/sirt_add_list_elements.R index 607e3b3..df1b62d 100644 --- a/R/sirt_add_list_elements.R +++ b/R/sirt_add_list_elements.R @@ -1,10 +1,10 @@ ## File Name: sirt_add_list_elements.R -## File Version: 0.02 +## File Version: 0.04 sirt_add_list_elements <- function(res, res2) { NR <- length(res2) - for (rr in 1:NR){ + for (rr in 1L:NR){ res[[ names(res2)[[rr]] ]] <- res2[[rr]] } return(res) diff --git a/R/sirt_logit_to_probs.R b/R/sirt_logit_to_probs.R index 84f7c12..d04e1f6 100644 --- a/R/sirt_logit_to_probs.R +++ b/R/sirt_logit_to_probs.R @@ -1,12 +1,12 @@ ## File Name: sirt_logit_to_probs.R -## File Version: 0.03 +## File Version: 0.04 sirt_logit_to_probs <- function(y) { K1 <- length(y) x <- rep(0,K1+1) - x[1:K1] <- y + x[1L:K1] <- y x <- exp(x) x <- x / sum(x) return(x) diff --git a/R/sirt_matrix_lower_to_upper.R b/R/sirt_matrix_lower_to_upper.R index 1104cb9..9630ada 100644 --- a/R/sirt_matrix_lower_to_upper.R +++ b/R/sirt_matrix_lower_to_upper.R @@ -1,10 +1,10 @@ ## File Name: sirt_matrix_lower_to_upper.R -## File Version: 0.02 +## File Version: 0.03 sirt_matrix_lower_to_upper <- function(x) { N <- ncol(x) - for (ii in 1:N){ + for (ii in 1L:N){ for (jj in ii:N){ x[ii,jj] <- x[jj,ii] } diff --git a/R/sirt_optimizer.R b/R/sirt_optimizer.R index c7eb06d..0238db0 100644 --- a/R/sirt_optimizer.R +++ b/R/sirt_optimizer.R @@ -1,5 +1,5 @@ ## File Name: sirt_optimizer.R -## File Version: 0.371 +## File Version: 0.372 sirt_optimizer <- function(optimizer, par, fn, grad=NULL, method="L-BFGS-B", hessian=TRUE, control=list(), ...) @@ -35,7 +35,7 @@ sirt_optimizer <- function(optimizer, par, fn, grad=NULL, method="L-BFGS-B", res0 <- do.call(what=optimx_fun, args=args) # res <- optimx::Rvmmin(par=par, fn=fn, gr=grad, control=control, ...) res <- list() - m1 <- as.vector(as.numeric(res0[1,1:NP])) + m1 <- as.vector(as.numeric(res0[1,1L:NP])) res$par <- m1 res$convergence <- res0$convcode[1] res$iter <- res0$feval[1] diff --git a/R/sirt_pem_collect_parameters.R b/R/sirt_pem_collect_parameters.R index 4523e66..8e8b1c1 100644 --- a/R/sirt_pem_collect_parameters.R +++ b/R/sirt_pem_collect_parameters.R @@ -1,5 +1,5 @@ ## File Name: sirt_pem_collect_parameters.R -## File Version: 0.06 +## File Version: 0.07 sirt_pem_collect_parameters <- function( parmlist, pem_parameter_index ) { @@ -7,7 +7,7 @@ sirt_pem_collect_parameters <- function( parmlist, pem_parameter_index ) parm <- rep(0, NV) NP <- pem_parameter_index[['__n_parameters']] NPL <- length(parmlist) - for (pp in 1:NPL){ + for (pp in 1L:NPL){ var_pp <- names(parmlist)[pp] parm_index_pp <- pem_parameter_index[[ var_pp ]] x <- as.vector(parmlist[[pp]]) diff --git a/R/sirt_pem_create_parameter_index.R b/R/sirt_pem_create_parameter_index.R index a7fa47c..5eb5d11 100644 --- a/R/sirt_pem_create_parameter_index.R +++ b/R/sirt_pem_create_parameter_index.R @@ -1,5 +1,5 @@ ## File Name: sirt_pem_create_parameter_index.R -## File Version: 0.11 +## File Version: 0.121 sirt_pem_create_parameter_index <- function( parmlist ) @@ -9,7 +9,7 @@ sirt_pem_create_parameter_index <- function( parmlist ) parmlist_names <- names(parmlist) parm_index_names <- rep('', NP) last_index <- 0 - for (pp in 1:NP){ + for (pp in 1L:NP){ x <- parmlist[[pp]] x_pp <- list() parm_index_names[pp] <- parmlist_names[pp] diff --git a/R/sirt_rbind_fill.R b/R/sirt_rbind_fill.R index cd3ff6d..27eb022 100644 --- a/R/sirt_rbind_fill.R +++ b/R/sirt_rbind_fill.R @@ -1,5 +1,5 @@ ## File Name: sirt_rbind_fill.R -## File Version: 0.05 +## File Version: 0.06 ## reimplementation of plyr::rbind.fill @@ -11,7 +11,7 @@ sirt_rbind_fill <- function( x, y ) z <- matrix( NA, nrow=nx+ny, ncol=length(vars) ) colnames(z) <- vars z <- as.data.frame(z) - z[ 1:nx, colnames(x) ] <- x - z[ nx + 1:ny, colnames(y) ] <- y + z[ 1L:nx, colnames(x) ] <- x + z[ nx + 1L:ny, colnames(y) ] <- y return(z) } diff --git a/R/sirt_rename_list_entry.R b/R/sirt_rename_list_entry.R index 4e3e001..04ca3ef 100644 --- a/R/sirt_rename_list_entry.R +++ b/R/sirt_rename_list_entry.R @@ -1,14 +1,14 @@ ## File Name: sirt_rename_list_entry.R -## File Version: 0.09 +## File Version: 0.101 sirt_rename_list_entry <- function(model, name1, name2) { model_temp <- model G <- length(model) - for (gg in 1:G){ + for (gg in 1L:G){ model_gg <- model_temp[[gg]] TP <- length(model_gg) - for (tt in 1:TP){ + for (tt in 1L:TP){ model_gg_tt <- model_gg[[tt]] ind <- which( names(model_gg_tt)==name1) if (length(ind)>0){ diff --git a/R/smirt_postproc.R b/R/smirt_postproc.R index 09b28e8..2380d03 100644 --- a/R/smirt_postproc.R +++ b/R/smirt_postproc.R @@ -1,29 +1,32 @@ ## File Name: smirt_postproc.R -## File Version: 0.11 +## File Version: 0.123 -################################################################# -# person parameter estimates +#*** person parameter estimates .smirt.person.parameters <- function( data, D, theta.k, - p.xi.aj, p.aj.xi, weights ){ - #************************** - person <- data.frame("case"=1:(nrow(data)), "M"=rowMeans( data, na.rm=TRUE) ) + p.xi.aj, p.aj.xi, weights ) +{ + + person <- data.frame(case=1L:(nrow(data)), M=rowMeans( data, na.rm=TRUE)) EAP.rel <- rep(0,D) names(EAP.rel) <- colnames(theta.k) nstudl <- rep(1,nrow(data)) - for (dd in 1:D){ #dd <- 1 + for (dd in 1L:D){ #dd <- 1 dd1 <- colnames(theta.k)[dd] person$EAP <- rowSums( p.aj.xi * outer( nstudl, theta.k[,dd] ) ) - person$SE.EAP <- sqrt(rowSums( p.aj.xi * outer( nstudl, theta.k[,dd]^2 ) ) - person$EAP^2) + person$SE.EAP <- sqrt(rowSums( p.aj.xi * outer( nstudl, theta.k[,dd]^2 ) ) - + person$EAP^2) EAP.variance <- stats::weighted.mean( person$EAP^2, weights ) - ( stats::weighted.mean( person$EAP, weights ) )^2 EAP.error <- stats::weighted.mean( person$SE.EAP^2, weights ) EAP.rel[dd] <- EAP.variance / ( EAP.variance + EAP.error ) - colnames(person)[ which( colnames(person)=="EAP" ) ] <- paste("EAP.", dd1, sep="") - colnames(person)[ which( colnames(person)=="SE.EAP" ) ] <- paste("SE.EAP.", dd1, sep="") + colnames(person)[ which( colnames(person)=="EAP" ) ] <- + paste("EAP.", dd1, sep="") + colnames(person)[ which( colnames(person)=="SE.EAP" ) ] <- + paste("SE.EAP.", dd1, sep="") } if ( is.null( names(EAP.rel) ) ){ - names(EAP.rel) <- paste0( "Dim", 1:(length(EAP.rel)) ) - } + names(EAP.rel) <- paste0( "Dim", 1L:(length(EAP.rel)) ) + } # MLE mle.est <- theta.k[ max.col( p.xi.aj ),, drop=FALSE] colnames(mle.est) <- paste0( "MLE.", names(EAP.rel)) @@ -33,7 +36,7 @@ colnames(mle.est) <- paste0( "MAP.", names(EAP.rel)) person <- cbind( person, mle.est ) # results - res <- list( "person"=person, "EAP.rel"=EAP.rel ) + res <- list( person=person, EAP.rel=EAP.rel ) return(res) - } -################################################################# +} + diff --git a/R/summary.gom.em.R b/R/summary.gom.em.R index b4bbfcd..57dcdda 100644 --- a/R/summary.gom.em.R +++ b/R/summary.gom.em.R @@ -1,5 +1,5 @@ ## File Name: summary.gom.em.R -## File Version: 0.177 +## File Version: 0.179 #--- summary for gom object @@ -48,7 +48,8 @@ summary.gom <- function( object, file=NULL, ...) cat( "Number of estimated parameters", "=", object$ic$np, "\n" ) cat( " Number of estimated item parameters", "=", object$ic$np.item, "\n" ) - cat( " Number of estimated distribution parameters", "=", object$ic$np.trait, "\n\n" ) + cat( " Number of estimated distribution parameters", "=", + object$ic$np.trait, "\n\n" ) #--- information criteria rm_summary_information_criteria(object=object) diff --git a/R/summary.lsem.permutationTest.R b/R/summary.lsem.permutationTest.R index 60c65f9..4b1f604 100644 --- a/R/summary.lsem.permutationTest.R +++ b/R/summary.lsem.permutationTest.R @@ -1,5 +1,5 @@ ## File Name: summary.lsem.permutationTest.R -## File Version: 0.297 +## File Version: 0.298 summary.lsem.permutationTest <- function( object, file=NULL, digits=3, ... ) @@ -42,7 +42,7 @@ summary.lsem.permutationTest <- function( object, file=NULL, digits=3, ... ) cat('Global Test Statistics\n\n') obji <- object$teststat VV <- ncol(obji) - for (vv in 2:VV){ + for (vv in 2L:VV){ obji[,vv] <- round( obji[,vv], digits ) } print(obji) diff --git a/R/summary_round_helper.R b/R/summary_round_helper.R index dfd17d1..6ea7ca0 100644 --- a/R/summary_round_helper.R +++ b/R/summary_round_helper.R @@ -1,10 +1,10 @@ ## File Name: summary_round_helper.R -## File Version: 0.09 +## File Version: 0.101 summary_round_helper <- function( obji, digits, exclude=NULL, print=TRUE) { NC <- ncol(obji) - ind <- 1:NC + ind <- 1L:NC if ( ! is.null(exclude) ){ ind2 <- which( colnames(obji) %in% exclude ) ind <- setdiff( ind, ind2 ) diff --git a/R/write.fwf2.R b/R/write.fwf2.R index 1a93a31..148a8f3 100644 --- a/R/write.fwf2.R +++ b/R/write.fwf2.R @@ -1,17 +1,17 @@ ## File Name: write.fwf2.R -## File Version: 1.121 +## File Version: 1.124 write.fwf2 <- function( dat, format.full, format.round, savename ) { if (is.null( colnames(dat) ) ){ - colnames(dat) <- paste( 'V', 1:( ncol(dat) ), sep='') + colnames(dat) <- paste( 'V', 1L:(ncol(dat)), sep='') } matr <- matrix( ' ', nrow=nrow(dat), ncol=ncol(dat)) ind1 <- which( format.round <=0 ) format.full[ ind1 ] <- format.full[ind1] format.round[ ind1 ] <- format.round[ind1] - for (vv in 1:( ncol(matr) ) ){ + for (vv in 1L:(ncol(matr)) ){ fvv <- format.round[vv] fff <- format.full[vv] matr[,vv] <- write.format2( vec1=dat[,vv], ff=fff, fr=fvv ) @@ -20,8 +20,8 @@ write.fwf2 <- function( dat, format.full, format.round, savename ) if ( is.vector(matr) ){ writeLines( matr, paste( savename, '.dat', sep='') ) } else { - utils::write.table( matr, paste( savename, '.dat', sep=''), - row.names=FALSE, col.names=FALSE) + utils::write.table( matr, paste(savename, '.dat', sep=''), + row.names=FALSE, col.names=FALSE) } dfr <- data.frame( variable=colnames(dat), begin=c( 1, cumsum( format.full )[ - ncol(dat) ] + 1 ), diff --git a/R/xxirt.R b/R/xxirt.R index ac8f9b2..f892f8e 100644 --- a/R/xxirt.R +++ b/R/xxirt.R @@ -1,5 +1,5 @@ ## File Name: xxirt.R -## File Version: 1.112 +## File Version: 1.159 #--- user specified item response model @@ -7,10 +7,9 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, partable=NULL, customTheta=NULL, group=NULL, weights=NULL, globconv=1E-6, conv=1E-4, maxit=1000, mstep_iter=4, mstep_reltol=1E-6, maxit_nr=0, optimizer_nr="nlminb", - control_nr=list(trace=1), - h=1E-4, use_grad=TRUE, - verbose=TRUE, penalty_fun_item=NULL, - np_fun_item=NULL, verbose_index=NULL, + estimator="ML", control_nr=list(trace=1), + h=1E-4, use_grad=TRUE, verbose=TRUE, penalty_fun_item=NULL, + np_fun_item=NULL, pml_args=NULL, verbose_index=NULL, cv_kfold=0, cv_maxit=10) { #*** preliminaries @@ -63,6 +62,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, item_index <- res$item_index dat <- as.matrix(dat) + # create item list item_list <- xxirt_createItemList( customItems=customItems, itemtype=itemtype, items=items, partable=partable ) @@ -97,14 +97,16 @@ 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 em_args <- list( maxit=maxit, verbose1=verbose1, disp=disp, item_list=item_list, items=items, Theta=Theta, ncat=ncat, partable=partable, partable_index=partable_index, dat=dat, - resp_index=resp_index, dat_resp=dat_resp, dat_resp_bool=dat_resp_bool, - dat1=dat1, dat1_resp=dat1_resp, customTheta=customTheta, G=G, + I=ncol(dat), resp_index=resp_index, dat_resp=dat_resp, + dat_resp_bool=dat_resp_bool, dat1=dat1, dat1_resp=dat1_resp, + customTheta=customTheta, G=G, par0=par0, maxK=maxK, group_index=group_index, weights=weights, mstep_iter=mstep_iter, eps=eps, mstep_reltol=mstep_reltol, mstep_method=mstep_method, item_index=item_index, h=h, @@ -119,8 +121,20 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, em_args$verbose3 <- FALSE } + #-- settings for PML estimation + if (estimator=='PML'){ + do_nr <- TRUE + maxit_nr <- maxit + pml_args <- xxirt_nr_pml_preproc_data(em_args=em_args, pml_args=pml_args) + em_args$pml_args <- pml_args + } + #--- run EM algorithm - em_out <- res <- do.call(what=xxirt_em_algorithm, args=em_args) + if (estimator=='ML'){ + em_out <- res <- do.call(what=xxirt_em_algorithm, args=em_args) + } else { + em_out <- res <- NULL + } #--- collect EM output if (em_count==1){ @@ -150,11 +164,11 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, cv_loglike <- NA if (do_cv){ N <- nrow(dat) - cv_zone <- ( ( 0:(N-1) ) %% cv_kfold ) + 1 + cv_zone <- ( ( 0L:(N-1) ) %% cv_kfold ) + 1 cv_loglike <- 0 em_args$partable <- partable em_args$customTheta <- customTheta - for (jj in 1:cv_kfold){ + for (jj in 1L:cv_kfold){ weights1 <- weights*(cv_zone!=jj ) weights2 <- weights*(cv_zone==jj ) em_args$weights <- weights1 @@ -181,16 +195,22 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, if (do_nr){ res <- xxirt_newton_raphson(em_out=em_out, em_args=em_args, maxit_nr=maxit_nr, optimizer_nr=optimizer_nr, control_nr=control_nr, - verbose=verbose) + verbose=verbose, estimator=estimator) res_opt_nr <- res$res_opt_nr iter <- res_opt_nr$iter + iter iter_nr <- res_opt_nr$iter opt_values_nr <- res$opt_values + em_args <- res$em_args + pml_args <- res$em_args$pml_args em_args$partable <- partable <- res$partable + probs_items <- res$probs_items em_args$customTheta <- customTheta <- res$customTheta } # end NR optimization em_count <- 2 do_nr <- FALSE + if (estimator=='PML' | maxit==0){ + em_iterate <- FALSE + } } # loop for final EM iteration @@ -210,18 +230,26 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, #-- information criteria ic <- xxirt_ic( dev=dev, N=W, par_items=par_items, - par_Theta=par_Theta, I=I, par_items_bounds=par_items_bounds, - np_item=np_item) + par_Theta=par_Theta, I=I, par_items_bounds=par_items_bounds, + np_item=np_item, estimator=estimator) #-- compute EAP - EAP <- xxirt_EAP(p.aj.xi=p.aj.xi, Theta=Theta ) + if (estimator=='ML'){ + EAP <- xxirt_EAP(p.aj.xi=p.aj.xi, Theta=Theta ) + loglike <- -dev/2 + } + if (estimator=='PML'){ + EAP <- NULL + dev <- NULL + loglike <- NULL + } #--- output s2 <- Sys.time() res <- list( partable=partable, par_items=par_items, par_items_summary=par_items_summary, par_items_bounds=par_items_bounds, par_Theta=par_Theta, Theta=Theta, probs_items=probs_items, - probs_Theta=prior_Theta, deviance=dev, loglike=-dev/2, + probs_Theta=prior_Theta, deviance=dev, loglike=loglike, opt_val=opt_val, pen_val=pen_val, ic=ic, item_list=item_list, customItems=customItems, customTheta=customTheta, cv_loglike=cv_loglike, @@ -235,6 +263,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, resp_index=resp_index, converged=converged, res_opt_nr=res_opt_nr, opt_values_nr=opt_values_nr, maxit_nr=maxit_nr, iter=iter-1, iter_em=iter_em-1, iter_nr=iter_nr, + estimator=estimator, pml_args=pml_args, em_args=em_args, CALL=CALL, s1=s1, s2=s2 ) class(res) <- 'xxirt' return(res) diff --git a/R/xxirt_EAP.R b/R/xxirt_EAP.R index eaa632b..f271267 100644 --- a/R/xxirt_EAP.R +++ b/R/xxirt_EAP.R @@ -1,5 +1,5 @@ ## File Name: xxirt_EAP.R -## File Version: 0.116 +## File Version: 0.117 #--- compute EAP and its standard deviation @@ -7,12 +7,12 @@ xxirt_EAP <- function(p.aj.xi, Theta ) { D <- ncol(Theta) e1 <- p.aj.xi %*% Theta - colnames(e1) <- paste0('EAP.Dim',1:D) + colnames(e1) <- paste0('EAP.Dim', 1L:D) e2 <- p.aj.xi %*% Theta^2 - colnames(e2) <- paste0('SD.EAP.Dim',1:D) + colnames(e2) <- paste0('SD.EAP.Dim', 1L:D) e2 <- e2 - e1^2 res <- cbind( e1, e2) - res <- res[, rep( c(0,D), D ) + 1:D ] + res <- res[, rep( c(0,D), D ) + 1L:D ] return(res) } diff --git a/R/xxirt_classprobs_lca.R b/R/xxirt_classprobs_lca.R index f25e771..75f5ace 100644 --- a/R/xxirt_classprobs_lca.R +++ b/R/xxirt_classprobs_lca.R @@ -1,12 +1,12 @@ ## File Name: xxirt_classprobs_lca.R -## File Version: 0.08 +## File Version: 0.09 xxirt_classprobs_lca <- function(par, Theta, G) { K <- nrow(Theta) probs <- matrix(NA, nrow=K, ncol=G) for (gg in 1L:G){ - logitprobs_gg <- par[ (K-1)*(gg-1) + ( 1:(K-1) ) ] + logitprobs_gg <- par[ (K-1)*(gg-1) + ( 1L:(K-1) ) ] probs[,gg] <- xxirt_classprobs_lca_compute_probs(logitprobs=logitprobs_gg) } return(probs) diff --git a/R/xxirt_classprobs_lca_init_par.R b/R/xxirt_classprobs_lca_init_par.R index a188a0c..a91195f 100644 --- a/R/xxirt_classprobs_lca_init_par.R +++ b/R/xxirt_classprobs_lca_init_par.R @@ -1,12 +1,12 @@ ## File Name: xxirt_classprobs_lca_init_par.R -## File Version: 0.091 +## File Version: 0.092 xxirt_classprobs_lca_init_par <- function(K, G, random_sd=0, par_Theta_init=NULL) { - if (G==1){ + if (G==1L){ par_Theta <- xxirt_classprobs_lca_init_par_create(K=K, random_sd=random_sd) - names(par_Theta) <- paste0('pi',1:(K-1) ) + names(par_Theta) <- paste0('pi',1L:(K-1) ) if ( !is.null(par_Theta_init) ){ par_Theta <- par_Theta_init } @@ -14,7 +14,7 @@ xxirt_classprobs_lca_init_par <- function(K, G, random_sd=0, par_Theta_init=NULL v1 <- NULL for (gg in 1L:G){ par_Theta <- xxirt_classprobs_lca_init_par_create(K=K, random_sd=random_sd) - names(par_Theta) <- paste0('pi',1:(K-1), '_Gr', gg) + names(par_Theta) <- paste0('pi',1L:(K-1), '_Gr', gg) v1 <- c(v1, par_Theta) } par_Theta <- v1 diff --git a/R/xxirt_compute_itemprobs.R b/R/xxirt_compute_itemprobs.R index 622affd..3f8605f 100644 --- a/R/xxirt_compute_itemprobs.R +++ b/R/xxirt_compute_itemprobs.R @@ -1,5 +1,5 @@ ## File Name: xxirt_compute_itemprobs.R -## File Version: 0.214 +## File Version: 0.215 # compute item probabilities @@ -10,19 +10,19 @@ xxirt_compute_itemprobs <- function( item_list, items, Theta, ncat, maxK <- max(ncat) if ( is.null(item_index) ){ I <- length(items) - item_index <- 1:I + item_index <- 1L:I } I <- length(item_index) # compute item probabilities as a function of theta probs <- array( 0, dim=c(I,maxK,TP) ) - for (jj in 1:I){ + for (jj in 1L:I){ ii <- item_index[jj] item_ii <- item_list[[ii]] par_ii <- partable[ partable_index[[ii]], 'value' ] ncat_ii <- ncat[ii] arg_ii <- list( par=par_ii, Theta=Theta, ncat=ncat_ii ) probs_ii <- do.call( item_ii$P, arg_ii ) - probs[ jj, 1:ncat_ii,] <- t(probs_ii) + probs[ jj, 1L:ncat_ii,] <- t(probs_ii) } return(probs) } diff --git a/R/xxirt_compute_posterior.R b/R/xxirt_compute_posterior.R index 68b54ca..fee7831 100644 --- a/R/xxirt_compute_posterior.R +++ b/R/xxirt_compute_posterior.R @@ -1,5 +1,5 @@ ## File Name: xxirt_compute_posterior.R -## File Version: 0.293 +## File Version: 0.294 ##-- xxirt: compute posterior @@ -21,12 +21,12 @@ xxirt_compute_posterior <- function( prior_Theta, p.xi.aj, group, n.ik <- array( 0, dim=c(I,maxK, TP,G) ) N.ik <- array( 0, dim=c(I,maxK, TP) ) pi.k <- matrix( 0, nrow=TP, ncol=G ) - for (gg in 1:G){ + for (gg in 1L:G){ ind_gg <- group_index[[gg]] p.aj.xi.gg <- as.matrix( p.aj.xi[ind_gg, ] ) dat1_resp_gg <- dat1_resp[ ind_gg,, ] weights_gg <- weights[ind_gg] - for (kk in 1:maxK){ + for (kk in 1L:maxK){ dat1_resp_gg_kk <- as.matrix(dat1_resp_gg[,,kk]) n.ik[,kk,,gg] <- sirt_rcpp_xxirt_compute_posterior_expected_counts( dat1_resp_gg=dat1_resp_gg_kk, p_aj_xi_gg=p.aj.xi.gg, diff --git a/R/xxirt_compute_prior_Theta_from_x.R b/R/xxirt_compute_prior_Theta_from_x.R index 19a0513..f5cd105 100644 --- a/R/xxirt_compute_prior_Theta_from_x.R +++ b/R/xxirt_compute_prior_Theta_from_x.R @@ -1,9 +1,8 @@ ## File Name: xxirt_compute_prior_Theta_from_x.R -## File Version: 0.03 +## File Version: 0.07 xxirt_compute_prior_Theta_from_x <- function(x, em_args) { - # include parameters customTheta <- xxirt_parTheta_include_freeParameters( customTheta=em_args$customTheta, diff --git a/R/xxirt_compute_prob_item_from_x.R b/R/xxirt_compute_prob_item_from_x.R index c14b745..8de8709 100644 --- a/R/xxirt_compute_prob_item_from_x.R +++ b/R/xxirt_compute_prob_item_from_x.R @@ -1,13 +1,14 @@ ## File Name: xxirt_compute_prob_item_from_x.R -## File Version: 0.03 +## File Version: 0.06 -xxirt_compute_prob_item_from_x <- function(x, em_args) +xxirt_compute_prob_item_from_x <- function(x, em_args, item_index=NULL) { partable <- xxirt_partable_include_freeParameters( partable=em_args$partable, x=x[ em_args$parindex_items ] ) probs_items <- xxirt_compute_itemprobs( item_list=em_args$item_list, items=em_args$items, Theta=em_args$Theta, ncat=em_args$ncat, partable=partable, - partable_index=em_args$partable_index ) + partable_index=em_args$partable_index, + item_index=item_index) return(probs_items) } diff --git a/R/xxirt_createItemList.R b/R/xxirt_createItemList.R index 805bae3..1d59e47 100644 --- a/R/xxirt_createItemList.R +++ b/R/xxirt_createItemList.R @@ -1,20 +1,20 @@ ## File Name: xxirt_createItemList.R -## File Version: 0.13 +## File Version: 0.141 -######################################################### -# create item list + +#--- create item list xxirt_createItemList <- function( customItems, itemtype,items, partable ) { #*** list with customized items I <- length(items) - item_list <- as.list( 1:I ) + item_list <- as.list( 1L:I ) names(item_list) <- items CI <- length(customItems) - for (ii in 1:I){ + for (ii in 1L:I){ itemtype_ii <- itemtype[ii] item_list[[ii]] <- NA partable_ii <- partable[ partable$itemnr==ii, ] - for (vv in 1:CI){ + for (vv in 1L:CI){ if ( paste(customItems[[vv]]$name)==itemtype_ii ){ item_list[[ii]] <- customItems[[vv]] } # end if @@ -27,4 +27,3 @@ xxirt_createItemList <- function( customItems, itemtype,items, partable ) } # end ii return(item_list) } -############################################################### diff --git a/R/xxirt_createParTable.R b/R/xxirt_createParTable.R index febdd2d..b8742cb 100644 --- a/R/xxirt_createParTable.R +++ b/R/xxirt_createParTable.R @@ -1,5 +1,5 @@ ## File Name: xxirt_createParTable.R -## File Version: 0.258 +## File Version: 0.259 #*** create parameter table xxirt_createParTable <- function( dat, itemtype, customItems=NULL ) @@ -12,10 +12,10 @@ xxirt_createParTable <- function( dat, itemtype, customItems=NULL ) } dfr <- NULL CI <- length(customItems) - for (ii in 1:I){ + for (ii in 1L:I){ type_ii <- itemtype[ii] item_ii <- NULL - for ( vv in 1:CI){ + for ( vv in 1L:CI){ ci_ii <- customItems[[vv]] if ( ci_ii$name==type_ii ){ item_ii <- ci_ii @@ -49,7 +49,7 @@ xxirt_createParTable <- function( dat, itemtype, customItems=NULL ) } #**** create parameter indices NP <- nrow(dfr) - dfr$rowindex <- 1:NP + dfr$rowindex <- 1L:NP # parameter index dfr$parindex <- cumsum( dfr$est ) #*** parameter label diff --git a/R/xxirt_data_proc.R b/R/xxirt_data_proc.R index 8f46796..e6ed01b 100644 --- a/R/xxirt_data_proc.R +++ b/R/xxirt_data_proc.R @@ -1,5 +1,5 @@ ## File Name: xxirt_data_proc.R -## File Version: 0.213 +## File Version: 0.215 #-- data processing xxirt xxirt_data_proc <- function(dat, group=NULL, weights=NULL ) @@ -21,16 +21,16 @@ xxirt_data_proc <- function(dat, group=NULL, weights=NULL ) group <- match( group0, groups_unique ) maxK <- max(ncat) #*** group_index - group_index <- as.list( 1:G ) - for (gg in 1:G){ + group_index <- as.list( 1L:G ) + for (gg in 1L:G){ group_index[[gg]] <- which( group==gg ) } #*** data with response indices dat_na <- is.na(dat) dat_resp_bool <- ! dat_na dat_resp <- 1 - dat_na - resp_index <- as.list( 1:I ) - for ( ii in 1:I){ + resp_index <- as.list( 1L:I ) + for ( ii in 1L:I){ resp_index[[ii]] <- which( dat_resp[,ii]==1 ) } dat1 <- as.matrix(dat) diff --git a/R/xxirt_em_algorithm.R b/R/xxirt_em_algorithm.R index 7fdea80..01c9d3c 100644 --- a/R/xxirt_em_algorithm.R +++ b/R/xxirt_em_algorithm.R @@ -1,12 +1,12 @@ ## File Name: xxirt_em_algorithm.R -## File Version: 0.088 +## File Version: 0.092 xxirt_em_algorithm <- function(maxit, verbose1, verbose2, verbose3, disp, item_list, items, Theta, ncat, partable, partable_index, dat, resp_index, dat_resp, dat_resp_bool, dat1, dat1_resp, group, customTheta, G, par0, maxK, group_index, weights, mstep_iter, eps, mstep_reltol, mstep_method, item_index, h, use_grad, penalty_fun_item, par1, globconv, conv, - verbose_index) + verbose_index, I) { iter <- 1 dev <- 1E100 diff --git a/R/xxirt_em_args_extract.R b/R/xxirt_em_args_extract.R new file mode 100644 index 0000000..c4f125f --- /dev/null +++ b/R/xxirt_em_args_extract.R @@ -0,0 +1,11 @@ +## File Name: xxirt_em_args_extract.R +## File Version: 0.02 + +xxirt_em_args_extract <- function(em_args, em_out, object) +{ + v1 <- em_out[[object]] + if (is.null(v1)){ + v1 <- em_args[[object]] + } + return(v1) +} diff --git a/R/xxirt_ic.R b/R/xxirt_ic.R index af2e654..b509a2d 100644 --- a/R/xxirt_ic.R +++ b/R/xxirt_ic.R @@ -1,9 +1,10 @@ ## File Name: xxirt_ic.R -## File Version: 0.191 +## File Version: 0.194 #-- information criteria xxirt -xxirt_ic <- function( dev, N, par_items, par_Theta, I, par_items_bounds, np_item=NULL ) +xxirt_ic <- function( dev, N, par_items, par_Theta, I, par_items_bounds, np_item=NULL, + estimator="ML") { # Information criteria ic <- list( deviance=dev, n=N, I=I ) @@ -12,7 +13,7 @@ xxirt_ic <- function( dev, N, par_items, par_Theta, I, par_items_bounds, np_item ic$np.items <- np_item } ic$np.Theta <- length(par_Theta) - ic <- xxirt_ic_compute_criteria(ic=ic) + ic <- xxirt_ic_compute_criteria(ic=ic, estimator=estimator) return(ic) } diff --git a/R/xxirt_ic_compute_criteria.R b/R/xxirt_ic_compute_criteria.R index 4d32a66..acaa0e4 100644 --- a/R/xxirt_ic_compute_criteria.R +++ b/R/xxirt_ic_compute_criteria.R @@ -1,16 +1,18 @@ ## File Name: xxirt_ic_compute_criteria.R -## File Version: 0.05 +## File Version: 0.07 -xxirt_ic_compute_criteria <- function(ic, compute_np=TRUE) +xxirt_ic_compute_criteria <- function(ic, compute_np=TRUE, estimator="ML") { dev <- ic$deviance if (compute_np){ ic$np <- ic$np.item + ic$np.Theta } - ic$AIC <- dev + 2*ic$np - log_n <- log(ic$n) - ic$BIC <- dev + log_n * ic$np - ic$CAIC <- dev + ( log_n + 1 )*ic$np - ic$AICc <- ic$AIC + 2*ic$np * ( ic$np + 1 ) / ( ic$n - ic$np - 1 ) + if (estimator=='ML'){ + ic$AIC <- dev + 2*ic$np + log_n <- log(ic$n) + ic$BIC <- dev + log_n * ic$np + ic$CAIC <- dev + ( log_n + 1 )*ic$np + ic$AICc <- ic$AIC + 2*ic$np * ( ic$np + 1 ) / ( ic$n - ic$np - 1 ) + } return(ic) } diff --git a/R/xxirt_irf_lca_init_par.R b/R/xxirt_irf_lca_init_par.R index 67f4030..0ff8889 100644 --- a/R/xxirt_irf_lca_init_par.R +++ b/R/xxirt_irf_lca_init_par.R @@ -1,5 +1,5 @@ ## File Name: xxirt_irf_lca_init_par.R -## File Version: 0.04 +## File Version: 0.05 xxirt_irf_lca_init_par <- function(K, ncat, random_sd=0, rg_val=1.5) @@ -7,12 +7,12 @@ xxirt_irf_lca_init_par <- function(K, ncat, random_sd=0, rg_val=1.5) v1 <- NULL for (cc in 2L:ncat){ par <- rg_val*seq( -1, 1, length=K ) - names(par) <- paste0('b',1:K, '_Cat',cc-1) + names(par) <- paste0('b',1L:K, '_Cat',cc-1) par <- par + stats::rnorm(K, sd=random_sd) v1 <- c(v1, par) } if (ncat==2){ - names(v1) <- paste0('b',1:K) + names(v1) <- paste0('b',1L:K) } return(v1) } diff --git a/R/xxirt_mstep_ThetaParameters.R b/R/xxirt_mstep_ThetaParameters.R index 33df996..3c7dd2d 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.217 +## File Version: 0.218 xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, @@ -15,7 +15,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 1L:NP){ if ( ! is.na( customTheta$prior[pp] ) ) { prior_pp <- customTheta$prior[pp] # val <- par1[ names(customTheta$prior) ] diff --git a/R/xxirt_mstep_itemParameters.R b/R/xxirt_mstep_itemParameters.R index be2be8d..492dbb5 100644 --- a/R/xxirt_mstep_itemParameters.R +++ b/R/xxirt_mstep_itemParameters.R @@ -1,5 +1,5 @@ ## File Name: xxirt_mstep_itemParameters.R -## File Version: 0.405 +## File Version: 0.406 #--- M-step item parameters @@ -43,7 +43,7 @@ xxirt_mstep_itemParameters <- function( partable, item_list, items, Theta, NP <- sum( partable$parfree==1) pen1 <- grad1 <- rep(0,NP) pen0 <- 0 - for (pp in 1:NP){ + for (pp in 1L:NP){ item_index_pp <- item_index[[pp]] xh <- x xh[pp] <- xh[pp] + h @@ -64,7 +64,7 @@ xxirt_mstep_itemParameters <- function( partable, item_list, items, Theta, #-- penalty function if (!is.null(penalty_fun_item)){ pen0 <- penalty_fun_item(x=x) - for (pp in 1:NP){ + for (pp in 1L:NP){ xh <- x xh[pp] <- xh[pp] + h pen1[pp] <- penalty_fun_item(x=xh) diff --git a/R/xxirt_mstep_itemParameters_evalPrior.R b/R/xxirt_mstep_itemParameters_evalPrior.R index 1b50fa8..7aa14eb 100644 --- a/R/xxirt_mstep_itemParameters_evalPrior.R +++ b/R/xxirt_mstep_itemParameters_evalPrior.R @@ -1,5 +1,5 @@ ## File Name: xxirt_mstep_itemParameters_evalPrior.R -## File Version: 0.153 +## File Version: 0.154 #*** evaluate prior in M-step @@ -12,7 +12,7 @@ xxirt_mstep_itemParameters_evalPrior <- function(partable, h=0) #*** evaluate prior distributions in partable pen <- rep(0,NP1) if (NP>0){ - for (pp in 1:NP){ + for (pp in 1L:NP){ if ( ! is.na( partable1[pp,'prior'] ) ) { prior_pp <- partable1[pp,'prior'] val <- partable1[pp,'value'] + h diff --git a/R/xxirt_newton_raphson.R b/R/xxirt_newton_raphson.R index e839b5c..0cfbe8b 100644 --- a/R/xxirt_newton_raphson.R +++ b/R/xxirt_newton_raphson.R @@ -1,13 +1,14 @@ ## File Name: xxirt_newton_raphson.R -## File Version: 0.216 +## File Version: 0.259 xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, - control_nr, verbose=TRUE) + control_nr, verbose=TRUE, estimator="ML") { - partable <- em_out$partable - customTheta <- em_out$customTheta - + partable <- xxirt_em_args_extract(em_args=em_args, em_out=em_out, + object='partable') + customTheta <- xxirt_em_args_extract(em_args=em_args, em_out=em_out, + object='customTheta') G <- em_args$G Theta <- em_args$Theta items <- em_args$items @@ -28,18 +29,19 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, NPT <- length(par_Theta) NP <- NPI+NPT - em_args$parindex_items <- 1:NPI + em_args$parindex_items <- 1L:NPI if (NPI==0){ em_args$parindex_items <- NULL } + em_args$parindex_Theta <- (NPI+1):(NPI+NPT) x <- par0 # define parameter tables for estimation I <- ncol(dat) partable$item_group_comp <- 0 - for (ii in 1:I){ + for (ii in 1L:I){ ind_ii <- which(partable$itemnr==ii) - partable$item_group_comp[ ind_ii ] <- 1:length(ind_ii) + partable$item_group_comp[ ind_ii ] <- 1L:length(ind_ii) } partable$item_group_comp[ partable$parfree==0 ] <- 0 t2 <- table(partable$parindex) @@ -59,7 +61,7 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, i2 <- stats::aggregate(partable$itemnr, list(partable$parindex), max ) if (NPI>0){ - free_pars_design <- data.frame( pid=1:NPI, type='item', parlabel=names(par1) ) + free_pars_design <- data.frame( pid=1L:NPI, type='item', parlabel=names(par1) ) free_pars_design$one_item <- i1[,2]==i2[,2] free_pars_design$itemnr <- ifelse(free_pars_design$one_item, i1[,2], -9 ) partable_free <- partable[ partable$parfree==1 & partable$est, ] @@ -73,11 +75,28 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, test <- TRUE test <- FALSE + + ##** select functions + if (estimator=='ML'){ + opt_fun <- xxirt_nr_optim_fun + grad_fun <- xxirt_nr_grad_fun_Rcpp + } + if (estimator=='PML'){ + opt_fun <- xxirt_nr_pml_opt_fun + res <- xxirt_nr_pml_more_arguments( em_args=em_args, + grad_fun=xxirt_nr_pml_grad_fun) + em_args <- res$em_args + grad_fun <- res$grad_fun + pml_args <- em_args$pml_args + } + + + if (test){ - ll1 <- xxirt_nr_optim_fun(x=x, em_args=em_args) - grad1 <- xxirt_nr_grad_fun_numapprox(x=x, em_args=em_args) - grad2 <- xxirt_nr_grad_fun_Rcpp(x=x, em_args=em_args) - # Revalpr_round('grad1-grad2',5) + ll1 <- opt_fun(x=x, em_args=em_args) + grad1 <- xxirt_nr_grad_fun_numapprox(x=x, em_args=em_args, opt_fun=opt_fun) + grad2 <- grad_fun(x=x, em_args=em_args) + Revalpr_round('grad1-grad2',5) } #-- optimize @@ -99,12 +118,13 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, utils::flush.console() } res_opt_nr <- sirt_optimizer(optimizer=optimizer_nr, - par=x, fn=xxirt_nr_optim_fun, - grad=xxirt_nr_grad_fun_Rcpp, hessian=FALSE, + par=x, fn=opt_fun, grad=grad_fun, hessian=FALSE, control=control_nr, em_args=em_args, lower=lower, upper=upper ) #-- collect output + res_opt_nr$opt_fun <- opt_fun + res_opt_nr$grad_fun <- grad_fun x <- res_opt_nr$par partable <- xxirt_partable_include_freeParameters( partable=em_args$partable, x=x[ em_args$parindex_items ] ) @@ -113,10 +133,16 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, x=x[ em_args$parindex_Theta ]) #-- optimization function values - opt_values <- xxirt_nr_optim_fun(x=x, em_args=em_args, output_all=TRUE) + opt_values <- opt_fun(x=x, em_args=em_args, output_all=TRUE) + + #-- item probabilities + probs_items <- xxirt_compute_itemprobs( item_list=item_list, + items=items, Theta=Theta, ncat=ncat, + partable=partable, partable_index=partable_index ) #--- output res <- list(res_opt_nr=res_opt_nr, partable=partable, customTheta=customTheta, - opt_values=opt_values) + opt_values=opt_values, probs_items=probs_items, + em_args=em_args, pml_args=pml_args) return(res) } diff --git a/R/xxirt_nr_grad_fun_R.R b/R/xxirt_nr_grad_fun_R.R index 1fd453b..c24f2fd 100644 --- a/R/xxirt_nr_grad_fun_R.R +++ b/R/xxirt_nr_grad_fun_R.R @@ -1,5 +1,5 @@ ## File Name: xxirt_nr_grad_fun_R.R -## File Version: 0.13 +## File Version: 0.141 xxirt_nr_grad_fun_R <- function(x, em_args, eps=1e-100) { @@ -29,7 +29,7 @@ xxirt_nr_grad_fun_R <- function(x, em_args, eps=1e-100) MIGC <- em_args$MIGC probs_items_der <- list() - for (mm in 1:MIGC){ + for (mm in 1L:MIGC){ free_pars_design_mm <- free_pars_design[ free_pars_design$item_group_comp==mm, ] x1 <- x[ em_args$parindex_items ] x1 <- sirt_add_increment(x=x1, pos=free_pars_design_mm$pid, value=h) @@ -43,7 +43,7 @@ xxirt_nr_grad_fun_R <- function(x, em_args, eps=1e-100) } #--- loop over item parameters - for (pp in 1:NPI){ + for (pp in 1L:NPI){ free_pars_design_pp <- free_pars_design[pp,] partable_pp <- partable[ partable$parlabel==free_pars_design_pp$parlabel, ] diff --git a/R/xxirt_nr_grad_fun_Rcpp.R b/R/xxirt_nr_grad_fun_Rcpp.R index 37b3f54..3fbd7c5 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.174 +## File Version: 0.175 xxirt_nr_grad_fun_Rcpp <- function(x, em_args, eps=1e-100) { @@ -78,7 +78,7 @@ xxirt_nr_grad_fun_Rcpp <- function(x, em_args, eps=1e-100) #*** prior distributions if (NPI > 0){ - index_items <- 1:NPI + index_items <- 1L:NPI x1 <- x[ em_args$parindex_items ] partable1 <- xxirt_partable_include_freeParameters( partable=em_args$partable, x=x1) @@ -93,7 +93,7 @@ xxirt_nr_grad_fun_Rcpp <- function(x, em_args, eps=1e-100) if (!is.null(em_args$penalty_fun_item)){ pen0 <- em_args$penalty_fun_item(x=x1) pen1 <- rep(0,NPI) - for (pp in 1:NPI){ + for (pp in 1L:NPI){ xh <- sirt_add_increment(x=x1, pos=pp, value=h) pen1[pp] <- em_args$penalty_fun_item(x=xh) } diff --git a/R/xxirt_nr_grad_fun_numapprox.R b/R/xxirt_nr_grad_fun_numapprox.R index c970503..a2dd1b9 100644 --- a/R/xxirt_nr_grad_fun_numapprox.R +++ b/R/xxirt_nr_grad_fun_numapprox.R @@ -1,15 +1,17 @@ ## File Name: xxirt_nr_grad_fun_numapprox.R -## File Version: 0.04 +## File Version: 0.08 -xxirt_nr_grad_fun_numapprox <- function(x, em_args) +xxirt_nr_grad_fun_numapprox <- function(x, em_args, opt_fun) { grad <- 0*x - ll1 <- xxirt_nr_optim_fun(x=x, em_args=em_args) + ll1 <- opt_fun(x=x, em_args=em_args) h <- em_args$h - for (pp in 1:em_args$NP){ + for (pp in 1L:em_args$NP){ x1 <- sirt_add_increment(x=x, pos=pp, value=h) - ll2 <- xxirt_nr_optim_fun(x=x1, em_args=em_args) - grad[pp] <- (ll2-ll1) / h + ll2 <- opt_fun(x=x1, em_args=em_args) + x1 <- sirt_add_increment(x=x, pos=pp, value=-h) + ll2b <- opt_fun(x=x1, em_args=em_args) + grad[pp] <- (ll2-ll2b) / (2*h) } return(grad) } diff --git a/R/xxirt_nr_grad_fun_pml_casewise.R b/R/xxirt_nr_grad_fun_pml_casewise.R new file mode 100644 index 0000000..b8df858 --- /dev/null +++ b/R/xxirt_nr_grad_fun_pml_casewise.R @@ -0,0 +1,21 @@ +## File Name: xxirt_nr_grad_fun_pml_casewise.R +## File Version: 0.05 + + +#-- case-wise optimization function for PML +xxirt_nr_grad_fun_pml_casewise <- function(x, args, opt_fun, h=1e-4) +{ + NP <- length(x) + f0 <- do.call(what=opt_fun, args=args) + N <- length(f0) + grad <- matrix(0, nrow=N, ncol=NP) + colnames(grad) <- names(x) + for (pp in 1L:NP){ + args1 <- args + args1$x <- sirt_add_increment(x=x, pos=pp, value=h) + f1 <- do.call(what=opt_fun, args=args1) + grad[,pp] <- (f1-f0)/h + } + #- output + return(grad) +} diff --git a/R/xxirt_nr_opt_fun_pml_casewise.R b/R/xxirt_nr_opt_fun_pml_casewise.R new file mode 100644 index 0000000..85dc889 --- /dev/null +++ b/R/xxirt_nr_opt_fun_pml_casewise.R @@ -0,0 +1,31 @@ +## File Name: xxirt_nr_opt_fun_pml_casewise.R +## File Version: 0.03 + + +#-- case-wise optimization function for PML +xxirt_nr_opt_fun_pml_casewise <- function(x, object, eps=1e-14) +{ + em_args <- object$em_args + pml_args <- em_args$pml_args + prior_Theta <- xxirt_compute_prior_Theta_from_x(x=x, em_args=em_args) + probs_items <- xxirt_compute_prob_item_from_x(x=x, em_args=em_args) + I <- pml_args$I + K <- pml_args$K + G <- pml_args$G + TP <- pml_args$TP + W1 <- pml_args$W1 + freq1 <- pml_args$freq1 + NI2 <- pml_args$NI2 + W2_long <- as.matrix(pml_args$W2_long) + freq2 <- pml_args$freq2 + group0 <- object$group - 1 + weights <- object$weights + dat1 <- object$dat1 + dat_resp <- object$dat_resp==1 + eval_fun_args <- list(prior_Theta=prior_Theta, probs_items=probs_items, + W1=W1, W2_long=W2_long, G=G, K=K, + I=I, TP=TP, NI2=NI2, eps=eps, group0=group0, + weights=weights, dat1=dat1, dat_resp=dat_resp) + res <- do.call( what=sirt_rcpp_xxirt_nr_pml_casewise_opt_fun, args=eval_fun_args) + return(-res$case_val) +} diff --git a/R/xxirt_nr_pml_grad_fun.R b/R/xxirt_nr_pml_grad_fun.R new file mode 100644 index 0000000..28e3e4d --- /dev/null +++ b/R/xxirt_nr_pml_grad_fun.R @@ -0,0 +1,82 @@ +## File Name: xxirt_nr_pml_grad_fun.R +## File Version: 0.129 + +xxirt_nr_pml_grad_fun <- function(x, em_args, output_all=FALSE, use_rcpp=TRUE) +{ + + Theta <- em_args$Theta + pml_args <- em_args$pml_args + h <- em_args$h + + #*** 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) + eps <- 1e-14 + + I <- pml_args$I + K <- pml_args$K + G <- pml_args$G + TP <- pml_args$TP + W1 <- pml_args$W1 + freq1 <- pml_args$freq1 + NI2 <- pml_args$NI2 + W2_long <- as.matrix(pml_args$W2_long) + freq2 <- pml_args$freq2 + + NP <- em_args$NP + NPI <- em_args$NPI + NPT <- em_args$NPT + + #- define more arguments + parindex_items <- em_args$parindex_items + parindex_Theta <- em_args$parindex_Theta + pml_args$in_par_Theta <- em_args$pml_args$in_par_Theta + + #- define gradient + eval_fun_args <- list(prior_Theta=prior_Theta, probs_items=probs_items, freq1=freq1, + freq2=freq2, W1=W1, W2_long=W2_long, G=G, K=K, + I=I, TP=TP, NI2=NI2, eps=eps ) + fun_pml_eval <- sirt_rcpp_xxirt_nr_pml_opt_fun + pml_fun0 <- do.call( what=fun_pml_eval, args=eval_fun_args) + + grad <- 0*x + eval_grad_args <- eval_fun_args + eval_grad_args$NP <- NP + eval_grad_args$val1 <- pml_fun0$val1 + eval_grad_args$val2 <- pml_fun0$val2 + + eval_grad_args$der_prior_Theta <- 0*prior_Theta + eval_grad_args$der_probs_items <- der_probs_items0 <- 0*probs_items + + for (tpp in 1L:NP){ + pp_Theta <- pml_args$in_par_Theta[tpp] + + x1 <- sirt_add_increment(x=x, pos=tpp, value=h) + if (pp_Theta){ + prior_Theta1 <- xxirt_compute_prior_Theta_from_x(x=x1, em_args=em_args) + der_prior_Theta <- ( (prior_Theta1-prior_Theta)/h ) + eval_grad_args$der_prior_Theta <- der_prior_Theta + } else { + itpp <- em_args$item_index[[tpp]] + probs_items1 <- xxirt_compute_prob_item_from_x(x=x1, em_args=em_args, + item_index=itpp) + der_probs_items <- der_probs_items0 + der_probs_items[itpp,,] <- ( (probs_items1-probs_items[itpp,,,drop=FALSE])/h ) + eval_grad_args$der_probs_items <- der_probs_items + } + + eval_grad_args$pp_Theta <- pp_Theta + eval_grad_args$index_freq1 <- em_args$pml_args$index_freq1[[tpp]] + eval_grad_args$index_freq2 <- em_args$pml_args$index_freq2[[tpp]] + + grad[tpp] <- do.call( what=sirt_rcpp_xxirt_nr_pml_grad_fun_eval, + args=eval_grad_args) + + } + + #-- output + return(grad) +} + diff --git a/R/xxirt_nr_pml_more_arguments.R b/R/xxirt_nr_pml_more_arguments.R new file mode 100644 index 0000000..f0a1e62 --- /dev/null +++ b/R/xxirt_nr_pml_more_arguments.R @@ -0,0 +1,47 @@ +## File Name: xxirt_nr_pml_more_arguments.R +## File Version: 0.129 + +xxirt_nr_pml_more_arguments <- function(em_args, grad_fun) +{ + parindex <- em_args$partable$parindex + parindex1 <- setdiff(parindex,0) + if( length(parindex1)>length(unique(parindex1))){ + grad_fun <- NULL + no_itempar_constraints <- FALSE + } else { + no_itempar_constraints <- TRUE + } + parindex_items <- em_args$parindex_items + parindex_Theta <- em_args$parindex_Theta + em_args$pml_args$in_par_Theta <- sapply(1L:em_args$NP, FUN=function(pp){ + pp %in% parindex_Theta }, + simplify=TRUE ) + #- index handling + I <- em_args$I + NPI <- em_args$NPI + NP <- em_args$NP + NI2 <- em_args$pml_args$NI2 + + index_freq2 <- index_freq1 <- as.list(1L:NP) + W2_long <- em_args$pml_args$W2_long + + for (pp in 1L:NP){ + if (pp<=NPI){ + item_pp <- em_args$item_index[[pp]] + index_freq1[[pp]] <- item_pp - 1 + v1 <- sort( union( which( W2_long[,1] %in% item_pp ), + which( W2_long[,2] %in% item_pp ) ) ) + index_freq2[[pp]] <- v1 - 1 + } else { + index_freq1[[pp]] <- (1L:I) - 1 + index_freq2[[pp]] <- (1L:NI2) - 1 + } + } + em_args$pml_args$index_freq1 <- index_freq1 + em_args$pml_args$index_freq2 <- index_freq2 + em_args$pml_args$no_itempar_constraints <- no_itempar_constraints + + #-- output + res <- list(em_args=em_args, grad_fun=grad_fun) + return(res) +} diff --git a/R/xxirt_nr_pml_opt_fun.R b/R/xxirt_nr_pml_opt_fun.R new file mode 100644 index 0000000..34439fb --- /dev/null +++ b/R/xxirt_nr_pml_opt_fun.R @@ -0,0 +1,43 @@ +## File Name: xxirt_nr_pml_opt_fun.R +## File Version: 0.127 + +xxirt_nr_pml_opt_fun <- function(x, em_args, output_all=FALSE, use_rcpp=TRUE) +{ + + Theta <- em_args$Theta + pml_args <- em_args$pml_args + + #*** 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) + eps <- 1e-14 + + I <- pml_args$I + K <- pml_args$K + G <- pml_args$G + TP <- pml_args$TP + W1 <- pml_args$W1 + freq1 <- pml_args$freq1 + NI2 <- pml_args$NI2 + W2_long <- as.matrix(pml_args$W2_long) + freq2 <- pml_args$freq2 + + eval_fun_args <- list(prior_Theta=prior_Theta, probs_items=probs_items, freq1=freq1, + freq2=freq2, W1=W1, W2_long=W2_long, G=G, K=K, + I=I, TP=TP, NI2=NI2, eps=eps ) + + if (use_rcpp){ + fun_pml_eval <- sirt_rcpp_xxirt_nr_pml_opt_fun + } else { + fun_pml_eval <- xxirt_nr_pml_opt_fun_R + } + + res <- do.call( what=fun_pml_eval, args=eval_fun_args) + res <- -res$val + + #-- output + return(res) +} + diff --git a/R/xxirt_nr_pml_opt_fun_R.R b/R/xxirt_nr_pml_opt_fun_R.R new file mode 100644 index 0000000..5170a68 --- /dev/null +++ b/R/xxirt_nr_pml_opt_fun_R.R @@ -0,0 +1,43 @@ +## File Name: xxirt_nr_pml_opt_fun_R.R +## File Version: 0.06 + + +xxirt_nr_pml_opt_fun_R <- function(prior_Theta, probs_items, freq1, freq2, W1, + W2_long, G, K, I, TP, NI2, eps ) +{ + val <- 0 + + #** first-order frequencies + for (gg in 1L:G){ + for (ii in 1L:I){ + t1 <- 0 + for (hh in 1L:K){ + p1 <- sum( prior_Theta[,gg] * probs_items[ii,hh,] ) + t1 <- t1 + freq1[ii,hh,gg]*log(p1+eps) + } + val <- val + W1[ii]*t1 + } + } + + #*** second-order frequencies + for (gg in 1L:G){ + for (nn in 1L:NI2){ + ii1 <- W2_long[nn,'item1'] + ii2 <- W2_long[nn,'item2'] + w2ii <- W2_long[nn,'w2'] + t1 <- 0 + for (hh in 1L:K){ + for (kk in 1L:K){ + p1 <- sum( prior_Theta[,gg] * probs_items[ii1,hh,] * + probs_items[ii2,kk,]) + t1 <- t1 + freq2[nn,hh,kk,gg]*log(p1+eps) + } + } + val <- val + w2ii*t1 + } + } + + #- output + res <- list(vcal=val) + return(res) +} diff --git a/R/xxirt_nr_pml_preproc_data.R b/R/xxirt_nr_pml_preproc_data.R new file mode 100644 index 0000000..d731184 --- /dev/null +++ b/R/xxirt_nr_pml_preproc_data.R @@ -0,0 +1,67 @@ +## File Name: xxirt_nr_pml_preproc_data.R +## File Version: 0.11 + + +xxirt_nr_pml_preproc_data <- function(em_args, pml_args) +{ + + I <- em_args$I + ncat <- em_args$ncat + items <- em_args$items + G <- em_args$G + group <- em_args$group + Theta <- em_args$Theta + TP <- nrow(Theta) + K <- max(ncat) + if (is.null(pml_args$W1)){ + pml_args$W1 <- rep(1/I, I) + } + if (is.null(pml_args$W2)){ + W2 <- matrix( 2/I/(I-1), nrow=I, ncol=I) + rownames(W2) <- colnames(W2) <- items + W2[ upper.tri(W2)] <- 0 + diag(W2) <- 0 + pml_args$W2 <- W2 + } + W2 <- pml_args$W2 + # convert W2 into a long format + W2_long <- data.frame( item1=rep(1L:I,each=I), item2=rep(1L:I, I), + w2=as.numeric(W2) ) + W2_long <- W2_long[ W2_long$w2 > 0, ] + NI2 <- nrow(W2_long) + dat1_resp <- em_args$dat1_resp + weights <- em_args$weights + + freq1 <- array(0, dim=c(I,K,G)) + for (gg in 1L:G){ + for (ii in 1L:I){ + for (hh in 1L:K){ + freq1[ii,hh,gg] <- sum( dat1_resp[,ii,hh] * weights * (group==gg) ) + } + } + } + freq2 <- array(0, dim=c(NI2, K, K,G) ) + nn <- 1 + for (gg in 1L:G){ + for (nn in 1L:NI2){ + for (hh in 1L:K){ + for (jj in 1L:K){ + freq2[nn,hh,jj,gg] <- sum( dat1_resp[, W2_long[nn,1], hh ] * + dat1_resp[, W2_long[nn,2], jj ] * weights * + ( group==gg) ) + } + } + } + } + pml_args$W2_long <- W2_long + pml_args$I <- I + pml_args$G <- G + pml_args$K <- K + pml_args$TP <- TP + pml_args$NI2 <- NI2 + pml_args$freq1 <- freq1 + pml_args$freq2 <- freq2 + + #-- output + return(pml_args) +} diff --git a/R/xxirt_partable_extract_freeParameters.R b/R/xxirt_partable_extract_freeParameters.R index af36fe1..0bf58de 100644 --- a/R/xxirt_partable_extract_freeParameters.R +++ b/R/xxirt_partable_extract_freeParameters.R @@ -1,5 +1,5 @@ ## File Name: xxirt_partable_extract_freeParameters.R -## File Version: 0.11 +## File Version: 0.121 xxirt_partable_extract_freeParameters <- function( partable ) diff --git a/R/xxirt_postproc_parameters.R b/R/xxirt_postproc_parameters.R index b09e017..be528e0 100644 --- a/R/xxirt_postproc_parameters.R +++ b/R/xxirt_postproc_parameters.R @@ -1,5 +1,5 @@ ## File Name: xxirt_postproc_parameters.R -## File Version: 0.237 +## File Version: 0.242 @@ -27,7 +27,6 @@ xxirt_postproc_parameters <- function( partable, customTheta, p1 <- partable[ ! duplicated(partable$item ), ] dfr <- data.frame( item=items, type=paste(p1$type), m1 ) rownames(dfr) <- NULL - #*** probabilities item parameters pi_dim <- dim(probs_items) dimnames(probs_items)[[1]] <- items @@ -38,7 +37,6 @@ xxirt_postproc_parameters <- function( partable, customTheta, p1$active <- 1 * ( p1$value > p1$lower ) p1$active <- p1$active * ( p1$value < p1$upper ) par_items_bounds <- p1 - np_item <- NULL if ( ! is.null(np_fun_item) ){ np_item <- np_fun_item(x=par_items) diff --git a/R/xxirt_sandwich_pml.R b/R/xxirt_sandwich_pml.R new file mode 100644 index 0000000..cea8628 --- /dev/null +++ b/R/xxirt_sandwich_pml.R @@ -0,0 +1,46 @@ +## File Name: xxirt_sandwich_pml.R +## File Version: 0.118 + +xxirt_sandwich_pml <- function(object, h=1e-4) +{ + requireNamespace('MASS') + weights <- object$weights + partable <- object$partable + customTheta <- object$customTheta + par1 <- xxirt_partable_extract_freeParameters( partable=partable ) + par2 <- xxirt_parTheta_extract_freeParameters( customTheta=customTheta ) + par0 <- par <- x <- c(par1, par2) + NP <- length(x) + parnames <- names(x) + + # compute gradient function + opt_fun <- xxirt_nr_opt_fun_pml_casewise + args <- list(x=x, object=object, eps=1e-14) + grad <- xxirt_nr_grad_fun_pml_casewise(x=x, args=args, opt_fun=opt_fun, h=h) + + #- compute variance matrix + G1 <- grad*sqrt(weights) + A <- crossprod(G1) + rownames(A) <- colnames(A) <- parnames + + #- compute hessian + hess <- matrix(NA, nrow=NP, ncol=NP, dimnames=list(parnames,parnames)) + + pp <- 1 + for (pp in 1L:NP){ + x1 <- sirt_add_increment(x=x, pos=pp, value=h) + grad1 <- xxirt_nr_grad_fun_pml_casewise(x=x1, args=args, opt_fun=opt_fun, h=h) + hess_pp <- (grad1-grad)/h + hess[,pp] <- colSums(hess_pp*weights) + } + hess <- ( hess + t(hess) ) / 2 + B1 <- MASS::ginv(hess) + + V <- B1 %*% A %*% t(B1) + se <- sqrt(diag(V)) + names(se) <- parnames + + #-- output + res <- list(V=V, A=A, hess=hess, se=se) + return(res) +} diff --git a/R/xxirt_simulate.R b/R/xxirt_simulate.R index 5b7deda..4577a5c 100644 --- a/R/xxirt_simulate.R +++ b/R/xxirt_simulate.R @@ -1,5 +1,5 @@ ## File Name: xxirt_simulate.R -## File Version: 0.144 +## File Version: 0.145 xxirt_simulate <- function(partable, customItems, Theta, customTheta, N=1e4, method="random") @@ -38,7 +38,7 @@ xxirt_simulate <- function(partable, customItems, Theta, customTheta, N=1e4, #---- simulate item responses I <- max(partable$itemnr) items <- unique(paste(partable$item)) - ncat <- aggregate(partable$ncat, list(partable$itemnr), max)[,2] + ncat <- stats::aggregate(partable$ncat, list(partable$itemnr), max)[,2] dat <- as.data.frame( matrix(NA, nrow=N, ncol=I) ) colnames(dat) <- items diff --git a/R/xxirt_summary_parts.R b/R/xxirt_summary_parts.R index 516ce76..c85cbf8 100644 --- a/R/xxirt_summary_parts.R +++ b/R/xxirt_summary_parts.R @@ -1,5 +1,5 @@ ## File Name: xxirt_summary_parts.R -## File Version: 0.06 +## File Version: 0.088 # xxirt_summary_parts <- function(object, digits, len_disp=66) { @@ -18,11 +18,21 @@ xxirt_summary_parts <- function(object, digits, len_disp=66) object$G, 'Group(s)', '\n') sirt_display_function(length=len_disp) - cat( 'Number of EM iterations','=', object$iter_em, '\n' ) + + est_ml <- object$estimator=='ML' + + cat( 'Estimator','=', object$estimator, '\n' ) + if (est_ml){ + cat( 'Number of EM iterations','=', object$iter_em, '\n' ) + } cat( 'Number of Newton-Raphson iterations','=', object$iter_nr, '\n' ) - cat( 'Deviance','=', round( object$deviance, 2 ), ' | ' ) - cat( 'Log Likelihood','=', round( -object$deviance/2, 2 ), '\n' ) - cat( 'Penalty function','=', round( object$pen_val, 4 ), '\n' ) + if (est_ml){ + cat( 'Deviance','=', round( object$deviance, 2 ), ' | ' ) + cat( 'Log Likelihood','=', round( -object$deviance/2, 2 ), '\n' ) + cat( 'Penalty function','=', round( object$pen_val, 4 ), '\n' ) + } else { + cat( 'Optimization function','=', round( object$res_opt_nr$objective, 2), '\n') + } cat( 'Number of persons','=', object$ic$n, '\n' ) cat( 'Number of estimated parameters','=', object$ic$np, '\n' ) @@ -30,7 +40,9 @@ xxirt_summary_parts <- function(object, digits, len_disp=66) cat( ' Number of estimated distribution parameters','=', object$ic$np.Theta, '\n\n' ) #--- information criteria - rm_summary_information_criteria(object=object) + if (est_ml){ + rm_summary_information_criteria(object=object) + } #- trait parameters sirt_display_function(length=len_disp) diff --git a/README.md b/README.md index 0fc8736..86df7db 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-40 (2024-03-19) +#### GitHub version `sirt` 4.2-51 (2024-04-12) -[![](https://img.shields.io/badge/github%20version-4.2--40-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   +[![](https://img.shields.io/badge/github%20version-4.2--51-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 6974add..5c45352 100644 --- a/docs/404.html +++ b/docs/404.html @@ -24,7 +24,7 @@ sirt - 4.2-40 + 4.2-51