diff --git a/DESCRIPTION b/DESCRIPTION index ad15315..438c3a8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sirt Type: Package Title: Supplementary Item Response Theory Models -Version: 4.2-73 -Date: 2024-09-07 13:27:24 +Version: 4.2-89 +Date: 2024-11-13 10:10:05 Author: Alexander Robitzsch [aut,cre] () Maintainer: Alexander Robitzsch Description: diff --git a/R/IRT.expectedCounts.mirt.R b/R/IRT.expectedCounts.mirt.R index 87e0508..d7eb625 100644 --- a/R/IRT.expectedCounts.mirt.R +++ b/R/IRT.expectedCounts.mirt.R @@ -1,5 +1,5 @@ ## File Name: IRT.expectedCounts.mirt.R -## File Version: 0.081 +## File Version: 0.082 # IRT.expectedCounts.mirt @@ -28,7 +28,7 @@ IRT.expectedCounts.MultipleGroupClass <- function( object, ... ) ind_group <- list() pweights <- list() type <- 'exp_counts' - for (gg in 1:G){ + for (gg in 1L:G){ object <- mirt.wrapper.posterior(mirt.obj=mobj, group=groups[gg]) if (gg==1){ theta <- object$theta.k @@ -39,10 +39,10 @@ IRT.expectedCounts.MultipleGroupClass <- function( object, ... ) ind_group[[gg]] <- object$ind_group pweights[[gg]] <- object$pweights } - dims <- dim(ll_list[[1]])[1:3] + dims <- dim(ll_list[[1]])[1L:3] ll <- array(NA, dim=c(dims, G)) ll_pw <- rep(NA, object$N_orig) - for (gg in 1:G){ + for (gg in 1L:G){ ll[,,,gg] <- ll_list[[gg]] ll_pw[ ind_group[[gg]] ] <- pweights[[gg]] } diff --git a/R/RcppExports.R b/R/RcppExports.R index 321d938..904e29f 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,5 +1,5 @@ ## File Name: RcppExports.R -## File Version: 4.002073 +## File Version: 4.002089 # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 diff --git a/R/attach.environment.sirt.R b/R/attach.environment.sirt.R index 7ee868e..881e0a9 100644 --- a/R/attach.environment.sirt.R +++ b/R/attach.environment.sirt.R @@ -1,11 +1,13 @@ ## File Name: attach.environment.sirt.R -## File Version: 0.07 -################################################## -# attach all elements of an object in an environment -.attach.environment.sirt <- function( res, envir ){ +## File Version: 0.081 + + +#--- attach all elements of an object in an environment +.attach.environment.sirt <- function( res, envir ) +{ CC <- length(res) - for (cc in 1:CC){ + for (cc in 1L:CC){ assign( names(res)[cc], res[[cc]], envir=envir ) - } - } -################################################## + } +} + diff --git a/R/automatic.recode.R b/R/automatic.recode.R index fec65e1..3aec455 100644 --- a/R/automatic.recode.R +++ b/R/automatic.recode.R @@ -1,5 +1,5 @@ ## File Name: automatic.recode.R -## File Version: 1.193 +## File Version: 1.194 #*** automatic recoding of a dataset automatic.recode <- function( data, exclude=NULL, pstart.min=.6, @@ -12,7 +12,7 @@ automatic.recode <- function( data, exclude=NULL, pstart.min=.6, # compute frequencies fstart <- TAM::tam.ctt3( data, allocate=allocate, progress=FALSE) I <- ncol(data) - prbar <- floor( 10 * ( 1:I ) / (I+1) ) + prbar <- floor( 10 * ( 1L:I ) / (I+1) ) prbar <- c(1,which( diff(prbar)==1 ) ) fstart1 <- fstart diff --git a/R/brm.irf.R b/R/brm.irf.R index dad56f7..0bf305a 100644 --- a/R/brm.irf.R +++ b/R/brm.irf.R @@ -1,5 +1,5 @@ ## File Name: brm.irf.R -## File Version: 0.11 +## File Version: 0.121 #--- item response function (discretized) beta response model brm.irf <- function( Theta, delta, tau, ncat, thdim=1, eps=1E-10 ) @@ -15,7 +15,7 @@ brm.irf <- function( Theta, delta, tau, ncat, thdim=1, eps=1E-10 ) m1 <- exp( m1 / 2 ) m2 <- - Theta[,thdim] + delta + tau m2 <- exp( m2 / 2 ) - for (cc in 1:ncat){ + for (cc in 1L:ncat){ probs[,cc] <- stats::dbeta( mp[cc], shape1=m1, shape2=m2 ) } probs <- probs + eps diff --git a/R/brm.sim.R b/R/brm.sim.R index aa76394..9d03a9e 100644 --- a/R/brm.sim.R +++ b/R/brm.sim.R @@ -1,5 +1,5 @@ ## File Name: brm.sim.R -## File Version: 0.121 +## File Version: 0.123 #-- brm.sim @@ -12,7 +12,7 @@ brm.sim <- function( theta, delta, tau, K=NULL) if ( ! is.null(K) ){ br <- seq( 0, 1, len=K+1 ) } - for (ii in 1:I){ + for (ii in 1L:I){ # ii <- 1 m1 <- exp( ( theta - delta[ii] + tau[ii] ) / 2 ) n1 <- exp( ( - theta + delta[ii] + tau[ii] ) / 2 ) diff --git a/R/categorize.R b/R/categorize.R index 9a8a974..e510121 100644 --- a/R/categorize.R +++ b/R/categorize.R @@ -1,5 +1,5 @@ ## File Name: categorize.R -## File Version: 0.192 +## File Version: 0.193 #--- categorize variables into classes @@ -12,7 +12,7 @@ categorize <- function( dat, categorical=NULL, quant=NULL, lowest=0) dfr <- NULL if (! is.null( categorical) ){ VV <- length(categorical) - for (vv in 1:VV){ + for (vv in 1L:VV){ var.vv <- categorical[vv] dat.vv <- dat[,var.vv] vals.vv <- sort( unique( dat.vv ) ) @@ -31,7 +31,7 @@ categorize <- function( dat, categorical=NULL, quant=NULL, lowest=0) vars <- names(quant) VV <- length(vars) - for (vv in 1:VV){ + for (vv in 1L:VV){ vars.vv <- vars[vv] q1 <- quant[ vars.vv ] quant.vv <- stats::quantile( dat[,vars.vv], na.rm=TRUE, diff --git a/R/conf.detect.R b/R/conf.detect.R index ec38487..f2f9e03 100644 --- a/R/conf.detect.R +++ b/R/conf.detect.R @@ -1,5 +1,5 @@ ## File Name: conf.detect.R -## File Version: 1.211 +## File Version: 1.212 # Confirmatory DETECT analysis @@ -38,7 +38,7 @@ conf.detect <- function( data, score, itemcluster, bwscale=1.1, progress=TRUE, } else { ccovtable.list <- list() args_ccov_np$progress <- FALSE - for (pp in 1:PP){ + for (pp in 1L:PP){ cat( paste( 'DETECT Calculation Score ', pp, '\n', sep='') ) ; utils::flush.console() args_ccov_np$score <- score[,pp] diff --git a/R/create.ccov.R b/R/create.ccov.R index 0652b08..69b46c7 100644 --- a/R/create.ccov.R +++ b/R/create.ccov.R @@ -1,5 +1,5 @@ ## File Name: create.ccov.R -## File Version: 1.08 +## File Version: 1.091 #**** auxiliary function for creating a covariance matrix @@ -10,7 +10,7 @@ create.ccov <- function( cc, data ) ccov.matrix <- matrix( 0, nrow=I, ncol=I) rownames(ccov.matrix) <- colnames(ccov.matrix) <- colnames(data) LL <- nrow(ccc) - for (ll in 1:LL){ + for (ll in 1L:LL){ ccov.matrix[ ccc$item1ID[ll], ccc$item2ID[ll] ] <- ccc$ccov[ll] ccov.matrix[ ccc$item2ID[ll], ccc$item1ID[ll] ] <- ccov.matrix[ ccc$item1ID[ll], ccc$item2ID[ll] ] diff --git a/R/data.prep.R b/R/data.prep.R index 645be4a..0554a5e 100644 --- a/R/data.prep.R +++ b/R/data.prep.R @@ -1,5 +1,5 @@ ## File Name: data.prep.R -## File Version: 1.149 +## File Version: 1.153 #----- data preparations for rasch.jml and rasch.mml data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE, @@ -25,7 +25,7 @@ data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE, freq.patt <- apply( dat.9, 1, FUN=function(ll){ paste(ll, collapse='' ) } ) # dat1 <- data.frame( table( freq.patt ) ) } else { - freq.patt <- paste('FP', 1000000 + 1:n, sep='') + freq.patt <- paste('FP', 1000000 + 1L:n, sep='') dat1 <- data.frame( freq.patt ) colnames(dat1)[1] <- 'freq.patt' } @@ -51,8 +51,8 @@ data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE, dat2[ dat2==9 ] <- 0 # mean right dat1$mean <- rowSums( dat2 * dat2.resp ) / rowSums( dat2.resp ) - freq.patt <- data.frame( freq.patt, rowMeans( dat, na.rm=TRUE ), 1:n ) - colnames(freq.patt)[2:3] <- c('mean', 'index' ) + freq.patt <- data.frame( freq.patt, rowMeans( dat, na.rm=TRUE ), 1L:n ) + colnames(freq.patt)[2L:3] <- c('mean', 'index' ) list( dat=dat, dat2=dat2, dat2.resp=dat2.resp, dat1=dat1, freq.patt=freq.patt, I=I, n=n, dat9=dat.9 ) } @@ -66,8 +66,10 @@ data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE, .prnum <- function( matr, digits ) { VV <- ncol(matr) - for (vv in 1:VV){ - if ( is.numeric( matr[,vv]) ){ matr[,vv] <- round( matr[,vv], digits ) } + for (vv in 1L:VV){ + if ( is.numeric( matr[,vv]) ){ + matr[,vv] <- round( matr[,vv], digits ) + } } print(matr) } @@ -79,7 +81,7 @@ resp.pattern2 <- function(x) { n <- nrow(x) p <- ncol(x) - mdp <- (x %*% (2^((1:ncol(x)) - 1))) + 1 + mdp <- (x %*% (2^((1L:ncol(x)) - 1))) + 1 misspattern <- mdp[,1] misspattern <- list( miss.pattern=mdp[,1], mp.index=match( mdp[,1], sort( unique(mdp[,1] ) ) ) ) diff --git a/R/data.recode.sirt.R b/R/data.recode.sirt.R index 3f42789..cd98bf2 100644 --- a/R/data.recode.sirt.R +++ b/R/data.recode.sirt.R @@ -1,5 +1,5 @@ ## File Name: data.recode.sirt.R -## File Version: 0.01 +## File Version: 0.02 #*** utility function for recoding a raw dataset @@ -9,7 +9,7 @@ data.recode.sirt <- function( data.raw, keys ) V <- ncol(data.raw) data.scored <- matrix( 0, nrow(data.raw), ncol(data.raw) ) colnames(data.scored) <- colnames(data.raw ) - for (vv in 1:V){ + for (vv in 1L:V){ data.scored[,vv] <- 1* ( paste(data.raw[,vv])== paste(item.stat[ item.stat$item==colnames(data.raw)[vv], 'key' ]) ) data.scored[ paste( data.raw[,vv] )=='NA', vv ] <- NA diff --git a/R/data.wide2long.R b/R/data.wide2long.R index b4b5d46..2ebabf6 100644 --- a/R/data.wide2long.R +++ b/R/data.wide2long.R @@ -1,5 +1,5 @@ ## File Name: data.wide2long.R -## File Version: 0.261 +## File Version: 0.263 #--- converts a data frame in wide format into long format @@ -10,27 +10,27 @@ data.wide2long <- function( dat, id=NULL, X=NULL, Q=NULL) N <- nrow(dat) dat1 <- matrix( t(dat[,items]), nrow=N*I, ncol=1, byrow=FALSE ) dat2 <- data.frame( dat1 ) - colnames(dat2) <- "resp" - dat1 <- data.frame( id_index=rep(1:N, each=I ) ) + colnames(dat2) <- 'resp' + dat1 <- data.frame( id_index=rep(1L:N, each=I ) ) if ( ! is.null(id) ){ dat1 <- cbind( dat1, rep( dat[, id], each=I ) ) colnames(dat1)[2] <- id } - dat1 <- cbind( dat1, item=rep( items, N ), item_index=rep(1:I, N), + dat1 <- cbind( dat1, item=rep( items, N ), item_index=rep(1L:I, N), resp=dat2$resp ) if ( ! is.null(X) ){ - dat1 <- cbind( dat1, X[ rep(1:N, each=I ), ] ) + dat1 <- cbind( dat1, X[ rep(1L:N, each=I ), ] ) } rownames(dat1) <- NULL if ( ! is.null(Q) ){ if ( is.null(colnames(Q) ) ){ - colnames(Q) <- paste0("q",1:ncol(Q) ) + colnames(Q) <- paste0('q',1L:ncol(Q) ) } - if ( sum( colnames(Q) %in% "item" )==0 ){ + if ( sum( colnames(Q) %in% 'item' )==0 ){ Q <- as.data.frame(Q) Q$item <- colnames(dat) } - dat1 <- merge( x=dat1, y=Q, by="item", all.x=TRUE ) + dat1 <- merge( x=dat1, y=Q, by='item', all.x=TRUE ) } dat1 <- dat1[ order( 1E9*dat1$id_index + dat1$item_index ), ] dat1 <- data.frame( rowindex=1:(N*I), dat1 ) diff --git a/R/detect.index.R b/R/detect.index.R index bf77820..64efcff 100644 --- a/R/detect.index.R +++ b/R/detect.index.R @@ -1,5 +1,5 @@ ## File Name: detect.index.R -## File Version: 0.351 +## File Version: 0.352 @@ -24,27 +24,27 @@ detect.index <- function( ccovtable, itemcluster ) ii <- 1 indizes[ii] <- 100*mean(ccov*delta) weighted.indizes[ii] <- 100*stats::weighted.mean( ccov * delta, sqrt_N ) - parnames[ii] <- "DETECT" + parnames[ii] <- 'DETECT' #--- ASSI ii <- 2 indizes[ii] <- mean( sign_ccov * delta ) weighted.indizes[ii] <- stats::weighted.mean( sign_ccov * delta, sqrt_N ) - parnames[ii] <- "ASSI" + parnames[ii] <- 'ASSI' #--- RATIO ii <- 3 indizes[ii] <- sum( ccov * delta ) / sum( abs_ccov ) weighted.indizes[ii] <- sum( ccov * delta * sqrt_N ) / sum( abs_ccov * sqrt_N ) - parnames[ii] <- "RATIO" + parnames[ii] <- 'RATIO' #--- MADCOV ii <- 4 indizes[ii] <- 100 * mean( abs_ccov ) weighted.indizes[ii] <- 100* stats::weighted.mean( abs_ccov, sqrt_N ) - parnames[ii] <- "MADCOV100" + parnames[ii] <- 'MADCOV100' #--- MCOV ii <- 5 indizes[ii] <- 100 * mean( ccov ) weighted.indizes[ii] <- 100* stats::weighted.mean( ccov, sqrt_N ) - parnames[ii] <- "MCOV100" + parnames[ii] <- 'MCOV100' #--- output res <- data.frame( unweighted=indizes, weighted=weighted.indizes ) diff --git a/R/dif.logisticregression.R b/R/dif.logisticregression.R index 07cd6ff..cee8330 100644 --- a/R/dif.logisticregression.R +++ b/R/dif.logisticregression.R @@ -1,5 +1,5 @@ ## File Name: dif.logisticregression.R -## File Version: 1.19 +## File Version: 1.204 #---------------------------------------------------------------------------------------## # This function performs itemwise DIF analysis by using logistic regression methods ## @@ -13,14 +13,16 @@ dif.logistic.regression <- function( dat, group, score, quant=1.645) I <- ncol(dat) matr <- NULL - cat("Items ") - for (ii in 1:I){ - dat.ii <- na.omit(data.frame( "y"=dat[,ii], "score"=score, "group"=group )) - mod1 <- stats::glm( y ~ score, data=dat.ii, family="binomial") - mod2 <- stats::glm( y ~ score + group, data=dat.ii, family="binomial") - mod3 <- stats::glm( y ~ score + group + score*group, data=dat.ii, family="binomial") + cat('Items ') + for (ii in 1L:I){ + dat.ii <- na.omit(data.frame( y=dat[,ii], score=score, group=group )) + mod1 <- stats::glm( y ~ score, data=dat.ii, family='binomial') + mod2 <- stats::glm( y ~ score + group, data=dat.ii, family='binomial') + mod3 <- stats::glm( y ~ score + group + score*group, data=dat.ii, + family='binomial') - h1 <- data.frame( "item"=colnames(dat)[ii], N=sum( 1- is.na( dat[,ii] ), na.rm=TRUE)) + h1 <- data.frame( item=colnames(dat)[ii], N=sum( 1- is.na( dat[,ii] ), + na.rm=TRUE)) h1$R <- min(group) h1$F <- max(group) h1$nR <- sum( ( 1- is.na( dat[,ii] ) )* (1-group), na.rm=T) @@ -34,54 +36,54 @@ dif.logistic.regression <- function( dat, group, score, quant=1.645) h1$uniformDIF <- mod2$coef[3] h1$se.uniformDIF <- sqrt( diag( stats::vcov(mod2)) )[3] h1$t.uniformDIF <- mod2$coef[3] / sqrt( diag( stats::vcov(mod2) ) )[3] - h1$sig.uniformDIF <- "" + h1$sig.uniformDIF <- '' if ( ! is.na( h1$t.uniformDIF ) ){ - if ( h1$t.uniformDIF > quant ){ h1$sig.uniformDIF <- "+" } - if ( h1$t.uniformDIF < - quant ){ h1$sig.uniformDIF <- "-" } + if ( h1$t.uniformDIF > quant ){ h1$sig.uniformDIF <- '+' } + if ( h1$t.uniformDIF < - quant ){ h1$sig.uniformDIF <- '-' } } - h1$DIF.ETS <- "" + h1$DIF.ETS <- '' #**** # nonuniform DIF h1$nonuniformDIF <- mod3$coef[4] h1$se.nonuniformDIF <- sqrt( diag( stats::vcov(mod3)) )[4] h1$t.nonuniformDIF <- mod3$coef[4] / sqrt( diag( stats::vcov(mod3) ) )[4] - h1$sig.nonuniformDIF <- "" + h1$sig.nonuniformDIF <- '' if ( ! is.na( h1$t.nonuniformDIF ) ){ - if ( h1$t.nonuniformDIF > quant ){ h1$sig.nonuniformDIF <- "+" } - if ( h1$t.nonuniformDIF < - quant ){ h1$sig.nonuniformDIF <- "-" } + if ( h1$t.nonuniformDIF > quant ){ h1$sig.nonuniformDIF <- '+' } + if ( h1$t.nonuniformDIF < - quant ){ h1$sig.nonuniformDIF <- '-' } } matr <- rbind( matr, h1 ) - cat( ii, " " ) + cat( ii, ' ' ) utils::flush.console() - if ( ii %% 15==0 ){ cat("\n") } + if ( ii %% 15==0 ){ cat('\n') } } - cat("\n") + cat('\n') # include variable of adjusted p values - matr[, "pdiff.adj" ] <- matr$pR - matr$pF - mean( matr$pR - matr$pF ) - ind1 <- grep( "ETS", colnames(matr) ) + matr[, 'pdiff.adj' ] <- matr$pR - matr$pF - mean( matr$pR - matr$pF ) + ind1 <- grep( 'ETS', colnames(matr) ) #*** # DIF ETS classifiaction stat <- abs( matr$uniformDIF ) stat.low <- stat - quant * matr$se.uniformDIF - matr[,"DIF.ETS"] <- "B" + matr[,'DIF.ETS'] <- 'B' # DIF classification C ind <- which( ( stat > .64 ) & ( stat.low > .43 ) ) - if (length(ind) > 0){ matr[ind, "DIF.ETS"] <- "C" } + if (length(ind) > 0){ matr[ind, 'DIF.ETS'] <- 'C' } # DIF classification A ind <- which( ( stat < .43 ) | ( stat.low < 0 ) ) if (length(ind) > 0){ - matr[ind, "DIF.ETS"] <- "A" + matr[ind, 'DIF.ETS'] <- 'A' } - matr$DIF.ETS <- paste0( matr$DIF.ETS, ifelse( matr$uniformDIF > 0, "+", "-" )) + matr$DIF.ETS <- paste0( matr$DIF.ETS, ifelse( matr$uniformDIF > 0, '+', '-' )) #**** calculation of DIF variance dif1 <- dif.variance( dif=matr$uniformDIF, se.dif=matr$se.uniformDIF ) matr <- data.frame( matr[, seq(1,ind1)], uniform.EBDIF=dif1$eb.dif, DIF.SD=dif1$unweighted.DIFSD, matr[, seq(ind1+1, ncol(matr)) ] ) - cat( paste0("\nDIF SD=", round(dif1$unweighted.DIFSD, 3 ) ), "\n") + cat( paste0('\nDIF SD=', round(dif1$unweighted.DIFSD, 3 ) ), '\n') # sorting of the items g1 <- rank( matr$uniformDIF ) - matr <- data.frame( itemnr=1:nrow(matr), sortDIFindex=g1, matr ) + matr <- data.frame( itemnr=1L:nrow(matr), sortDIFindex=g1, matr ) return(matr) } #------------------------------------------------------------------------------ diff --git a/R/equating.rasch.jackknife.R b/R/equating.rasch.jackknife.R index 7481dc0..218a758 100644 --- a/R/equating.rasch.jackknife.R +++ b/R/equating.rasch.jackknife.R @@ -1,5 +1,5 @@ ## File Name: equating.rasch.jackknife.R -## File Version: 0.157 +## File Version: 0.158 @@ -11,7 +11,7 @@ equating.rasch.jackknife <- function( pars.data, display=TRUE, se.linkerror=FALS itemunits <- unique( pars.data[,1] ) N.units <- length( itemunits ) N.items <- nrow( pars.data ) - pars.data[,4] <- paste('I', 1:N.items,sep='') + pars.data[,4] <- paste('I', 1L:N.items,sep='') # display if (display){ cat( paste( 'Jackknife Equating Procedure (Stocking-Lord)\n', @@ -22,7 +22,7 @@ equating.rasch.jackknife <- function( pars.data, display=TRUE, se.linkerror=FALS res1 <- data.frame( unit=itemunits, shift=0, SD=0, linkerror=0) # perform jackknife - for (nn in 1:N.units){ + for (nn in 1L:N.units){ pars.data1 <- pars.data[ pars.data[,1] !=itemunits[nn], ] mod.nn <- equating.rasch( x=pars.data1[, c(4,2) ], y=pars.data1[, c(4,3) ] ) res1[ nn, 'shift' ] <- mod.nn$B.est$Stocking.Lord diff --git a/R/expl.detect.R b/R/expl.detect.R index e90d737..76cadd3 100644 --- a/R/expl.detect.R +++ b/R/expl.detect.R @@ -1,5 +1,5 @@ ## File Name: expl.detect.R -## File Version: 1.314 +## File Version: 1.317 #**** Exploratory DETECT analysis @@ -26,10 +26,10 @@ expl.detect <- function( data, score, nclusters, N.est=NULL, seed=NULL, N.est <- floor(N/2) } if (is.null(estsample)){ - estsample <- sort( sample( 1:N, floor( N.est ) ) ) + estsample <- sort( sample( 1L:N, floor( N.est ) ) ) } # validation sample - valsample <- setdiff( 1:N, estsample ) + valsample <- setdiff( 1L:N, estsample ) #--- Maximizing DETECT index # nonparametric estimation of conditional covariance @@ -44,11 +44,11 @@ expl.detect <- function( data, score, nclusters, N.est=NULL, seed=NULL, d <- stats::as.dist(cc1) fit <- stats::hclust(d, method=hclust_method) # hierarchical cluster analysis clusterfit <- fit - itemcluster <- data.frame( matrix( 0, I, nclusters ) ) + itemcluster <- data.frame( matrix( 0, nrow=I, ncol=nclusters ) ) itemcluster[,1] <- colnames(data) - colnames(itemcluster) <- c( 'item', paste( 'cluster', 2:nclusters, sep='') ) + colnames(itemcluster) <- c( 'item', paste( 'cluster', 2L:nclusters, sep='') ) detect.unweighted <- detect.weighted <- NULL - for (k in 2:nclusters){ + for (k in 2L:nclusters){ itemcluster[,k] <- stats::cutree( fit, k=k ) h1 <- detect.index( ccovtable=cc, itemcluster=itemcluster[,k] ) detect.unweighted <- rbind( detect.unweighted, h1$unweighted ) @@ -57,11 +57,11 @@ expl.detect <- function( data, score, nclusters, N.est=NULL, seed=NULL, parnames <- c( 'DETECT', 'ASSI', 'RATIO', 'MADCOV100', 'MCOV100') colnames(detect.unweighted) <- paste( parnames, '.est', sep='') colnames(detect.weighted) <- paste( parnames, '.est', sep='') - dfr1 <- data.frame( N.Cluster=2:nclusters ) + dfr1 <- data.frame( N.Cluster=2L:nclusters ) dfr1$N.items <- I dfr1$N.est <- N.est dfr1$N.val <- length(valsample) - dfr1$size.cluster <- sapply( 2:nclusters, FUN=function(tt){ + dfr1$size.cluster <- sapply( 2L:nclusters, FUN=function(tt){ paste( table( itemcluster[,tt] ), collapse='-' ) } ) detu <- data.frame( dfr1, detect.unweighted ) @@ -73,7 +73,7 @@ expl.detect <- function( data, score, nclusters, N.est=NULL, seed=NULL, ccov_np_args$score=score[valsample] cc <- do.call( what=ccov.np, args=ccov_np_args) detect.unweighted <- detect.weighted <- NULL - for (k in 2:nclusters){ + for (k in 2L:nclusters){ h1 <- detect.index( ccovtable=cc, itemcluster=itemcluster[,k] ) detect.unweighted <- rbind( detect.unweighted, h1$unweighted ) detect.weighted <- rbind( detect.weighted, h1$weighted ) @@ -83,13 +83,13 @@ expl.detect <- function( data, score, nclusters, N.est=NULL, seed=NULL, detu <- data.frame( detu, detect.unweighted ) detw <- data.frame( detw, detect.weighted ) } - rownames(detect.unweighted) <- paste0('Cl', 2:nclusters) + rownames(detect.unweighted) <- paste0('Cl', 2L:nclusters) rownames(detect.weighted) <- rownames(detect.unweighted) cat('\n\nDETECT (unweighted)\n\n') clopt <- which.max( detu$DETECT.est ) + 1 cat('Optimal Cluster Size is ', clopt, ' (Maximum of DETECT Index)\n\n' ) detu1 <- detu - for (vv in 6:ncol(detu)){ + for (vv in 6L:ncol(detu)){ detu1[,vv] <- round( detu1[,vv], 3) } print(detu1) diff --git a/R/f1d.irt.R b/R/f1d.irt.R index 710db67..9231c7a 100644 --- a/R/f1d.irt.R +++ b/R/f1d.irt.R @@ -1,5 +1,5 @@ ## File Name: f1d.irt.R -## File Version: 1.295 +## File Version: 1.296 #--- Functional unidimensional model (Ip et al., 2013) f1d.irt <- function( dat=NULL, nnormal=1000, nfactors=3, @@ -51,7 +51,7 @@ f1d.irt <- function( dat=NULL, nnormal=1000, nfactors=3, D <- ncol(A) # a_i ' theta_p Zpi <- matrix( 0, TP, I ) - for (dd in 1:D){ + for (dd in 1L:D){ Zpi <- Zpi + theta[,dd] * matrix( A[,dd], TP, I, byrow=TRUE ) } # Z_pi=a_i theta_p + d_i diff --git a/R/fit.adisop.R b/R/fit.adisop.R index 28cc34e..e1f4b32 100644 --- a/R/fit.adisop.R +++ b/R/fit.adisop.R @@ -1,5 +1,5 @@ ## File Name: fit.adisop.R -## File Version: 2.255 +## File Version: 2.256 #---- Fit ADISOP model @@ -14,11 +14,11 @@ fit.adisop <- function( freq.correct, wgt, conv=.0001, RR <- nrow(freq.correct) CC <- ncol(freq.correct) # auxiliary variables - scores <- 0:(RR-1) + scores <- 0L:(RR-1) item.psx <- colSums( freq.correct * wgt ) / colSums( wgt ) # initialization - dfr0 <- data.frame( stud.index=rep(1:RR, CC), - item.index=rep(1:CC,each=RR), + dfr0 <- data.frame( stud.index=rep(1L:RR, CC), + item.index=rep(1L:CC,each=RR), stud.p=rep( scores,CC), item.p=rep( item.psx, each=RR ), wgt2=matrix( as.matrix( wgt2 ), RR*CC, 1 ), @@ -50,8 +50,8 @@ fit.adisop <- function( freq.correct, wgt, conv=.0001, # calculate sum Y <- outer( X1, X2, '+' ) # preparation isotonic regression - dfr0 <- data.frame( stud.index=rep(1:RR, CC), - item.index=rep(1:CC,each=RR), + dfr0 <- data.frame( stud.index=rep(1L:RR, CC), + item.index=rep(1L:CC,each=RR), stud.X=rep( X1,CC), item.X=rep( X2, each=RR ), wgt2=matrix( as.matrix( wgt2 ), RR*CC, 1 ), @@ -75,8 +75,8 @@ fit.adisop <- function( freq.correct, wgt, conv=.0001, #-------------- end algorithm #**** calculate link function - dfr0 <- data.frame( stud.index=rep(1:RR, CC), - item.index=rep(1:CC,each=RR), + dfr0 <- data.frame( stud.index=rep(1L:RR, CC), + item.index=rep(1L:CC,each=RR), stud.p=rep( X1,CC), item.p=rep( X2, each=RR ), wgt=matrix( as.matrix(wgt), RR*CC, 1 ), diff --git a/R/fit.gradedresponse.R b/R/fit.gradedresponse.R index 5338fa7..1327b25 100644 --- a/R/fit.gradedresponse.R +++ b/R/fit.gradedresponse.R @@ -1,14 +1,12 @@ ## File Name: fit.gradedresponse.R -## File Version: 1.15 +## File Version: 1.162 -################################################### -# fit logistic model +#--- fit logistic model fit.gradedresponse <- function( freq.categories, SC, I, K, conv=.0001, maxit=100, progress=TRUE ) { - #************************* - if (progress){ cat("\n*******Graded Response Model***********\n") } + if (progress){ cat('\n*******Graded Response Model***********\n') } theta <- stats::qlogis( seq( .5, SC-1, len=SC ) / SC ) # item parameters b <- rep(0,I) @@ -28,23 +26,23 @@ fit.gradedresponse <- function( freq.categories, SC, I, K, max.increment <- max.increment * increment.factor^(-iter) # update theta res <- .update.theta.grm( theta, b, b.cat, freq.categories, - numdiff.parm, max.increment=max.increment) + numdiff.parm, max.increment=max.increment) theta <- res$theta ll <- res$ll prob <- res$prob # update b res <- .update.b.grm( theta, b, b.cat, freq.categories, - numdiff.parm, max.increment) + numdiff.parm, max.increment) b <- res$b b <- b - mean(b) # update b.cat res <- .update.bcat.grm( theta, b, b.cat, freq.categories, - numdiff.parm, max.increment) + numdiff.parm, max.increment) b.cat <- res$b.cat b.cat <- b.cat - mean(b.cat ) deviation <- max( abs( c( theta-theta0, b - b0, b.cat - b.cat0 )) ) if (progress){ - cat( "Iteration", iter, "- Deviation=", round( deviation, 6 ), "\n") + cat( 'Iteration', iter, '- Deviation=', round( deviation, 6 ), '\n') utils::flush.console() } iter <- iter + 1 @@ -52,9 +50,8 @@ fit.gradedresponse <- function( freq.categories, SC, I, K, #************************* llcase.grm <- ll / mean( rowSums(colSums( freq.categories ))) # output - res <- list( item.sc=b, cat.sc=b.cat, - person.sc=theta, ll=ll, llcase=llcase.grm, - prob=prob ) + res <- list( item.sc=b, cat.sc=b.cat, person.sc=theta, ll=ll, llcase=llcase.grm, + prob=prob ) return(res) } -#################################################################### + diff --git a/R/fit.gradedresponse_alg.R b/R/fit.gradedresponse_alg.R index 4a8b1d6..4e1e113 100644 --- a/R/fit.gradedresponse_alg.R +++ b/R/fit.gradedresponse_alg.R @@ -1,5 +1,5 @@ ## File Name: fit.gradedresponse_alg.R -## File Version: 1.174 +## File Version: 1.175 ##################################################################### @@ -11,7 +11,7 @@ I <- length(b) K <- length(b.cat) prob1 <- prob <- array( 1, dim=c(TP,I,K+1) ) - for (kk in 1:K){ + for (kk in 1L:K){ prob1[,,kk+1] <- stats::plogis( theta + matrix( b, nrow=TP, ncol=I, byrow=TRUE) + b.cat[kk] ) prob[,,kk] <- prob1[,,kk]-prob1[,,kk+1] diff --git a/R/gom_em_calc_theta.R b/R/gom_em_calc_theta.R index 6682e4e..20eef7b 100644 --- a/R/gom_em_calc_theta.R +++ b/R/gom_em_calc_theta.R @@ -1,5 +1,5 @@ ## File Name: gom_em_calc_theta.R -## File Version: 0.17 +## File Version: 0.181 #--- calculate theta grid @@ -9,9 +9,9 @@ gom_em_calc_theta <- function( K, problevels, eps=1e-5 ) if ( ! is.matrix(problevels) ){ PL <- length(problevels) m1 <- matrix(problevels, PL, 1 ) - for (kk in 2:K){ + for (kk in 2L:K){ NM <- nrow(m1) - m1 <- cbind( m1[ rep( 1:NM, PL), ], rep( problevels, each=NM) ) + m1 <- cbind( m1[ rep( 1L:NM, PL), ], rep( problevels, each=NM) ) m1 <- m1[ rowSums(m1) <=1, ] } } diff --git a/R/gom_em_est_b.R b/R/gom_em_est_b.R index efb831f..e690f77 100644 --- a/R/gom_em_est_b.R +++ b/R/gom_em_est_b.R @@ -1,5 +1,5 @@ ## File Name: gom_em_est_b.R -## File Version: 0.06 +## File Version: 0.071 @@ -8,9 +8,9 @@ gom_em_est_b <- function( lambda, I, K, n.ik, b, theta0.k, numdiff.parm=.001, max.increment, theta.k, msteps, mstepconv, eps=.001, progress=progress ) { h <- numdiff.parm - diffindex <- 1:I + diffindex <- 1L:I if (progress){ - cat(" M steps b parameter |") + cat(' M steps b parameter |') } an.ik <- aperm( n.ik, c(2,3,1) ) it <- 0 @@ -31,11 +31,11 @@ gom_em_est_b <- function( lambda, I, K, n.ik, b, theta0.k, numdiff.parm=.001, conv1 <- max( abs( b - b0 ) ) it <- it+1 if (progress){ - cat("-") + cat('-') } } if (progress){ - cat(" ", it, "Step(s) \n") + cat(' ', it, 'Step(s) \n') } res <- list(b=b, se.b=sqrt(-1/res$d2), ll=sum(res$ll0) ) return(res) diff --git a/R/gom_em_est_covariance.R b/R/gom_em_est_covariance.R index 7b2d93a..4ffcc70 100644 --- a/R/gom_em_est_covariance.R +++ b/R/gom_em_est_covariance.R @@ -1,5 +1,5 @@ ## File Name: gom_em_est_covariance.R -## File Version: 0.04 +## File Version: 0.05 #-- GOM: estimation of covariance @@ -17,7 +17,7 @@ gom_em_est_covariance <- function( f.qk.yi, Sigma, theta.kM, N ) theta.k.adj <- theta.k - matrix( mu, nrow=nrow(theta.k), ncol=ncol(theta.k), byrow=TRUE) D <- 2 - for (dd1 in 1:D){ + for (dd1 in 1L:D){ for (dd2 in dd1:D){ tk <- theta.k.adj[,dd1]*theta.k.adj[,dd2] h1 <- ( hwt %*% tk ) * delta.theta diff --git a/R/gom_em_est_lambda.R b/R/gom_em_est_lambda.R index 1c61231..42f0e74 100644 --- a/R/gom_em_est_lambda.R +++ b/R/gom_em_est_lambda.R @@ -1,5 +1,5 @@ ## File Name: gom_em_est_lambda.R -## File Version: 0.206 +## File Version: 0.207 # estimation of lambda parameters @@ -8,12 +8,12 @@ gom_em_est_lambda <- function( lambda, I, K, n.ik, numdiff.parm=.001, progress=TRUE, lambda_partable) { h <- numdiff.parm - diffindex <- 1:I + diffindex <- 1L:I Q0 <- 0*lambda se.lambda <- Q0 an.ik <- aperm( n.ik, c(2,3,1) ) if (progress){ - cat(" M steps lambda parameter |") + cat(' M steps lambda parameter |') } it <- 0 conv1 <- 1000 @@ -22,7 +22,7 @@ gom_em_est_lambda <- function( lambda, I, K, n.ik, numdiff.parm=.001, while( ( it < msteps ) & ( conv1 > mstepconv ) ){ lambda0 <- lambda - for (kk in 1:K){ + for (kk in 1L:K){ Q1 <- Q0 Q1[,kk] <- 1 parm_temp <- lambda0 @@ -46,18 +46,18 @@ gom_em_est_lambda <- function( lambda, I, K, n.ik, numdiff.parm=.001, increment <- matrix( increment[ lambda_partable$par_index ], nrow=I, ncol=K) lambda <- lambda + increment d2 <- matrix(d2a[ lambda_partable$par_index ], nrow=I, ncol=K) - for (kk in 1:K){ + for (kk in 1L:K){ lambda[,kk] <- sirt_squeeze(lambda[,kk], lower=eps, upper=1-eps) se.lambda[,kk] <- sqrt(-1/d2[,kk]) } conv1 <- max( abs( lambda - lambda0 ) ) it <- it+1 if (progress){ - cat("-") + cat('-') } } if (progress){ - cat(" ", it, "Step(s) \n") + cat(' ', it, 'Step(s) \n') } res <- list(lambda=lambda, se.lambda=se.lambda, ll=sum(res$ll0) ) return(res) diff --git a/R/linking.haberman.lq.R b/R/linking.haberman.lq.R index 25b193f..eb904d7 100644 --- a/R/linking.haberman.lq.R +++ b/R/linking.haberman.lq.R @@ -1,12 +1,16 @@ ## File Name: linking.haberman.lq.R -## File Version: 0.202 +## File Version: 0.255 linking.haberman.lq <- function(itempars, pow=2, eps=1e-3, a_log=TRUE, - use_nu=FALSE, est_pow=FALSE, lower_pow=.1, upper_pow=3) + use_nu=FALSE, est_pow=FALSE, lower_pow=.1, upper_pow=3, method="joint", + le=FALSE, vcov_list=NULL) { CALL <- match.call() s1 <- Sys.time() + use_pw <- method %in% c('pw1','pw2') + if ( (!use_pw) | (pow!=2) | (!a_log) ){ le <- FALSE } + #--- process item parameters res <- linking_proc_itempars(itempars=itempars) itempars <- res$itempars @@ -49,6 +53,16 @@ linking.haberman.lq <- function(itempars, pow=2, eps=1e-3, a_log=TRUE, for (ii in 1L:I){ X[ ind_items==ii, ii+G-1] <- 1 } + + if ( use_pw ){ + res <- linking_haberman_lq_pw_create_design(y=y, ind_studies=ind_studies, + ind_items=ind_items, method=method) + y <- res$y + X <- res$X + w <- res$w + des_pw_slopes <- res + } + #- fit res_optim$slopes <- mod0 <- lq_fit(y=y, X=X, w=w, pow=pow, eps=eps, eps_vec=eps_vec, est_pow=est_pow, lower_pow=lower_pow, @@ -57,18 +71,30 @@ linking.haberman.lq <- function(itempars, pow=2, eps=1e-3, a_log=TRUE, pow_slopes <- mod0$pow ind_groups <- 1L:(G-1) coef0_A <- coef0[ind_groups] - a_joint <- coef0[-c(ind_groups)] - ar <- y - X %*% coef0 + if (!use_pw){ + a_joint <- coef0[-c(ind_groups)] + } else { + a_joint <- NULL + } + if (!use_pw){ + ar <- y - X %*% coef0 + } else { + ar <- NULL + } if (a_log){ coef0_A <- exp(c(0,coef0_A)) - a_joint <- exp(a_joint) - ar <- a_joint[ind_items]*( exp(ar) - 1 ) + if (!use_pw){ + a_joint <- exp(a_joint) + ar <- a_joint[ind_items]*( exp(ar) - 1 ) + } } else { coef0_A <- 1+c(0,coef0_A) } resid <- itempars1 - resid$adif <- ar + if (!use_pw){ + resid$adif <- ar + } At <- coef0_A transf.personpars <- data.frame(study=studies, A_theta=coef0_A, se_A_theta=NA) @@ -83,28 +109,53 @@ linking.haberman.lq <- function(itempars, pow=2, eps=1e-3, a_log=TRUE, X[ ind_gg, gg-1] <- itempars1[ ind_gg,3]/coef0_A[gg] } } + if ( use_pw ){ + res <- linking_haberman_lq_pw_create_design(y=y, ind_studies=ind_studies, + ind_items=ind_items, method=method) + y <- res$y + X <- res$X + w <- res$w + } #- fit res_optim$intercepts <- mod0 <- lq_fit(y=y, X=X, w=w, pow=pow, eps=eps, eps_vec=eps_vec, est_pow=est_pow, lower_pow=lower_pow, upper_pow=upper_pow) coef0 <- mod0$coefficients + pow_intercepts <- mod0$pow coef0_B <- -coef0[ind_groups] - b_joint <- coef0[-c(ind_groups)] + + if (!use_pw){ + b_joint <- coef0[-c(ind_groups)] + } if (use_nu){ coef0_B <- -coef0_B - b_joint <- -b_joint + if (!use_pw){ + b_joint <- -b_joint + } + } + if (!use_pw){ + resid$b_resid <- y - X %*% coef0 } - resid$b_resid <- y - X %*% coef0 Bt <- coef0_B <- c(0, coef0_B) transf.personpars$B_theta <- coef0_B transf.personpars$se_B_theta <- NA rownames(transf.personpars) <- NULL - item$b <- b_joint + if (!use_pw){ + item$b <- b_joint + } + + res_vcov <- NULL + if (le){ + res_vcov <- linking_haberman_lq_pw_le(des=des_pw_slopes, res_optim=res_optim, + vcov_list=vcov_list) + } # end if le==TRUE #- include joint item parameters - resid$a_joint <- item[ind_items, 'a'] - resid$b_joint <- item[ind_items, 'b'] + if (!use_pw){ + resid$a_joint <- item[ind_items, 'a'] + resid$b_joint <- item[ind_items, 'b'] + } # transformation for item parameters transf.itempars <- data.frame( study=studies, At=1/At, se_At=NA, @@ -120,7 +171,7 @@ linking.haberman.lq <- function(itempars, pow=2, eps=1e-3, a_log=TRUE, pow=pow, eps=eps, item=item, resid=resid, description=description, converged=converged, a_log=a_log, use_nu=use_nu, est_pow=est_pow, pow_slopes=pow_slopes, pow_intercepts=pow_intercepts, - CALL=CALL, time=time) + vcov=res_vcov, method=method, CALL=CALL, time=time) class(res) <- 'linking.haberman.lq' return(res) diff --git a/R/linking_haberman_lq_pw_create_design.R b/R/linking_haberman_lq_pw_create_design.R new file mode 100644 index 0000000..db4e04b --- /dev/null +++ b/R/linking_haberman_lq_pw_create_design.R @@ -0,0 +1,53 @@ +## File Name: linking_haberman_lq_pw_create_design.R +## File Version: 0.141 + +linking_haberman_lq_pw_create_design <- function(y, ind_studies, ind_items, method) +{ + G <- max(ind_studies) + I <- max(ind_items) + + des <- data.frame(study=ind_studies, item=ind_items) + dfr <- NULL + y_long <- NULL + X <- NULL + w <- NULL + + for (gg in 1L:(G-1)){ + for (hh in (gg+1):G){ + items_sel <- intersect( des[ des$study==gg, 'item' ], + des[ des$study==hh, 'item' ] ) + if (length(items_sel)>0){ + dfr1 <- data.frame(study1=gg, study2=hh, item=items_sel) + I1 <- nrow(dfr1) + X1 <- matrix(0, I1, G-1) + if (gg>1){ + X1[,gg-1] <- 1 + } + X1[,hh-1] <- -1 + for (ii in items_sel){ + i1 <- which( ind_studies==gg & ind_items==ii ) + i2 <- which( ind_studies==hh & ind_items==ii ) + y1 <- y[i1] - y[i2] + y_long <- c(y_long, y1) + } + X <- rbind(X, X1) + dfr <- rbind(dfr, dfr1) + } + } # end hh + } # end gg + + #- define weights + w <- rep(1, nrow(X) ) + if (method=='pw2'){ + a1 <- stats::aggregate(1+0*des$item, list(des$item), sum) + a2 <- stats::aggregate(1+0*dfr$item, list(dfr$item), sum) + a1 <- a1[ a1[,1] %in% a2[,1], ] + w0 <- a1[,2] / a2[,2] + w <- w0[ match( dfr$item, a1[,1]) ] + } + + #-- output + res <- list(y=y_long, X=X, w=w, ind_studies=ind_studies, ind_items=ind_items, + design=dfr, G=G, I=I) + return(res) +} diff --git a/R/linking_haberman_lq_pw_le.R b/R/linking_haberman_lq_pw_le.R new file mode 100644 index 0000000..57a0121 --- /dev/null +++ b/R/linking_haberman_lq_pw_le.R @@ -0,0 +1,86 @@ +## File Name: linking_haberman_lq_pw_le.R +## File Version: 0.132 + +linking_haberman_lq_pw_le <- function(des, res_optim, vcov_list=NULL, symm_hess=FALSE) +{ + requireNamespace('MASS') + ind_studies <- des$ind_studies + ind_items <- des$ind_items + design <- des$design + I <- des$I + G <- des$G + par_delta <- c( exp(res_optim$slopes$coefficients), + res_optim$intercepts$coefficients ) + parnames <- c( paste0('sig',2L:G), paste0('mu',2L:G) ) + names(par_delta) <- parnames + par_gamma <- rep(NA, 2*I*G) + for (gg in 1L:G){ + v1 <- paste0( rep( paste0('I', 1L:I), each=2), '_', rep(c('a','b'), I), '_G', gg ) + names(par_gamma)[ 2*I*(gg-1)+seq(1,2*I) ] <- v1 + } + NIP <- nrow(itempars) + for (hh in 1L:NIP){ + ind <- 2*I*(ind_studies[hh]-1)+ 2*(ind_items[hh]-1)+1 + par_gamma[ ind ] <- itempars[hh,3] + par_gamma[ ind+1 ] <- itempars[hh,4] + } + + #** arrange variance matrix of item parameters + Vgamma <- linking_haberman_lq_pw_le_arrange_Vgamma( vcov_list=vcov_list, + par_gamma=par_gamma, I=I, G=G, ind_items=ind_items, + ind_studies=ind_studies ) + + #*** compute derivatives + args1 <- list(par_delta=par_delta, par_gamma=par_gamma, des=des) + + # H_delta + grad_delta <- do.call(what=linking_haberman_lq_pw_le_grad, args=args1) + # H_delta_delta + hess_delta_delta <- do.call(what=linking_haberman_lq_pw_le_hess_delta, args=args1) + if (symm_hess){ + hess_delta_delta <- ( hess_delta_delta + t(hess_delta_delta) ) / 2 + } + HDD <- MASS::ginv(X=hess_delta_delta) + + # H_delta_gamma + hess_delta_gamma <- do.call(what=linking_haberman_lq_pw_le_hess_gamma, args=args1) + + #-- compute standard errors + pn2 <- list(parnames, parnames) + A <- HDD %*% hess_delta_gamma + V_SE <- A %*% Vgamma %*% t(A) + dimnames(V_SE) <- pn2 + SE <- sqrt_diag_positive(x=V_SE) + + #-- compute linking errors + B <- crossprod( grad_delta ) + V_LE <- I/(I-1) * HDD %*% B %*% t(HDD) + dimnames(V_LE) <- pn2 + LE <- sqrt_diag_positive(x=V_LE) + + #-- compute bias-corrected linking errors + V_LEbc <- 0*V_LE + for (ii in 1L:I){ + ind_ii <- rep(2*I*((1L:G)-1),each=2) + rep( 2*(ii-1)+1L:2, G) + V_gammai <- Vgamma[ind_ii, ind_ii] + H_delta_gammai <- hess_delta_gamma[,ind_ii] + V_LEbc <- V_LEbc + H_delta_gammai %*% V_gammai %*% t(H_delta_gammai) + } + V_LEbc <- I/(I-1) * HDD %*% V_LEbc %*% t(HDD) + dimnames(V_LEbc) <- pn2 + V_LEbc <- V_LE - V_LEbc + LEbc <- sqrt_diag_positive(x=V_LEbc) + + #-- compute total errors + V_TE <- V_SE + V_LE + dimnames(V_TE) <- pn2 + TE <- sqrt_diag_positive(x=V_TE) + TEbc <- sqrt( SE^2 + LEbc^2 ) + + #-- output + res <- list(V_SE=V_SE, SE=SE, V_LE=V_LE, LE=LE, V_TE=V_TE, TE=TE, + V_LEbc=V_LEbc, LEbc=LEbc, TEbc=TEbc, + grad_delta=grad_delta, hess_delta_delta=hess_delta_delta, + hess_delta_gamma=hess_delta_gamma, I=I, G=G) + return(res) +} diff --git a/R/linking_haberman_lq_pw_le_arrange_Vgamma.R b/R/linking_haberman_lq_pw_le_arrange_Vgamma.R new file mode 100644 index 0000000..4c2f33c --- /dev/null +++ b/R/linking_haberman_lq_pw_le_arrange_Vgamma.R @@ -0,0 +1,22 @@ +## File Name: linking_haberman_lq_pw_le_arrange_Vgamma.R +## File Version: 0.04 + + +linking_haberman_lq_pw_le_arrange_Vgamma <- function(vcov_list, par_gamma, I, G, + ind_items, ind_studies) +{ + + NPG <- 2*I*G + Vgamma <- matrix(0, nrow=NPG, ncol=NPG, + dimnames=list(names(par_gamma), names(par_gamma) ) ) + if (!is.null(vcov_list)){ + for (gg in 1L:G){ + items_gg <- ind_items[ ind_studies==gg ] + Igg <- length(items_gg) + items_gg <- rep(items_gg, each=2) + ind2 <- 2*I*(gg-1) + 2*(items_gg-1)+rep(1L:2, Igg) + Vgamma[ind2, ind2] <- vcov_list[[gg]] + } + } + return(Vgamma) +} diff --git a/R/linking_haberman_lq_pw_le_grad.R b/R/linking_haberman_lq_pw_le_grad.R new file mode 100644 index 0000000..76f2a1c --- /dev/null +++ b/R/linking_haberman_lq_pw_le_grad.R @@ -0,0 +1,51 @@ +## File Name: linking_haberman_lq_pw_le_grad.R +## File Version: 0.101 + + +# define analytical gradient +linking_haberman_lq_pw_le_grad <- function(par_delta, par_gamma, des, h=1e-4) +{ + I <- des$I + G <- des$G + G1 <- G-1 + val <- matrix(0,nrow=I,ncol=2) + design <- des$design + w <- des$w + ND <- nrow(design) + coef_A <- c(1, par_delta[1L:(G-1)] ) + coef_B <- c(0, par_delta[(G-1)+1L:(G-1)] ) + NPD <- length(par_delta) + grad <- matrix(0, nrow=I, ncol=NPD) + + for (dd in 1L:ND){ + gg_dd <- design$study1[dd] + hh_dd <- design$study2[dd] + ii_dd <- design$item[dd] + + #-- item discriminations + w1 <- log( par_gamma[ 2*I*(gg_dd-1) + 2*(ii_dd-1)+1] ) + w2 <- log( par_gamma[ 2*I*(hh_dd-1) + 2*(ii_dd-1)+1] ) + e1 <- w[dd]*(w1-w2-log(coef_A[gg_dd])+log(coef_A[hh_dd]) ) + + if (gg_dd>1){ + grad[ii_dd,gg_dd-1] <- grad[ii_dd,gg_dd-1] - e1 + } + if (hh_dd>1){ + grad[ii_dd,hh_dd-1] <- grad[ii_dd,hh_dd-1] + e1 + } + + #-- item intercepts + y1 <- par_gamma[ 2*I*(gg_dd-1) + 2*(ii_dd-1)+2] * coef_A[ gg_dd ] + y2 <- par_gamma[ 2*I*(hh_dd-1) + 2*(ii_dd-1)+2] * coef_A[ hh_dd ] + e2 <- w[dd]*(y1-y2-coef_B[gg_dd]+coef_B[hh_dd]) + if (gg_dd>1){ + grad[ii_dd,G1+gg_dd-1] <- grad[ii_dd,G1+gg_dd-1] - e2 + } + if (hh_dd>1){ + grad[ii_dd,G1+hh_dd-1] <- grad[ii_dd,G1+hh_dd-1] + e2 + } + } + # grad <- grad/I + + return(grad) +} diff --git a/R/linking_haberman_lq_pw_le_hess_delta.R b/R/linking_haberman_lq_pw_le_hess_delta.R new file mode 100644 index 0000000..58e0223 --- /dev/null +++ b/R/linking_haberman_lq_pw_le_hess_delta.R @@ -0,0 +1,22 @@ +## File Name: linking_haberman_lq_pw_le_hess_delta.R +## File Version: 0.06 + + +linking_haberman_lq_pw_le_hess_delta <- function(par_delta, par_gamma, des, h=1e-4) +{ + I <- des$I + G <- des$G + NPD <- length(par_delta) + hess <- matrix(0, nrow=NPD, ncol=NPD) + rownames(hess) <- colnames(hess) <- names(par_delta) + grad_fun <- linking_haberman_lq_pw_le_grad + for (pp in 1L:NPD){ + par_delta2 <- mgsem_add_increment(x=par_delta, h=h, i1=pp) + args <- list(par_delta=par_delta2, par_gamma=par_gamma, des=des,h=h) + opt1 <- do.call(what=grad_fun, args=args) + args$par_delta <- mgsem_add_increment(x=par_delta, h=-h, i1=pp) + opt2 <- do.call(what=grad_fun, args=args) + hess[,pp] <- (colSums(opt1)-colSums(opt2))/(2*h) + } + return(hess) +} diff --git a/R/linking_haberman_lq_pw_le_hess_gamma.R b/R/linking_haberman_lq_pw_le_hess_gamma.R new file mode 100644 index 0000000..b984fa6 --- /dev/null +++ b/R/linking_haberman_lq_pw_le_hess_gamma.R @@ -0,0 +1,34 @@ +## File Name: linking_haberman_lq_pw_le_hess_gamma.R +## File Version: 0.04 + + +linking_haberman_lq_pw_le_hess_gamma <- function(par_delta, par_gamma, des, h=1e-4) +{ + I <- des$I + G <- des$G + NPD <- length(par_delta) + NPG <- length(par_gamma) + hess <- matrix(0, nrow=NPD, ncol=NPG) + rownames(hess) <- names(par_delta) + colnames(hess) <- names(par_gamma) + grad_fun <- linking_haberman_lq_pw_le_grad + pgi <- is.na(par_gamma) + if (sum(pgi)>0){ + par_gamma[ pgi ] <- 0 + } + + for (gg in 1L:G){ + for (ss in 1L:2){ + ind_gg <- 2*I*(gg-1)+seq(ss,2*I,by=2) + par_gamma2 <- mgsem_add_increment(x=par_gamma, h=h, i1=ind_gg) + args <- list(par_delta=par_delta, par_gamma=par_gamma2, des=des,h=h) + opt1 <- do.call(what=grad_fun, args=args) + args$par_gamma <- mgsem_add_increment(x=par_gamma, h=-h, i1=ind_gg) + opt2 <- do.call(what=grad_fun, args=args) + hess1 <- (opt1-opt2)/(2*h) + hess[,ind_gg] <- t(hess1) + } + } + + return(hess) +} diff --git a/R/lq_fit.R b/R/lq_fit.R index 58a981f..219cbe8 100644 --- a/R/lq_fit.R +++ b/R/lq_fit.R @@ -1,5 +1,5 @@ ## File Name: lq_fit.R -## File Version: 0.154 +## File Version: 0.163 lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL, est_pow=FALSE, optimizer="optim", eps_vec=10^seq(0,-10, by=-.5), @@ -24,11 +24,13 @@ lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL, fct_optim <- function(x, y, X, pow, eps, w, Xs=NULL) { beta <- x - # e <- y - X %*% beta - # e <- sirt_rcpp_lq_fit_matrix_mult( Z=Xs, y=y, beta=beta) - # val <- sum( w*(e^2 + eps )^(pow/2) ) - val <- sirt_rcpp_lq_fit_fct_optim(Z=Xs, y=y, beta=beta, pow=pow, w=w, + if (pow==0){ + e <- y - X %*% beta + val <- sum( e^2 / ( e^2 + eps ) ) + } else { + val <- sirt_rcpp_lq_fit_fct_optim(Z=Xs, y=y, beta=beta, pow=pow, w=w, eps=eps) + } return(val) } @@ -45,6 +47,20 @@ lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL, return(der) } + + if (pow==0){ + grad_optim <- function(x, y, X, pow, eps, w, Xs=NULL) + { + NP <- length(x) + grad <- rep(NA, NP) + e <- y - X %*% x + for (pp in 1L:NP){ + grad[pp] <- sum(-2*e*eps/(e^2+eps)^2*X[,pp]) + } + return(grad) + } + } + #- define epsilon sequence eps_vec <- sirt_define_eps_sequence(eps=eps, eps_vec=eps_vec) pow00 <- pow diff --git a/R/mcmc_WaldTest.R b/R/mcmc_WaldTest.R index 5d36700..04ee08c 100644 --- a/R/mcmc_WaldTest.R +++ b/R/mcmc_WaldTest.R @@ -1,5 +1,5 @@ ## File Name: mcmc_WaldTest.R -## File Version: 0.25 +## File Version: 0.261 #** Wald Test for a set of hypotheses @@ -20,10 +20,10 @@ mcmc_WaldTest <- function( mcmcobj, hypotheses ) eps <- 1E-10 NH <- sum( v1_svd$d > eps ) W <- t(c1) %*% v1a %*% c1 - stat <- c( "chi2"=W, "df"=NH) - stat["p"] <- 1 - stats::pchisq( W, df=stat["df"]) + stat <- c( chi2=W, df=NH) + stat['p'] <- 1 - stats::pchisq( W, df=stat['df']) res <- list( hypotheses_summary=s1, chisq_stat=stat) - class(res) <- "mcmc_WaldTest" + class(res) <- 'mcmc_WaldTest' return(res) } diff --git a/R/mcmc_summary.R b/R/mcmc_summary.R index 8cffa01..0e435d8 100644 --- a/R/mcmc_summary.R +++ b/R/mcmc_summary.R @@ -1,9 +1,9 @@ ## File Name: mcmc_summary.R -## File Version: 0.21 +## File Version: 0.223 -######################################################## -# mcmclist descriptives + +#**** mcmclist descriptives mcmc_summary <- function( mcmcobj, quantiles=c(.025,.05,.50,.95,.975) ) { summary.mcmcobj <- summary(mcmcobj, quantile=quantiles) @@ -29,7 +29,7 @@ mcmc_summary <- function( mcmcobj, quantiles=c(.025,.05,.50,.95,.975) ) to=max(dat.vv) ) MAP[ii] <- m1$x[ which( m1$y==max( m1$y) ) ] } - res <- data.frame( "MAP"=MAP, "Rhat"=Rhat ) + res <- data.frame( MAP=MAP, Rhat=Rhat ) rownames(res) <- vars smc3 <- res smc2 <- summary.mcmcobj$statistics @@ -42,11 +42,10 @@ mcmc_summary <- function( mcmcobj, quantiles=c(.025,.05,.50,.95,.975) ) apply( as.matrix(mcmcobj), 2, skewness.sirt ), statis[,c(3,4) ] ) colnames(statis)[3:4] <- c("MAD", "skewness" ) - dfr <- data.frame( "parameter"=rownames(smc3), - statis, smc3, "SERatio"=smc2[,4] / smc2[,2], - "sampSize"=nrow(as.matrix(mcmcobj)), "effSize"=effSize, - summary.mcmcobj$quantiles ) + dfr <- data.frame( parameter=rownames(smc3), statis, smc3, + SERatio=smc2[,4] / smc2[,2], sampSize=nrow(as.matrix(mcmcobj)), + effSize=effSize, summary.mcmcobj$quantiles ) rownames(dfr) <- NULL return(dfr) } -########################################### + diff --git a/R/mcmc_summary_print_information_criteria.R b/R/mcmc_summary_print_information_criteria.R index aac2ffe..676601b 100644 --- a/R/mcmc_summary_print_information_criteria.R +++ b/R/mcmc_summary_print_information_criteria.R @@ -1,16 +1,16 @@ ## File Name: mcmc_summary_print_information_criteria.R -## File Version: 0.13 +## File Version: 0.141 mcmc_summary_print_information_criteria <- function(object) { - cat( "Number of iterations","=", object$iter, "\n" ) - cat( "Number of burnin iterations","=", object$burnin, "\n\n" ) + cat( 'Number of iterations','=', object$iter, '\n' ) + cat( 'Number of burnin iterations','=', object$burnin, '\n\n' ) if ( ! is.null(object$description) ){ - cat( object$description, "\n") + cat( object$description, '\n') } - cat("-----------------------------------------------------------------\n") - cat( "Dbar","=", round( object$ic$Dbar, 2 ), "\n" )#, " | " ) - cat( "Dhat","=", round( object$ic$Dhat, 2 ), "\n" )#, " | " ) - cat( "pD","=", round( object$ic$pD, 2 ), " | pD=Dbar - Dhat \n" ) - cat( "DIC","=", round( object$ic$DIC, 2 ), " | DIC=Dhat + pD \n\n" ) + cat('-----------------------------------------------------------------\n') + cat( 'Dbar','=', round( object$ic$Dbar, 2 ), '\n' )#, ' | ' ) + cat( 'Dhat','=', round( object$ic$Dhat, 2 ), '\n' )#, ' | ' ) + cat( 'pD','=', round( object$ic$pD, 2 ), ' | pD=Dbar - Dhat \n' ) + cat( 'DIC','=', round( object$ic$DIC, 2 ), ' | DIC=Dhat + pD \n\n' ) } diff --git a/R/noharm.sirt.R b/R/noharm.sirt.R index 258dfcb..f15ff13 100644 --- a/R/noharm.sirt.R +++ b/R/noharm.sirt.R @@ -1,5 +1,5 @@ ## File Name: noharm.sirt.R -## File Version: 0.935 +## File Version: 0.940 ######################################## diff --git a/R/rm_facets_pp_mle.R b/R/rm_facets_pp_mle.R index 96122b5..b732055 100644 --- a/R/rm_facets_pp_mle.R +++ b/R/rm_facets_pp_mle.R @@ -1,5 +1,5 @@ ## File Name: rm_facets_pp_mle.R -## File Version: 0.184 +## File Version: 0.187 #*** person parameter estimation in partial credit model @@ -46,7 +46,8 @@ rm_facets_pp_mle <- function( data, a, b, theta, WLE=FALSE, theta <- ifelse( abs(theta) > maxval, sign(theta)*maxval, theta ) conv <- max( abs( theta - theta0) ) if (progress){ - cat("* Iteration", iter, ":", "maximum parameter change", round( conv, 5), "\n") + cat('* Iteration', iter, ':', 'maximum parameter change', + round( conv, 5), '\n') utils::flush.console() } iter <- iter + 1 diff --git a/R/rm_facets_pp_mle_calc_ll.R b/R/rm_facets_pp_mle_calc_ll.R index 404ab63..460a1a8 100644 --- a/R/rm_facets_pp_mle_calc_ll.R +++ b/R/rm_facets_pp_mle_calc_ll.R @@ -1,6 +1,5 @@ ## File Name: rm_facets_pp_mle_calc_ll.R -## File Version: 0.14 - +## File Version: 0.151 # calculate individual likelihood for item ii @@ -8,7 +7,7 @@ rm_facets_pp_mle_calc_ll <- function( probs, data, ii, eps=1e-20 ) { N <- nrow(data) probs <- log(probs) - m1 <- matrix(1:N, nrow=N, ncol=2) + m1 <- matrix(1L:N, nrow=N, ncol=2) m1[,2] <- data[,ii] + 1 h1 <- probs[ m1 ] h1[ is.na(h1) ] <- 0 diff --git a/R/rm_numdiff_index.R b/R/rm_numdiff_index.R index 6964c64..81a8181 100644 --- a/R/rm_numdiff_index.R +++ b/R/rm_numdiff_index.R @@ -1,5 +1,5 @@ ## File Name: rm_numdiff_index.R -## File Version: 0.24 +## File Version: 0.252 #################################################################### @@ -14,9 +14,12 @@ rm_numdiff_index <- function( pjk, pjk1, pjk2, n.ik, diffindex, # [items, categories, nodes] #--- evaluate expected likelihood - ll0 <- rm_grouped_expected_likelihood(pjk=pjk, n.ik=an.ik, diffindex=diffindex, eps=eps) - ll1 <- rm_grouped_expected_likelihood(pjk=pjk1, n.ik=an.ik, diffindex=diffindex, eps=eps) - ll2 <- rm_grouped_expected_likelihood(pjk=pjk2, n.ik=an.ik, diffindex=diffindex, eps=eps) + ll0 <- rm_grouped_expected_likelihood(pjk=pjk, n.ik=an.ik, + diffindex=diffindex, eps=eps) + ll1 <- rm_grouped_expected_likelihood(pjk=pjk1, n.ik=an.ik, + diffindex=diffindex, eps=eps) + ll2 <- rm_grouped_expected_likelihood(pjk=pjk2, n.ik=an.ik, + diffindex=diffindex, eps=eps) if (! is.null(prior)){ if (length(value) > length(ll0)){ @@ -40,7 +43,8 @@ rm_numdiff_index <- function( pjk, pjk1, pjk2, n.ik, diffindex, #-- trim increment - increment <- rm_numdiff_trim_increment( increment=increment, max.increment=max.increment, eps2=eps2 ) + increment <- rm_numdiff_trim_increment( increment=increment, + max.increment=max.increment, eps2=eps2 ) #--- output res <- list(increment=increment, d1=d1, d2=d2, ll0=ll0, h=h) return(res) diff --git a/R/rm_pcm_calcprobs.R b/R/rm_pcm_calcprobs.R index 983ce56..1f44311 100644 --- a/R/rm_pcm_calcprobs.R +++ b/R/rm_pcm_calcprobs.R @@ -1,5 +1,5 @@ ## File Name: rm_pcm_calcprobs.R -## File Version: 0.13 +## File Version: 0.141 @@ -7,17 +7,17 @@ rm_pcm_calcprobs <- function( a, b, Qmatrix, theta.k, I, K, TP ) { probs <- array( 0, dim=c(I,K+1,TP) ) # categories 0, ..., K - for (kk in 1:K){ + for (kk in 1L:K){ l0 <- matrix( - b[,kk], nrow=I,ncol=TP) l0 <- l0 + TAM::tam_outer( a * Qmatrix[, kk], theta.k ) probs[,kk+1,] <- l0 } probs <- exp(probs) probs1 <- probs[,1,] - for (kk in 2:(K+1)){ + for (kk in 2L:(K+1)){ probs1 <- probs1 + probs[,kk,] } - for (kk in 1:(K+1)){ + for (kk in 1L:(K+1)){ probs[,kk,] <- probs[,kk,] / probs1 } return(probs) diff --git a/R/rm_posterior.R b/R/rm_posterior.R index 951ed1b..973c734 100644 --- a/R/rm_posterior.R +++ b/R/rm_posterior.R @@ -1,5 +1,5 @@ ## File Name: rm_posterior.R -## File Version: 0.17 +## File Version: 0.181 ####################################################### @@ -18,7 +18,7 @@ rm_posterior <- function( dat2, dat2.resp, TP, pi.k, #--- expected counts n.ik <- array( 0, dim=c(TP, I, K+1 ) ) N.ik <- array( 0, dim=c(TP, I ) ) - for (kk in 1:(K+1) ){ + for (kk in 1L:(K+1) ){ n.ik[,,kk] <- crossprod( f.qk.yi, dat2.ind.resp[,,kk] ) N.ik <- N.ik + n.ik[,,kk] } diff --git a/R/rm_proc_create_pseudoraters.R b/R/rm_proc_create_pseudoraters.R index ab91b57..c8bcdc5 100644 --- a/R/rm_proc_create_pseudoraters.R +++ b/R/rm_proc_create_pseudoraters.R @@ -1,5 +1,5 @@ ## File Name: rm_proc_create_pseudoraters.R -## File Version: 0.12 +## File Version: 0.133 rm_proc_create_pseudoraters <- function( dat, rater, pid, reference_rater=NULL ) { @@ -15,17 +15,17 @@ rm_proc_create_pseudoraters <- function( dat, rater, pid, reference_rater=NULL ) m0 <- as.data.frame( matrix(NA, nrow=nrow(dat0), ncol=I) ) colnames(m0) <- colnames(dat0) - for (ii in 1:I){ + for (ii in 1L:I){ dat_ii <- m0 dat_ii[, ii ] <- dat0[,ii] - rater_ii <- paste0( rater0, "-", items[ii] ) + rater_ii <- paste0( rater0, '-', items[ii] ) rater <- c( paste(rater), paste(rater_ii) ) dat <- rbind( dat, dat_ii) pid <- c( pid, pid0) } if ( ! is.null(reference_rater) ){ - reference_rater <- paste0(reference_rater, "-", items ) + reference_rater <- paste0(reference_rater, '-', items ) } #-- output diff --git a/R/rm_trim_increments_mstep.R b/R/rm_trim_increments_mstep.R index 27e77bf..c67bfa5 100644 --- a/R/rm_trim_increments_mstep.R +++ b/R/rm_trim_increments_mstep.R @@ -1,10 +1,11 @@ ## File Name: rm_trim_increments_mstep.R -## File Version: 0.05 +## File Version: 0.062 rm_trim_increments_mstep <- function( parm, parm0, max.increment ) { increment <- parm - parm0 - increment <- rm_numdiff_trim_increment( increment=increment, max.increment=max.increment, eps2=1E-10) + increment <- rm_numdiff_trim_increment( increment=increment, + max.increment=max.increment, eps2=1E-10) parm <- parm0 + increment return(parm) } diff --git a/R/sqrt_diag_positive.R b/R/sqrt_diag_positive.R new file mode 100644 index 0000000..95b1aa5 --- /dev/null +++ b/R/sqrt_diag_positive.R @@ -0,0 +1,11 @@ +## File Name: sqrt_diag_positive.R +## File Version: 0.02 + +sqrt_diag_positive <- function(x) +{ + y <- diag(x) + y <- ifelse(y<0,0,y) + names(y) <- rownames(x) + res <- sqrt(y) + return(res) +} diff --git a/R/tetrachoric2.R b/R/tetrachoric2.R index 1e3814d..0f8632c 100644 --- a/R/tetrachoric2.R +++ b/R/tetrachoric2.R @@ -1,15 +1,15 @@ ## File Name: tetrachoric2.R -## File Version: 1.343 +## File Version: 1.346 tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, cor.smooth=TRUE, progress=TRUE) { # process data dat <- as.matrix(dat) - if (method=="Ol"){ + if (method=='Ol'){ res <- polychoric2( dat=dat, cor.smooth=cor.smooth, maxiter=maxit) } - if ( method !="Ol" ){ + if ( method !='Ol' ){ # calculate tau tau <- - stats::qnorm( colMeans(dat,na.rm=TRUE ) ) @@ -17,7 +17,7 @@ tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, dat[ dat.resp==0] <- 9 I <- ncol(dat) # calculate frequencies - dfr <- data.frame( item1=rep(1:I,I), item2=rep(1:I, each=I ) ) + dfr <- data.frame( item1=rep(1L:I,I), item2=rep(1L:I, each=I ) ) h1 <- t( ( dat==1 ) ) %*% ( dat==0 ) dfr$f11 <- matrix( t( ( dat==1 ) ) %*% ( dat==1 ), ncol=1, byrow=TRUE ) + .5 dfr$f10 <- matrix( t( ( dat==1 ) ) %*% ( dat==0 ), ncol=1, byrow=TRUE ) + .5 @@ -36,7 +36,7 @@ tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, dfr$qi2 <- stats::qnorm( dfr$pi2) #****************** # method of Bonett - if ( method %in% c("Bo","Di") ){ + if ( method %in% c('Bo','Di') ){ dfr$pmin <- ifelse( dfr$pi1 < dfr$pi2, dfr$pi1, dfr$pi2 ) dfr$c <- ( 1 - abs( dfr$pi1 - dfr$pi2 ) / 5 - ( 0.5 - dfr$pmin)^2 ) / 2 dfr$omega <- ( dfr$f00 * dfr$f11 ) / ( dfr$f01 * dfr$f10) @@ -44,17 +44,23 @@ tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, } #***** # method of Divgi - if ( method=="Di"){ + if ( method=='Di'){ dfr2 <- as.matrix(dfr) numdiffparm <- 1e-5 - dfr$r0 <- sirt_rcpp_tetrachoric2_rcpp( dfr=dfr2, h=numdiffparm, maxiter=maxit ) + dfr$r0 <- sirt_rcpp_tetrachoric2_rcpp( dfr=dfr2, h=numdiffparm, + maxiter=maxit ) } #****************** # method of Tucker - if ( method=="Tu"){ + if ( method=='Tu'){ # functions defined by Cengiz Zopluoglu - L <- function(r,h,k) {(1/(2*pi*sqrt(1-r^2)))*exp(-((h^2-2*h*k*r+k^2)/(2*(1-r^2))))} - S <- function(x,r,h,k) { x-(L(r,h,k)*delta) } + L <- function(r,h,k) + { + (1/(2*pi*sqrt(1-r^2)))*exp(-((h^2-2*h*k*r+k^2)/(2*(1-r^2)))) + } + S <- function(x,r,h,k) + { x-(L(r,h,k)*delta) + } A0 <- dfr$A0 <- dfr$p11 - dfr$pi1 * dfr$pi2 dfr$r0 <- delta / 2 dfr$iter <- 0 @@ -63,8 +69,8 @@ tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, dfr$conv <- 0 ii <- 0 - vars <- c("A0","r0","iter","conv") - if (progress){ cat(" 0 : " ) } + vars <- c('A0','r0','iter','conv') + if (progress){ cat(' 0 : ' ) } while( ( mean( dfr$conv) < 1 ) & ( ii < maxit ) ){ # iterations dfr0 <- dfr @@ -82,22 +88,22 @@ tetrachoric2 <- function( dat, method="Ol", delta=.007, maxit=1000000, i2 <- which( dfr$conv==1 & dfr0$conv==1 ) if (length(i2) > 0){ dfr[ i2, vars] <- dfr0[ i2, vars] } if (progress){ - if (ii %% 100==0 ){ cat(".") ; utils::flush.console() } - if (ii %% 1000==0 ){ cat("\n",ii, ": ") } + if (ii %% 100==0 ){ cat('.') ; utils::flush.console() } + if (ii %% 1000==0 ){ cat('\n',ii, ': ') } } } #****************** - cat("\n") + cat('\n') } TC <- matrix(NA, I, I ) diag(TC) <- 1 - TC[ as.matrix(dfr[, c("item1","item2") ] ) ] <- dfr$r0 - TC[ as.matrix(dfr[, c("item2","item1") ] ) ] <- dfr$r0 + TC[ as.matrix(dfr[, c('item1','item2') ] ) ] <- dfr$r0 + TC[ as.matrix(dfr[, c('item2','item1') ] ) ] <- dfr$r0 if (cor.smooth){ TC <- sirt_import_psych_cor.smooth(x=TC) } rownames(TC) <- colnames(TC) <- colnames(dat) res <- list(tau=tau, rho=TC ) - } # method !="Ol" + } # method !='Ol' return(res) } diff --git a/R/truescore.irt.R b/R/truescore.irt.R index 4060d25..621e90e 100644 --- a/R/truescore.irt.R +++ b/R/truescore.irt.R @@ -1,5 +1,5 @@ ## File Name: truescore.irt.R -## File Version: 0.251 +## File Version: 0.252 #--- true score item response function @@ -13,7 +13,7 @@ truescore.irt <- function( A, B, c=NULL, d=NULL, theta=seq(-3,3,len=21), I <- nrow(B) if ( is.null(c) ){ c <- rep(0,I ) } if ( is.null(d) ){ d <- rep(1,I ) } - if ( is.null(pid) ){ pid <- 1:(length(theta)) } + if ( is.null(pid) ){ pid <- 1L:(length(theta)) } # true score truescore <- truescore_irt_irf( A=A, B=B, c=c, d=d, theta=theta ) # calculate standard error of true score diff --git a/R/truescore_irt_irf.R b/R/truescore_irt_irf.R index c09cd0b..5a641c9 100644 --- a/R/truescore_irt_irf.R +++ b/R/truescore_irt_irf.R @@ -1,5 +1,5 @@ ## File Name: truescore_irt_irf.R -## File Version: 0.03 +## File Version: 0.04 #---- true score IRT - item response function @@ -11,8 +11,8 @@ truescore_irt_irf <- function( A, B, c, d, theta ) nB <- ncol(B) maxK <- nB - rowSums( is.na( B )) I <- nrow(B) - scoreM <- matrix( 1:max(maxK), nrow=TP, ncol=max(maxK), byrow=TRUE ) - for (ii in 1:I){ + scoreM <- matrix( 1L:max(maxK), nrow=TP, ncol=max(maxK), byrow=TRUE ) + for (ii in 1L:I){ prob.ii <- matrix( NA, TP, maxK[ii] ) for (kk in seq(1, maxK[ii] ) ){ prob.ii[,kk] <- exp( theta * A[ii,kk] + B[ii,kk] ) diff --git a/R/xxirt.R b/R/xxirt.R index f892f8e..3a06782 100644 --- a/R/xxirt.R +++ b/R/xxirt.R @@ -1,5 +1,5 @@ ## File Name: xxirt.R -## File Version: 1.159 +## File Version: 1.172 #--- user specified item response model @@ -9,8 +9,8 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, mstep_reltol=1E-6, maxit_nr=0, optimizer_nr="nlminb", 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) + np_fun_item=NULL, penalty_fun_theta=NULL, np_fun_theta=NULL, + pml_args=NULL, verbose_index=NULL, cv_kfold=0, cv_maxit=10) { #*** preliminaries CALL <- match.call() @@ -110,7 +110,8 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, 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, - use_grad=use_grad, penalty_fun_item=penalty_fun_item, group=group, + use_grad=use_grad, penalty_fun_item=penalty_fun_item, + penalty_fun_theta=penalty_fun_theta, group=group, par1=par1, globconv=globconv, conv=conv, verbose2=verbose2, verbose3=FALSE, verbose_index=verbose_index) @@ -205,6 +206,7 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, 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 @@ -220,18 +222,20 @@ xxirt <- function( dat, Theta=NULL, itemtype=NULL, customItems=NULL, #-- collect parameters res <- xxirt_postproc_parameters( partable=partable, customTheta=customTheta, - items=items, probs_items=probs_items, np_fun_item=np_fun_item ) + items=items, probs_items=probs_items, np_fun_item=np_fun_item, + np_fun_theta=np_fun_theta) par_items <- res$par_items par_Theta <- res$par_Theta probs_items <- res$probs_items par_items_summary <- res$par_items_summary par_items_bounds <- res$par_items_bounds np_item <- res$np_item + np_theta <- res$np_theta #-- 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, estimator=estimator) + np_item=np_item, np_theta=np_theta, estimator=estimator) #-- compute EAP if (estimator=='ML'){ diff --git a/R/xxirt_em_algorithm.R b/R/xxirt_em_algorithm.R index 01c9d3c..3c6d4e5 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.092 +## File Version: 0.103 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, I) + item_index, h, use_grad, penalty_fun_item, penalty_fun_theta, par1, + globconv, conv, verbose_index, I) { iter <- 1 dev <- 1E100 @@ -64,10 +64,15 @@ xxirt_em_algorithm <- function(maxit, verbose1, verbose2, verbose3, disp, item_l par10 <- par1 res <- xxirt_mstep_ThetaParameters( customTheta=customTheta, G=G, eps=eps, mstep_iter=mstep_iter, N.k=N.k, par1=par1, - mstep_reltol=mstep_reltol, Theta=Theta ) + mstep_reltol=mstep_reltol, Theta=Theta, + penalty_fun_theta=penalty_fun_theta) ll2 <- res$ll2 customTheta <- res$customTheta par1 <- res$par1 + pen_theta <- res$pen_theta + + #- sum of penalty values + pen_val <- pen_val + pen_theta #*** compute deviance ll_case <- log( rowSums( post_unnorm ) ) diff --git a/R/xxirt_ic.R b/R/xxirt_ic.R index b509a2d..5241a90 100644 --- a/R/xxirt_ic.R +++ b/R/xxirt_ic.R @@ -1,10 +1,10 @@ ## File Name: xxirt_ic.R -## File Version: 0.194 +## File Version: 0.195 #-- information criteria xxirt xxirt_ic <- function( dev, N, par_items, par_Theta, I, par_items_bounds, np_item=NULL, - estimator="ML") + np_theta=NULL, estimator="ML") { # Information criteria ic <- list( deviance=dev, n=N, I=I ) @@ -12,7 +12,11 @@ xxirt_ic <- function( dev, N, par_items, par_Theta, I, par_items_bounds, np_item if ( ! is.null(np_item) ){ ic$np.items <- np_item } - ic$np.Theta <- length(par_Theta) + NPT <- length(par_Theta) + if (!is.null(np_theta)){ + NPT <- np_theta + } + ic$np.Theta <- NPT ic <- xxirt_ic_compute_criteria(ic=ic, estimator=estimator) return(ic) } diff --git a/R/xxirt_mstep_ThetaParameters.R b/R/xxirt_mstep_ThetaParameters.R index 3c7dd2d..30fea40 100644 --- a/R/xxirt_mstep_ThetaParameters.R +++ b/R/xxirt_mstep_ThetaParameters.R @@ -1,9 +1,9 @@ ## File Name: xxirt_mstep_ThetaParameters.R -## File Version: 0.218 +## File Version: 0.227 xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, - mstep_iter, N.k, par1, mstep_reltol, Theta ) + mstep_iter, N.k, par1, mstep_reltol, Theta, penalty_fun_theta=NULL ) { like_Theta <- function( x, ... ) { @@ -28,6 +28,10 @@ xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, } # end pp } ll2 <- ll2 + pen + if (!is.null(penalty_fun_theta)){ + pen <- penalty_fun_theta(x=x) + ll2 <- ll2+pen + } return(2*ll2) } #----- end definition likelihood function @@ -57,6 +61,13 @@ xxirt_mstep_ThetaParameters <- function( customTheta, G, eps, ll2 <- mod$value customTheta$par[ customTheta$est ] <- par1 par1 <- xxirt_ThetaDistribution_extract_freeParameters( customTheta=customTheta ) - res <- list( par1=par1, ll2=ll2, customTheta=customTheta) + + pen_theta <- 0 + if (!is.null(penalty_fun_theta)){ + pen_theta <- penalty_fun_theta(x=par1) + } + + #--- output + res <- list( par1=par1, ll2=ll2, customTheta=customTheta, pen_theta=pen_theta) return(res) } diff --git a/R/xxirt_newton_raphson.R b/R/xxirt_newton_raphson.R index cc3a313..480a965 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.261 +## File Version: 0.262 xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, @@ -90,8 +90,6 @@ xxirt_newton_raphson <- function(em_out, em_args, maxit_nr, optimizer_nr, pml_args <- em_args$pml_args } - - if (test){ 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) diff --git a/R/xxirt_nr_optim_fun.R b/R/xxirt_nr_optim_fun.R index 75e02d2..b09695d 100644 --- a/R/xxirt_nr_optim_fun.R +++ b/R/xxirt_nr_optim_fun.R @@ -1,5 +1,5 @@ ## File Name: xxirt_nr_optim_fun.R -## File Version: 0.124 +## File Version: 0.132 xxirt_nr_optim_fun <- function(x, em_args, output_all=FALSE) { @@ -31,6 +31,14 @@ xxirt_nr_optim_fun <- function(x, em_args, output_all=FALSE) pen_val <- em_args$penalty_fun_item(x=x1) ll <- ll + pen_val } + + x2 <- x[ em_args$parindex_Theta ] + if (!is.null(em_args$penalty_fun_theta)){ + pen2 <- em_args$penalty_fun_theta(x=x2) + ll <- ll + pen2 + pen_val <- pen_val + pen2 + } + res <- ll if (output_all){ res <- list(opt_fun=ll, pen_val=pen_val, par_prior=par_prior, ll=ll0) diff --git a/R/xxirt_postproc_parameters.R b/R/xxirt_postproc_parameters.R index be528e0..c5fb1df 100644 --- a/R/xxirt_postproc_parameters.R +++ b/R/xxirt_postproc_parameters.R @@ -1,10 +1,9 @@ ## File Name: xxirt_postproc_parameters.R -## File Version: 0.242 - +## File Version: 0.249 xxirt_postproc_parameters <- function( partable, customTheta, - items, probs_items, np_fun_item=NULL ) + items, probs_items, np_fun_item=NULL, np_fun_theta=NULL ) { #**** item parameters p1 <- partable[ partable$parfree==1, ] @@ -38,14 +37,19 @@ xxirt_postproc_parameters <- function( partable, customTheta, p1$active <- p1$active * ( p1$value < p1$upper ) par_items_bounds <- p1 np_item <- NULL + np_theta <- NULL if ( ! is.null(np_fun_item) ){ np_item <- np_fun_item(x=par_items) } + if ( ! is.null(np_fun_theta) ){ + np_theta <- np_fun_theta(x=customTheta$par) + } #*** output res <- list( par_items=par_items, par_Theta=par_Theta, probs_items=probs_items, par_items_summary=dfr, - par_items_bounds=par_items_bounds, np_item=np_item ) + par_items_bounds=par_items_bounds, np_item=np_item, + np_theta=np_theta) return(res) } diff --git a/README.md b/README.md index a95c2ff..836febb 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-73 (2024-09-07) +#### GitHub version `sirt` 4.2-89 (2024-11-13) -[![](https://img.shields.io/badge/github%20version-4.2--73-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   +[![](https://img.shields.io/badge/github%20version-4.2--89-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 aee2b65..4fbcb99 100644 --- a/docs/404.html +++ b/docs/404.html @@ -24,7 +24,7 @@ sirt - 4.2-73 + 4.2-89