From 30f928bfc31676bb3a10128c97e46152dec4f216 Mon Sep 17 00:00:00 2001 From: alexanderrobitzsch Date: Wed, 27 Feb 2019 22:48:11 +0100 Subject: [PATCH] 3.3-6 --- DESCRIPTION | 4 +- R/MAR_normalized.R | 10 --- R/RcppExports.R | 2 +- R/Rhat_sirt.R | 39 ++++---- R/data.wide2long.R | 14 ++- R/invariance_alignment_cfa_config.R | 11 ++- R/invariance_alignment_cfa_config_estimate.R | 12 ++- R/invgamma2.R | 23 ++--- R/linking_haberman_als_residual_weights.R | 4 +- R/mad_normalized.R | 10 +++ R/read.fwf2.R | 17 ++-- R/sirt_trim_increment.R | 10 +++ R/truescore.irt.R | 59 ++++--------- R/truescore_irt_irf.R | 29 ++++++ R/wle.rasch.R | 93 +++++++++----------- R/write.format2.R | 15 ++-- R/write.fwf2.R | 61 ++++++------- README.md | 4 +- inst/NEWS | 10 ++- man/invariance.alignment.Rd | 16 +++- src/RcppExports.cpp | 2 +- 21 files changed, 227 insertions(+), 218 deletions(-) delete mode 100644 R/MAR_normalized.R create mode 100644 R/mad_normalized.R create mode 100644 R/sirt_trim_increment.R create mode 100644 R/truescore_irt_irf.R diff --git a/DESCRIPTION b/DESCRIPTION index e81f019e..8ac5fba0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sirt Type: Package Title: Supplementary Item Response Theory Models -Version: 3.3-1 -Date: 2019-02-08 14:12:25 +Version: 3.3-6 +Date: 2019-02-27 22:42:09 Author: Alexander Robitzsch [aut,cre] () Maintainer: Alexander Robitzsch Description: diff --git a/R/MAR_normalized.R b/R/MAR_normalized.R deleted file mode 100644 index 3299e561..00000000 --- a/R/MAR_normalized.R +++ /dev/null @@ -1,10 +0,0 @@ -## File Name: MAR_normalized.R -## File Version: 0.02 - -MAR_normalized <- function(x) -{ - x <- stats::na.omit(x) - MAR <- stats::median( abs( x - stats::median(x) ) ) - MAR <- MAR/0.6745 - return(MAR) -} diff --git a/R/RcppExports.R b/R/RcppExports.R index 5ee68a19..f5f824e7 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,5 +1,5 @@ ## File Name: RcppExports.R -## File Version: 3.003001 +## File Version: 3.003006 # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 diff --git a/R/Rhat_sirt.R b/R/Rhat_sirt.R index b6cb01b0..898ed201 100644 --- a/R/Rhat_sirt.R +++ b/R/Rhat_sirt.R @@ -1,19 +1,20 @@ ## File Name: Rhat_sirt.R -## File Version: 1.08 +## File Version: 1.09 #################################### # auxiliary functions for Rhat statistic - ############################################################################ - # Code from rube package - # Source: http://www.stat.cmu.edu/~hseltman/rube/rube0.2-16/R/Rhat.R - # Inference from Iterative Simulation Using Multiple Sequences - # Author(s): Andrew Gelman and Donald B. Rubin - # Source: Statistical Science, Vol. 7, No. 4 (Nov., 1992), pp. 457-472 - # Stable URL: http://www.jstor.org/stable/2246093 - ## Matches gelman.diag() from package "coda", but not WinBUGS() "summary" component. - ## Better than gelman.diag() because multivariate stat is not bothered to be calculated -Rhat1 <- function(mat) { + ############################################################################ + # Code from rube package + # Source: http://www.stat.cmu.edu/~hseltman/rube/rube0.2-16/R/Rhat.R + # Inference from Iterative Simulation Using Multiple Sequences + # Author(s): Andrew Gelman and Donald B. Rubin + # Source: Statistical Science, Vol. 7, No. 4 (Nov., 1992), pp. 457-472 + # Stable URL: http://www.jstor.org/stable/2246093 + ## Matches gelman.diag() from package "coda", but not WinBUGS() "summary" component. + ## Better than gelman.diag() because multivariate stat is not bothered to be calculated +Rhat1 <- function(mat) +{ m <- ncol(mat) n <- nrow(mat) b <- apply(mat,2,mean) @@ -32,11 +33,11 @@ Rhat1 <- function(mat) { } - -Rhat <- function(arr) { - dm <- dim(arr) - if (length(dm)==2) return(Rhat1(arr)) - if (dm[2]==1) return(NULL) - if (dm[3]==1) return(Rhat1(arr[,,1])) - return(apply(arr,3,Rhat1)) - } +Rhat <- function(arr) +{ + dm <- dim(arr) + if (length(dm)==2) return(Rhat1(arr)) + if (dm[2]==1) return(NULL) + if (dm[3]==1) return(Rhat1(arr[,,1])) + return(apply(arr,3,Rhat1)) +} diff --git a/R/data.wide2long.R b/R/data.wide2long.R index 0e76e647..0a11aa1b 100644 --- a/R/data.wide2long.R +++ b/R/data.wide2long.R @@ -1,8 +1,8 @@ ## File Name: data.wide2long.R -## File Version: 0.24 +## File Version: 0.25 -##################################################################### -# converts a data frame in wide format into long format + +#--- converts a data frame in wide format into long format data.wide2long <- function( dat, id=NULL, X=NULL, Q=NULL) { items <- setdiff( colnames(dat), id ) @@ -12,14 +12,12 @@ data.wide2long <- function( dat, id=NULL, X=NULL, Q=NULL) dat2 <- data.frame( dat1 ) colnames(dat2) <- "resp" dat1 <- data.frame( "id_index"=rep( 1: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 ), - "resp"=dat2$resp ) + dat1 <- cbind( dat1, item=rep( items, N ), item_index=rep(1:I, N), + resp=dat2$resp ) if ( ! is.null(X) ){ dat1 <- cbind( dat1, X[ rep(1:N, each=I ), ] ) } @@ -38,4 +36,4 @@ data.wide2long <- function( dat, id=NULL, X=NULL, Q=NULL) dat1 <- data.frame( "rowindex"=1:(N*I), dat1 ) return(dat1) } -##################################################################### + diff --git a/R/invariance_alignment_cfa_config.R b/R/invariance_alignment_cfa_config.R index 2ac1b078..cd5eebdd 100644 --- a/R/invariance_alignment_cfa_config.R +++ b/R/invariance_alignment_cfa_config.R @@ -1,8 +1,8 @@ ## File Name: invariance_alignment_cfa_config.R -## File Version: 0.12 +## File Version: 0.16 -invariance_alignment_cfa_config <- function(dat, group, ...) +invariance_alignment_cfa_config <- function(dat, group, weights=NULL, ...) { groups <- unique(group) G <- length(groups) @@ -15,12 +15,17 @@ invariance_alignment_cfa_config <- function(dat, group, ...) colnames(nu) <- items lambda <- nu err_var <- nu + weights_gg <- NULL for (gg in 1:G){ dat_gg <- dat[ group==groups[gg], ] dat_gg <- dat_gg[, colMeans(is.na(dat_gg)) < 1 ] items_gg <- colnames(dat_gg) ind_gg <- match(items_gg, items) - res <- invariance_alignment_cfa_config_estimate(dat_gg=dat_gg, ...) + if (!is.null(weights)){ + weights_gg <- weights[ind_gg] + } + res <- invariance_alignment_cfa_config_estimate(dat_gg=dat_gg, + weights=weights_gg, ...) nu[gg, ind_gg] <- res$nu lambda[gg, ind_gg] <- res$lambda err_var[gg, ind_gg] <- res$err_var diff --git a/R/invariance_alignment_cfa_config_estimate.R b/R/invariance_alignment_cfa_config_estimate.R index 99187ea1..365985cb 100644 --- a/R/invariance_alignment_cfa_config_estimate.R +++ b/R/invariance_alignment_cfa_config_estimate.R @@ -1,13 +1,19 @@ ## File Name: invariance_alignment_cfa_config_estimate.R -## File Version: 0.04 +## File Version: 0.08 -invariance_alignment_cfa_config_estimate <- function(dat_gg, ...) +invariance_alignment_cfa_config_estimate <- function(dat_gg, weights_gg=NULL, ...) { I_gg <- ncol(dat_gg) items_gg <- colnames(dat_gg) lavmodel <- paste0("F=~", paste0(items_gg, collapse="+") ) + if (is.null(weights_gg)){ + weights_name <- NULL + } else { + dat_gg$weights <- weights_gg + weights_name <- "weights" + } mod <- lavaan::cfa(data=dat_gg, model=lavmodel, std.lv=TRUE, - meanstructure=TRUE, ...) + meanstructure=TRUE, sampling.weights=weights_name, ...) partable <- lavaan::parameterTable(object=mod) lambda <- partable[ partable$op=="=~", "est"] nu <- partable[ partable$op=="~1", "est"][1:I_gg] diff --git a/R/invgamma2.R b/R/invgamma2.R index f8a188f4..293b3eca 100644 --- a/R/invgamma2.R +++ b/R/invgamma2.R @@ -1,26 +1,20 @@ ## File Name: invgamma2.R -## File Version: 0.18 +## File Version: 0.19 -################################################################### -# inverse gamma distribution for variance -rinvgamma2 <- function( n, n0, var0 ){ - # INPUT: - # N ... number of random draws - # n0 ... sample size prior - # var0 ... prior variance -# res <- 1/ stats::rgamma( N, n0 / 2, n0 * var0 / 2 ) + +#---- inverse gamma distribution for variance +rinvgamma2 <- function( n, n0, var0 ) +{ res <- sirt_rinvgamma( n, shape=n0 / 2, scale=n0 * var0 / 2 ) return(res) } -##################################################################### dinvgamma2 <- function( x, n0, var0 ){ res <- sirt_dinvgamma( x, shape=n0 / 2, scale=n0 * var0 / 2 ) return(res) } -#----------------------------------------------------------------------- -# copied from MCMCpack package +#--- copied from MCMCpack package sirt_dinvgamma <- function (x, shape, scale=1) { if (shape <=0 | scale <=0) { @@ -28,12 +22,11 @@ sirt_dinvgamma <- function (x, shape, scale=1) } alpha <- shape beta <- scale - log.density <- alpha * log(beta) - lgamma(alpha) - (alpha + - 1) * log(x) - (beta/x) + log.density <- alpha * log(beta) - lgamma(alpha) - (alpha + 1) * log(x) - (beta/x) return(exp(log.density)) } sirt_rinvgamma <- function (n, shape, scale=1) { - return(1/rgamma(n=n, shape=shape, rate=scale)) + return(1/stats::rgamma(n=n, shape=shape, rate=scale)) } diff --git a/R/linking_haberman_als_residual_weights.R b/R/linking_haberman_als_residual_weights.R index 88a46152..76c7a833 100644 --- a/R/linking_haberman_als_residual_weights.R +++ b/R/linking_haberman_als_residual_weights.R @@ -1,5 +1,5 @@ ## File Name: linking_haberman_als_residual_weights.R -## File Version: 0.370 +## File Version: 0.371 linking_haberman_als_residual_weights <- function( logaj, logaAt, @@ -28,7 +28,7 @@ linking_haberman_als_residual_weights <- function( logaj, logaAt, #-- estimation BSQ or HUB if (estimation %in% c("BSQ","HUB") ){ if (cutoff==Inf ){ - k <- k_fac*MAR_normalized(x=loga_resid) + k <- k_fac*mad_normalized(x=loga_resid) k_estimate <- TRUE cutoff <- k + eps_k } diff --git a/R/mad_normalized.R b/R/mad_normalized.R new file mode 100644 index 00000000..debe1b0e --- /dev/null +++ b/R/mad_normalized.R @@ -0,0 +1,10 @@ +## File Name: mad_normalized.R +## File Version: 0.03 + +mad_normalized <- function(x) +{ + x <- stats::na.omit(x) + res <- stats::median( abs( x - stats::median(x) ) ) + res <- res/0.6745 + return(res) +} diff --git a/R/read.fwf2.R b/R/read.fwf2.R index 56adcaa7..eff376a0 100644 --- a/R/read.fwf2.R +++ b/R/read.fwf2.R @@ -1,11 +1,10 @@ ## File Name: read.fwf2.R -## File Version: 1.08 +## File Version: 1.11 - -############################################################################ -# This function reads fwf files # -read.fwf2 <- function( file, format, variables=NULL){ +# This function reads fwf files +read.fwf2 <- function( file, format, variables=NULL) +{ ff <- readLines( file ) ind.ff1 <- c( 1, cumsum(format)[- length(format) ] + 1 ) ind.ff2 <- cumsum(format) @@ -14,10 +13,10 @@ read.fwf2 <- function( file, format, variables=NULL){ dfr <- data.frame( matrix(0, nrow=n, ncol=I ) ) for (ii in 1:I){ dfr[,ii ] <- as.numeric( substring( ff, ind.ff1[ii], ind.ff2[ii] ) ) - } + } if (!is.null(variables)){ colnames(dfr) <- variables - } - return(dfr) } -############################################################################ + return(dfr) +} + diff --git a/R/sirt_trim_increment.R b/R/sirt_trim_increment.R new file mode 100644 index 00000000..8dcaa6af --- /dev/null +++ b/R/sirt_trim_increment.R @@ -0,0 +1,10 @@ +## File Name: sirt_trim_increment.R +## File Version: 0.02 + +sirt_trim_increment <- function(increment, max_increment) +{ + eps <- 1E-10 + ci <- ceiling( abs(increment) / ( max_increment + eps ) ) + change <- ifelse( abs(increment) > max_increment, increment/(2*ci), increment ) + return(change) +} diff --git a/R/truescore.irt.R b/R/truescore.irt.R index c5839c03..57755792 100644 --- a/R/truescore.irt.R +++ b/R/truescore.irt.R @@ -1,11 +1,11 @@ ## File Name: truescore.irt.R -## File Version: 0.22 +## File Version: 0.24 -############################################################### -truescore.irt <- function( A, B, c=NULL, d=NULL, - theta=seq(-3,3,len=21), error=NULL, pid=NULL, - h=.001 ){ - #********** + +#--- true score item response function +truescore.irt <- function( A, B, c=NULL, d=NULL, theta=seq(-3,3,len=21), + error=NULL, pid=NULL, h=.001 ) +{ if ( is.vector(B) ){ B <- matrix( B, ncol=1 ) } if ( is.vector(A) ){ A <- matrix( A, ncol=1 ) } nB <- ncol(B) @@ -15,18 +15,18 @@ truescore.irt <- function( A, B, c=NULL, d=NULL, if ( is.null(d) ){ d <- rep(1,I ) } if ( is.null(pid) ){ pid <- 1:(length(theta)) } # true score - truescore <- .ts.irf( A, B, c, d, theta ) + truescore <- truescore_irt_irf( A=A, B=B, c=c, d=d, theta=theta ) # calculate standard error of true score if ( is.null( error ) ){ truescore.error <- NULL } else { - truescore1 <- .ts.irf( A, B, c, d, theta + h ) + truescore1 <- truescore_irt_irf( A=A, B=B, c=c, d=d, theta=theta + h ) truescore.error <- sqrt( ( ( truescore1 - truescore ) / h )^2 ) * error - } + } percscore <- truescore / sum( maxK ) percscore.error <- NULL if ( ! is.null( error ) ){ percscore.error <- truescore.error / sum( maxK ) - } + } # optimization function values ind <- intersect( which( !is.na( theta ) ), which( !is.na(percscore) ) ) x0 <- theta[ind] @@ -39,44 +39,15 @@ truescore.irt <- function( A, B, c=NULL, d=NULL, # define fit function sum( ( y0 - minf - (maxf - minf) * stats::plogis( a*x0 + b ) )^2 ) } -# m1 <- glm( y0 ~ x0, family="binomial", control=list(maxit=4) ) -# m1c <- coef(m1) -# h1 <- optim( c(m1c[2], -m1c[1] ), critfct) h1 <- stats::optim( c( .5, 0 ), critfct ) -# fitval <- c( h1$par[1], h1$par[2] ) -# names(fitval) <- c("a", "b" ) - #***** - # OUTPUT - res <- data.frame( "pid"=pid,"truescore"=truescore ) + + #--- OUTPUT + res <- data.frame( "pid"=pid, "truescore"=truescore ) res$truescore.error <- truescore.error res2 <- data.frame( "percscore"=percscore ) res2$percscore.error <- percscore.error res3 <- data.frame( "lower"=minf, "upper"=maxf, "a"=h1$par[1], "b"=h1$par[2] ) res <- cbind( res, res2, res3 ) return(res) - } -#################################################### -.ts.irf <- function( A, B, c, d, theta ){ - TP <- length(theta) - truescore <- rep(0, length(theta)) - # extract maximum item categories from B parameters - 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){ - # ii <- 21 - prob.ii <- matrix( NA, TP, maxK[ii] ) - for (kk in seq(1, maxK[ii] ) ){ - # kk <- 1 - prob.ii[,kk] <- exp( theta * A[ii,kk] + B[ii,kk] ) - } - prob.ii <- prob.ii / ( rowSums( prob.ii ) + 1 ) - if ( ( maxK[ii]==1 ) & ( abs(c[ii])+abs(1-d[ii]) > 0 ) ){ - prob.ii <- c[ii] + ( d[ii] - c[ii] ) * prob.ii - } - truescore <- truescore + rowSums( prob.ii * scoreM[, seq( 1, maxK[ii] ), drop=FALSE] ) - } - return(truescore) - } -############################################################ +} + diff --git a/R/truescore_irt_irf.R b/R/truescore_irt_irf.R new file mode 100644 index 00000000..ce5d22ec --- /dev/null +++ b/R/truescore_irt_irf.R @@ -0,0 +1,29 @@ +## File Name: truescore_irt_irf.R +## File Version: 0.01 + + + +#---- true score IRT - item response function +truescore_irt_irf <- function( A, B, c, d, theta ) +{ + TP <- length(theta) + truescore <- rep(0, length(theta)) + # extract maximum item categories from B parameters + 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){ + 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] ) + } + prob.ii <- prob.ii / ( rowSums( prob.ii ) + 1 ) + if ( ( maxK[ii]==1 ) & ( abs(c[ii])+abs(1-d[ii]) > 0 ) ){ + prob.ii <- c[ii] + ( d[ii] - c[ii] ) * prob.ii + } + truescore <- truescore + rowSums( prob.ii * scoreM[, seq( 1, maxK[ii] ), drop=FALSE] ) + } + return(truescore) +} + diff --git a/R/wle.rasch.R b/R/wle.rasch.R index 3cf44bd2..c9618f37 100644 --- a/R/wle.rasch.R +++ b/R/wle.rasch.R @@ -1,88 +1,77 @@ ## File Name: wle.rasch.R -## File Version: 1.20 +## File Version: 1.294 -#---------------------------------------------------------- -# Function for WLE ability estimation -wle.rasch <- function( dat, dat.resp=NULL, b, # a=1 + 0*b, c=0*b, - itemweights=1+0*b, - theta=rep(0, nrow(dat)), conv=.001, maxit=200, - wle.adj=0, progress=FALSE) + +#-- WLE ability estimation +wle.rasch <- function( dat, dat.resp=NULL, b, itemweights=1+0*b, + theta=rep(0, nrow(dat)), conv=.001, maxit=200, + wle.adj=0, progress=FALSE) { - #-------------------------------------------------------# - # INPUT: # - # dat ... data frame with item response pattern # - # dat.resp ... response pattern (non-missing items) # - # theta ... initial values for ability estimation # - # itemweights itemweights for likelihood estimation # - # b ... item difficulty parameters # - # a ... item discriminations # - # c ... guessing parameter # - # conv ... convergence criterion # - #-------------------------------------------------------# theta.change <- 1 iter <- 0 -# if ( sum(theta)==0 ){ -# theta <- qnorm( (rowMeans( dat, na.rm=T ) + .01 ) / 1.02 ) -# } + N <- nrow(dat) if ( is.null(dat.resp) ){ - dat.resp <- 1 - is.na( dat) - # multiply response matrix with itemweights - dat.resp <- dat.resp * outer( rep(1,nrow(dat)), itemweights ) - } - dat[ is.na( dat) ] <- 9 + dat.resp <- 1 - is.na( dat) + # multiply response matrix with itemweights + dat.resp <- dat.resp * sirt_matrix2(itemweights, nrow=N) + } + dat[ is.na( dat) ] <- 9 sufftheta <- rowSums( dat.resp * dat ) if (wle.adj>0){ nr <- rowSums( dat.resp ) i1 <- which( nr==sufftheta ) - if ( length(i1) > 0 ){ sufftheta[ i1 ] <- nr[ i1 ] - wle.adj } + if ( length(i1) > 0 ){ + sufftheta[ i1 ] <- nr[ i1 ] - wle.adj + } i1 <- which( sufftheta==0 ) - if ( length(i1) > 0 ){ sufftheta[ i1 ] <- wle.adj } - } + if ( length(i1) > 0 ){ + sufftheta[ i1 ] <- wle.adj + } + } if ( progress){ cat("\n WLE estimation |" ) } old_increment <- rep( 5, length(theta)) + eps <- 1E-10 while( ( max(abs( theta.change )) > conv ) & ( iter < maxit ) ){ # calculate P and Q - p.ia <- stats::plogis( outer( theta, b, "-" ) ) + p.ia <- stats::plogis( theta - sirt_matrix2(b, nrow=N) ) q.ia <- 1 - p.ia # Log Likelihood (for every subject) - l1 <- sufftheta - rowSums( dat.resp * p.ia ) + M0 <- dat.resp * p.ia + l1 <- sufftheta - rowSums(M0) # I and J terms - M1 <- dat.resp * p.ia * q.ia - I.wle <- rowSums( M1 ) - I.wle <- rowSums( dat.resp * p.ia * q.ia ) - J.wle <- rowSums( M1 * (q.ia - p.ia ) ) + M1 <- M0*q.ia + I.wle <- rowSums(M1) + J.wle <- rowSums(M1*(q.ia - p.ia)) I1.wle <- J.wle - J1.wle <- rowSums( M1* ( 6 * p.ia^2 - 6 * p.ia + 1 ) ) + J1.wle <- rowSums( M1*( 6 * p.ia^2 - 6 * p.ia + 1 ) ) # calculate objective function - f.obj <- l1 + J.wle / ( 2 * I.wle ) + f.obj <- l1 + J.wle / ( 2 * I.wle + eps) # derivative of the objective function - f1.obj <- - I.wle + ( J1.wle * I.wle - J.wle * I1.wle )/ ( 2 * I.wle^2 ) + f1.obj <- - I.wle + ( J1.wle * I.wle - J.wle * I1.wle )/ ( 2 * I.wle^2 + eps) # theta change - increment <- theta.change <- - f.obj / f1.obj + theta.change <- - f.obj / ( f1.obj + eps ) # define damped increment!! - ci <- ceiling( abs(increment) / ( abs( old_increment) + 1E-10 ) ) - theta.change <- ifelse( abs( increment) > abs(old_increment), increment/(2*ci), increment ) + max_increment <- abs(old_increment) + theta.change <- sirt_trim_increment(increment=theta.change, max_increment=max_increment) old_increment <- abs(theta.change) theta <- theta + theta.change iter <- iter + 1 if ( any( is.nan( theta.change ) ) ){ - stop( "Numerical problems occur during WLE estimation procedure.") - } - if ( progress){ cat("-") } + stop( "Numerical problems occur during WLE estimation procedure.") + } + if ( progress){ + cat("-") + } } - res <- list( "theta"=theta, "se.theta"=1 / sqrt( abs( f1.obj ) ), - "dat.resp"=dat.resp, "p.ia"=p.ia ) + res <- list( theta=theta, se.theta=1/sqrt(abs(f1.obj)), + dat.resp=dat.resp, p.ia=p.ia ) #*** compute WLE reliability v1 <- stats::var(res$theta) v2 <- mean(res$se.theta^2 ) wle.rel <- ( v1 - v2 ) / v1 - cat("WLE Reliability=", round(wle.rel,3), "\n") + cat("\nWLE Reliability=", round(wle.rel,3), "\n") res$wle.rel <- wle.rel return(res) - } -#--------------------------------------------------------------------------------------------------------# - - - +} diff --git a/R/write.format2.R b/R/write.format2.R index ef089b8a..abce6402 100644 --- a/R/write.format2.R +++ b/R/write.format2.R @@ -1,16 +1,16 @@ ## File Name: write.format2.R -## File Version: 1.09 +## File Version: 1.12 -#---------------------------------------------------------------- -# utility function for formatting output in write.fwf2 -write.format2 <- function( vec1, ff, fr ){ +#*** utility function for formatting output in write.fwf2 +write.format2 <- function( vec1, ff, fr ) +{ if (fr==0){ vec2 <- round( vec1, fr ) blank.vv <- paste( rep( " ", ff ), collapse="" ) vec2 <- paste( substring( blank.vv, 1, ff - nchar(vec2) ), vec2, sep="") - } else { + } else { d.vv <- round( vec1, fr ) + 10^(-(fr+1)) # generate blank blank.vv <- paste( rep( " ", ff+1 ), collapse="" ) @@ -19,9 +19,8 @@ write.format2 <- function( vec1, ff, fr ){ d.vv[ g.vv ] <- ifelse( ff > 1, gsub( "NA", " .", d.vv[g.vv] ), gsub( "NA", ".", d.vv[g.vv] ) ) vec2 <- substring( d.vv, 1, ff ) vec2 - } - return(vec2 ) } -#--------------------------------------------------------------- + return(vec2) +} .write.format2 <- write.format2 diff --git a/R/write.fwf2.R b/R/write.fwf2.R index e853c940..6d0622de 100644 --- a/R/write.fwf2.R +++ b/R/write.fwf2.R @@ -1,39 +1,34 @@ ## File Name: write.fwf2.R -## File Version: 1.10 +## File Version: 1.11 -############################################################################## -write.fwf2 <- function( dat, format.full, format.round, savename ){ - if (is.null( colnames(dat) ) ){ - colnames(dat) <- paste( "V", 1:( 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) ) ){ - fvv <- format.round[vv] - fff <- format.full[vv] - matr[,vv] <- .write.format2( vec1=dat[,vv], ff=fff, fr=fvv ) - } - matr <- apply( matr, 1, FUN=function(ll){ paste( ll, - collapse="" ) } ) - if ( is.vector(matr) ){ - writeLines( matr, paste( savename, ".dat", sep="") ) - } else { - utils::write.table( matr, paste( savename, ".dat", sep=""), - row.names=F, col.names=F) - } - dfr <- data.frame( "variable"=colnames(dat), - "begin"=c( 1, cumsum( format.full )[ - ncol(dat) ] - + 1 ), - "end"=cumsum( format.full ), - "length"=format.full - ) - utils::write.table( dfr, paste( savename, "__LEGEND.txt",sep=""), +write.fwf2 <- function( dat, format.full, format.round, savename ) +{ + if (is.null( colnames(dat) ) ){ + colnames(dat) <- paste( "V", 1:( 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) ) ){ + fvv <- format.round[vv] + fff <- format.full[vv] + matr[,vv] <- write.format2( vec1=dat[,vv], ff=fff, fr=fvv ) + } + matr <- apply( matr, 1, FUN=function(ll){ paste( ll, collapse="" ) } ) + 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) + } + dfr <- data.frame( variable=colnames(dat), + begin=c( 1, cumsum( format.full )[ - ncol(dat) ] + 1 ), + end=cumsum( format.full), length=format.full ) + utils::write.table( dfr, paste( savename, "__LEGEND.txt",sep=""), quote=FALSE, row.names=FALSE, col.names=TRUE) + return(dfr) +} - return(dfr) - } -############################################################################## diff --git a/README.md b/README.md index 4f86477f..ca191d58 100644 --- a/README.md +++ b/README.md @@ -18,9 +18,9 @@ The CRAN version can be installed from within R using: utils::install.packages("sirt") ``` -#### GitHub version `sirt` 3.3-1 (2019-02-08) +#### GitHub version `sirt` 3.3-6 (2019-02-27) -[![](https://img.shields.io/badge/github%20version-3.3--1-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   +[![](https://img.shields.io/badge/github%20version-3.3--6-orange.svg)](https://github.com/alexanderrobitzsch/sirt)   The version hosted [here](https://github.com/alexanderrobitzsch/sirt) is the development version of `sirt`. The GitHub version can be installed using `devtools` as: diff --git a/inst/NEWS b/inst/NEWS index b9d445a2..8f6a4942 100644 --- a/inst/NEWS +++ b/inst/NEWS @@ -36,13 +36,17 @@ CHANGELOG sirt ------------------------------------------------------------------------ -VERSIONS sirt 3.3 | 2019-02-08 | Last: sirt 3.3-1 +VERSIONS sirt 3.3 | 2019-02-27 | Last: sirt 3.3-6 ------------------------------------------------------------------------ -xxx * ... +ADDED * included sampling weights for invariance alignment in function + invariance_alignment_cfa_config() +FIXED * fixed numerical instabilities in wle.rasch() + (thanks to Tom Tinkler) + DATA * included/modified datasets: --- -EXAMP * included/modified examples: --- +EXAMP * included/modified examples: invariance.alignment (7) diff --git a/man/invariance.alignment.Rd b/man/invariance.alignment.Rd index 6aef737b..85855a88 100644 --- a/man/invariance.alignment.Rd +++ b/man/invariance.alignment.Rd @@ -1,5 +1,5 @@ %% File Name: invariance.alignment.Rd -%% File Version: 1.603 +%% File Version: 1.606 \name{invariance.alignment} \alias{invariance.alignment} @@ -32,7 +32,8 @@ item intercepts and item slopes (see Example 6). The function \code{invariance_alignment_cfa_config} estimates one-factor models separately for each group as a preliminary step for invariance -alignment (see Example 6). +alignment (see Example 6). Sampling weights are accommodated by the +argument \code{weights}. } %% \emph{Note that this function is in an early experimental stage.} @@ -51,7 +52,7 @@ invariance_alignment_constraints(model, lambda_parm_tol, nu_parm_tol ) invariance_alignment_simulate(nu, lambda, err_var, mu, sigma, N) -invariance_alignment_cfa_config(dat, group, ...) +invariance_alignment_cfa_config(dat, group, weights=NULL, ...) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -640,6 +641,15 @@ mod1 <- sirt::invariance.alignment( lambda=pars$lambda, nu=pars$nu, eps=.01, cmod1 <- sirt::invariance_alignment_constraints(model=mod1, nu_parm_tol=.4, lambda_parm_tol=.2 ) summary(cmod1) + +#--- estimate CFA models with sampling weights + +#* simulate weights +weights <- stats::runif(sum(N), 0, 2) +#* estimate models +pars2 <- sirt::invariance_alignment_cfa_config(dat[,-1], group=dat$group, weights=weights) +print(pars2$nu) +print(pars$nu) } } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 30ef54c9..880d8429 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,5 +1,5 @@ //// File Name: RcppExports.cpp -//// File Version: 3.003001 +//// File Version: 3.003006 // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393