3.10-79

alexanderrobitzsch · Mar 9, 2021 · 213afad · 213afad
1 parent 6bc98a1
commit 213afad
Show file tree

Hide file tree

Showing 38 changed files with 486 additions and 399 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.10-70
-Date: 2020-09-05 19:52:03
+Version: 3.10-79
+Date: 2021-03-09 13:01:29
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

diff --git a/NAMESPACE b/NAMESPACE
@@ -106,6 +106,7 @@ importFrom(CDM, IRT.posterior)
 importFrom(CDM, IRT.predict)
 importFrom(CDM, IRT.repDesign)
 importFrom(CDM, IRT.se)
+importFrom(CDM, IRT_RMSD_calc_rmsd)
 importFrom(CDM, modelfit.cor2)
 importFrom(CDM, numerical_gradient)
 importFrom(CDM, numerical_Hessian)
@@ -406,6 +407,7 @@ S3method(coef, rasch.evm.pcm)
 S3method(coef, xxirt)
 S3method(confint, xxirt)
 S3method(IRT.expectedCounts, MultipleGroupClass)
+S3method(IRT.expectedCounts, rasch.mml)
 S3method(IRT.expectedCounts, SingleGroupClass)
 S3method(IRT.expectedCounts, xxirt)
 S3method(IRT.factor.scores, rm.facets)

diff --git a/R/IRT.expectedCounts_sirt.R b/R/IRT.expectedCounts_sirt.R
@@ -1,9 +1,9 @@
 ## File Name: IRT.expectedCounts_sirt.R
-## File Version: 0.07
+## File Version: 0.16
 
 
-###########################################################
-# object of class xxirt
+
+#*** object of class xxirt
 IRT.expectedCounts.xxirt <- function( object, ... )
 {
     ll <- object$n.ik
@@ -12,4 +12,20 @@ IRT.expectedCounts.xxirt <- function( object, ... )
     attr(ll,"G") <- object$G
     return(ll)
 }
-###########################################################
+
+#*** object of class rasch.mml
+IRT.expectedCounts.rasch.mml <- function( object, ... )
+{
+    njk <- object$n.jk
+    rjk <- object$r.jk
+    dims <- dim(njk)
+    ll <- array(0, dim=c(dims[1],2,dims[2],dims[3]))
+    ll[,2,,] <- rjk
+    ll[,1,,] <- njk-rjk
+    attr(ll,"theta") <- object$theta.k
+    attr(ll,"prob.theta") <- object$pi.k
+    attr(ll,"G") <- object$G
+    attr(ll,"dimnames") <- list()
+    attr(ll,"dimnames")[[1]] <- colnames(object$dat)
+    return(ll)
+}
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.010070
+## File Version: 3.010079
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 

diff --git a/R/data.prep.R b/R/data.prep.R
@@ -1,67 +1,62 @@
 ## File Name: data.prep.R
-## File Version: 1.13
+## File Version: 1.144
 
-#--------------------------------------------------------------------------
-# -----------------------------------------------------------
-# data preparations for rasch.jml and rasch.mml
-.data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE ){
-        #-------------------
-        # freq.patt ... should frequency pattern be taken into account?
-        #               default=TRUE
-        #-------------------
-        # should items being excluded?
-# a0 <- Sys.time()
-        item.means <- colMeans( dat, na.rm=T )
-        item.elim <- which( item.means %in% c(0,1))
-        if ( length( item.elim ) > 0 ){
-                stop( cat( paste( "There are", length(item.elim), "Items with no variance!") ) )
-                    }
-        if ( any( is.na( item.means )) ){ stop( "There are items which contains only missings!") }
-         n <- nrow(dat)
-        I <- ncol(dat)
-        if( is.null(weights) ){  weights <- rep( 1, n ) }
-        # indicator for nonmissing response
-        dat.9 <- dat
-        dat.9[ is.na(dat) ] <- 9
-         # pattern
-        if ( use.freqpatt ){         #
-            freq.patt <- apply( dat.9, 1, FUN=function(ll){ paste(ll, collapse="" ) } )  #
-#        freq.patt <- dat.9[,1]
-#        for (ii in 2:I){ freq.patt <- paste0( freq.patt, dat.9[,ii] ) }
-        # frequency pattern with frequencies
-                    dat1 <- data.frame( table( freq.patt ) )
-# cat("a220") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
-                        } else {
-                    freq.patt <- paste("FP", 1000000 + 1:n, sep="")
-                    dat1 <- data.frame( freq.patt )
-                    colnames(dat1)[1] <- "freq.patt"
-                            }
-         # weighting the frequencies if survey weights are supplied
-        if( !is.null(weights) ){
-            # standardize weights
-            weights <- weights / sum( weights ) * n
-            if ( use.freqpatt ){
-                    dat1[,2] <- stats::aggregate( weights, list( freq.patt), sum )[,2]
-                        } else {
-                    dat1[,"Freq"] <- weights
-                        }
-                    }
-         # item pattern corresponding to frequency pattern
-        if ( use.freqpatt){
-            dat2 <- matrix( as.numeric( unlist( strsplit( paste(dat1[,1]), "" ) ) ), ncol=ncol(dat), byrow=T)
-                    } else {
-                    dat2 <- dat.9 }
+#----- data preparations for rasch.jml and rasch.mml
+data.prep <- function( dat, weights=NULL, use.freqpatt=TRUE,
+        standardize_weights=TRUE)
+{
+    item.means <- colMeans( dat, na.rm=TRUE )
+    item.elim <- which( item.means %in% c(0,1))
+    if ( length( item.elim ) > 0 ){
+        stop( cat( paste( "There are", length(item.elim), "Items with no variance!") ) )
+    }
+    if ( any( is.na(item.means)) ){
+        stop( "There are items which contains only missings!")
+    }
+    n <- nrow(dat)
+    I <- ncol(dat)
+    if( is.null(weights) ){  weights <- rep( 1, n ) }
+    # indicator for nonmissing response
+    dat.9 <- dat
+    dat.9[ is.na(dat) ] <- 9
+
+    #* pattern
+    if ( use.freqpatt ){
+        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="")
+        dat1 <- data.frame( freq.patt )
+        colnames(dat1)[1] <- "freq.patt"
+    }
+    # weighting the frequencies if survey weights are supplied
+    if( !is.null(weights) ){
+        # standardize weights
+        if (standardize_weights){
+            weights <- weights / sum(weights) * n
+        }
+        if ( use.freqpatt ){
+            dat1[,2] <- stats::aggregate( weights, list( freq.patt), sum )[,2]
+        } else {
+            dat1[,"Freq"] <- weights
+        }
+    }
+    # item pattern corresponding to frequency pattern
+    if ( use.freqpatt){
+        dat2 <- matrix( as.numeric( unlist( strsplit( paste(dat1[,1]), "" ) ) ), ncol=ncol(dat), byrow=T)
+    } else {
+        dat2 <- dat.9 }
         dat2.resp <- 1 * ( dat2 !=9 )
         dat2[ dat2==9 ] <- 0
-         # mean right
+        # 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" )
         list( "dat"=dat, "dat2"=dat2, "dat2.resp"=dat2.resp, "dat1"=dat1,
-                "freq.patt"=freq.patt, "I"=I, "n"=n,
-                "dat9"=dat.9  )
-        }
-    #*******************
+                "freq.patt"=freq.patt, "I"=I, "n"=n, "dat9"=dat.9 )
+}
+
+#*******************
     #   OUTPUT:
     #   dat         ... original data
     #   dat2        ... reduced original data. Each different item resonse is represented by one row.
@@ -71,33 +66,35 @@
     #   I           ... number of items
     #   n           ... number of subjects
     #   dat9        ... This is the original data exact from the fact that missings are recoded by 9.
-#-----------------------------------------------------------------
-# -------------------------------------------------------------------
 
-#--------------------------------------------------------
-# Small function which helps for printing purposes
-.prnum <- function( matr, digits ){
-        VV <- ncol(matr)
-        for (vv in 1:VV){
-        # vv <- 1
+.data.prep <- data.prep
+
+
+
+#*** Small function which helps for printing purposes
+.prnum <- function( matr, digits )
+{
+    VV <- ncol(matr)
+    for (vv in 1:VV){
         if ( is.numeric( matr[,vv]) ){ matr[,vv] <- round( matr[,vv], digits ) }
-                        }
-        print(matr)
-                   }
-#--------------------------------------------------------
+    }
+    print(matr)
+}
+
 
 
 
-############################################################
-# Function for calculation of a response pattern
+
+#-- Function for calculation of a response pattern
 #   for dichotomous responses
-resp.pattern2 <- function( x ){
+resp.pattern2 <- function(x)
+{
     n <- nrow(x)
     p <- ncol(x)
     mdp <- (x %*% (2^((1:ncol(x)) - 1))) + 1
     misspattern <- mdp[,1]
-    misspattern <- list( "miss.pattern"=mdp[,1],
-                "mp.index"=match( mdp[,1], sort( unique(mdp[,1] ) ) ) )
-    return( misspattern )
-        }
-########################################################
+    misspattern <- list( miss.pattern=mdp[,1],
+                        mp.index=match( mdp[,1], sort( unique(mdp[,1] ) ) ) )
+    return(misspattern)
+}
+
diff --git a/R/lq_fit.R b/R/lq_fit.R
@@ -1,5 +1,5 @@
 ## File Name: lq_fit.R
-## File Version: 0.147
+## File Version: 0.152
 
 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),
@@ -11,7 +11,7 @@ lq_fit <- function(y, X, w=NULL, pow=2, eps=1e-3, beta_init=NULL,
     }
     ind <- ! is.na(y)
     y <- y[ind]
-    X <- X[ind,]
+    X <- X[ind,,drop=FALSE]
     w <- w[ind]
     if (is.null(beta_init)){
         mod <- stats::lm.wfit(y=y, x=X, w=w)

diff --git a/R/mle.rasch.R b/R/mle.rasch.R
@@ -0,0 +1,46 @@
+## File Name: mle.rasch.R
+## File Version: 1.03
+
+
+
+#---------------------------------------------------------
+# Maximum Likelihood Estimation (Rasch model)
+
+mle.rasch <- function( dat, dat.resp=1-is.na(dat), b, theta, conv=.001,
+            progress=FALSE, prior_sd=NULL)
+{
+    theta.change <- 1
+    if ( progress){ cat("\n  MLE estimation  |" ) }
+    I <- length(b)
+    bM <- matrix( b, nrow=length(theta), length(b), byrow=TRUE )
+
+    # prior for ability
+    if (! is.null(prior_sd) ){
+        is_prior <- TRUE
+    } else {
+        is_prior <- FALSE
+    }
+
+    while( max( abs( theta.change) > conv )){
+        # calculate P and Q
+        p.ia <- stats::plogis( theta - bM )
+        q.ia <- 1 - p.ia
+        # Likelihood
+        l1 <- rowSums( dat.resp* ( dat - p.ia ) )
+        # derivative of the objective function
+        f1.obj <- rowSums( - dat.resp * p.ia * q.ia  )
+        # add prior
+        if (is_prior){
+            l1 <- l1 - theta / prior_sd^2
+            f1.obj <- f1.obj - 1 / prior_sd^2
+        }
+
+        # theta change
+        theta.change <- - l1 / f1.obj
+        theta <- theta + theta.change
+        if ( progress){  cat("-") }
+    }
+    res <- list( "theta"=theta, "p.ia"=p.ia )
+    return(res)
+}
+
diff --git a/R/mml_raschtype_counts.R b/R/mml_raschtype_counts.R
@@ -1,5 +1,5 @@
 ## File Name: mml_raschtype_counts.R
-## File Version: 0.01
+## File Version: 0.08