3.4-20

alexanderrobitzsch · Mar 4, 2020 · 88fcfb2 · 88fcfb2
1 parent fd08c13
commit 88fcfb2
Show file tree

Hide file tree

Showing 23 changed files with 198 additions and 97 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: TAM
 Type: Package
 Title: Test Analysis Modules
-Version: 3.4-17
-Date: 2020-02-14 19:15:05
+Version: 3.4-20
+Date: 2020-03-04 13:59:55
 Author:
     Alexander Robitzsch [aut, cre], Thomas Kiefer [aut], Margaret Wu [aut]
 Maintainer: Alexander Robitzsch <[email protected]>

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.004017
+## File Version: 3.004020
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 

diff --git a/R/summary.tam.jml.R b/R/summary.tam.jml.R
@@ -1,14 +1,15 @@
 ## File Name: summary.tam.jml.R
-## File Version: 9.253
+## File Version: 9.2573
 
 
 #***** summary for tam object
 summary.tam.jml <- function( object, file=NULL, ...)
 {
     #* open sink
     tam_osink(file=file)
+    sdisplay <- tam_summary_display()
 
-    cat("------------------------------------------------------------\n")
+    cat(sdisplay)
     cat( tam_packageinfo("TAM"), "\n" )
     cat( tam_rsessinfo(), "\n\n")
 
@@ -22,23 +23,38 @@ summary.tam.jml <- function( object, file=NULL, ...)
     # print Call
     tam_print_call(object$CALL)
 
-    cat("------------------------------------------------------------\n")
+    cat(sdisplay)
     cat( "Number of iterations", "=", object$iter, "\n\n" )
     cat( "Deviance", "=", round( object$deviance, 2 ), " | " )
     cat( "Log Likelihood", "=", round( -object$deviance/2, 2 ), "\n" )
     cat( "Number of persons", "=", object$nstud, "\n" )
-
     if( ! is.null( object$formulaA)  ){
         cat( "Number of generalized items", "=", object$nitems, "\n" )
         cat( "Number of items", "=", ncol(object$resp_orig), "\n" )
     } else {
         cat( "Number of items", "=", object$nitems, "\n" )
     }
+    cat( "constraint", "=", object$constraint, "\n" )
+    cat( "bias", "=", object$bias, "\n" )
+
+    obji <- object$theta_summary
+    if (obji$ndim==1){
+        cat(sdisplay)
+        cat("Person Parameters xsi\n")
+        cat( "M", "=", round( obji$M, 2 ), "\n" )
+        cat( "SD", "=", round( obji$SD, 2 ), "\n" )
+    }
 
-    cat("\nItem Parameters xsi\n")
+    cat(sdisplay)
+    cat("Item Parameters xsi\n")
     obji <- object$item
     tam_round_data_frame_print(obji=obji, digits=3, from=2)
 
+    cat(sdisplay)
+    cat("Item Parameters -A*Xsi\n")
+    obji <- object$item1
+    tam_round_data_frame_print(obji=obji, from=2, to=ncol(obji), digits=3, rownames_null=TRUE)
+
     #** close sink
     tam_csink(file=file)
 }
diff --git a/R/tam.jml.R b/R/tam.jml.R
@@ -1,11 +1,12 @@
 ## File Name: tam.jml.R
-## File Version: 9.355
+## File Version: 9.3575
 
 
 tam.jml <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
             bias=TRUE, xsi.fixed=NULL,  xsi.inits=NULL,  theta.fixed=NULL,
             A=NULL, B=NULL, Q=NULL, ndim=1,
-            pweights=NULL, verbose=TRUE, control=list(), version=2 )
+            pweights=NULL, constraint="cases",
+            verbose=TRUE, control=list(), version=2 )
 {
     CALL <- match.call()
     #**** handle verbose argument
@@ -16,20 +17,41 @@ tam.jml <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     if ( ! is.null(theta.fixed) ){
         version <- 1
     }
+    if (! is.null(B)){
+        if (dim(B)[3]>1){
+            version <- 1
+        }
+    }
+
     #**** version=1
     if (version==1){
+        constraint <- "cases"
         res <- tam_jml_version1( resp=resp, group=group, adj=adj,
                     disattenuate=disattenuate, bias=bias, xsi.fixed=xsi.fixed,
                     xsi.inits=xsi.inits, A=A, B=B, Q=Q, ndim=ndim, theta.fixed=theta.fixed,
-                    pweights=pweights, control=control )
+                    pweights=pweights, control=control )                    
     }
     #**** version=2
     if (version==2){
         res <- tam_jml_version2( resp=resp, group=group, adj=adj,
-                    disattenuate=disattenuate, bias=bias, xsi.fixed=xsi.fixed,  xsi.inits=xsi.inits,
-                    A=A, B=B, Q=Q, ndim=ndim, pweights=pweights, control=control  )
+                    disattenuate=disattenuate, bias=bias, xsi.fixed=xsi.fixed,
+                    xsi.inits=xsi.inits, A=A, B=B, Q=Q, ndim=ndim,
+                    pweights=pweights, control=control, constraint=constraint )
     }
+
+    #- process item parameters
+    res$AXsi <- tam_jml_compute_Axsi(A=res$A, xsi=res$xsi, resp=resp)
+    #- item parameter table
+    res$item1 <- tam_jml_itempartable( resp=resp, maxK=res$maxK, AXsi=res$AXsi,
+                    B=res$B, resp.ind=res$resp.ind)
+    #- theta summary
+    res$theta_summary <- tam_jml_proc_abilities(theta=res$theta,
+                                pweights=res$pweights, B=res$B)
+
+    #- output
     res$CALL <- CALL
     res$resp <- resp
+    res$constraint <- constraint
+    res$bias <- bias
     return(res)
 }
diff --git a/R/tam.mml.R b/R/tam.mml.R
@@ -1,5 +1,5 @@
 ## File Name: tam.mml.R
-## File Version: 9.797
+## File Version: 9.8002
 
 tam.mml <- function( resp, Y=NULL, group=NULL,  irtmodel="1PL",
             formulaY=NULL, dataY=NULL,
@@ -76,7 +76,9 @@ tam.mml <- function( resp, Y=NULL, group=NULL,  irtmodel="1PL",
       A <- .A.PCM2( resp, constraint=constraint, Q=Q  )
     }
 
-    if ( !is.null(con$seed)){ set.seed( con$seed )     }
+    if ( !is.null(con$seed)){
+        set.seed( con$seed )
+    }
 
     nitems <- ncol(resp)       # number of items
     if (is.null(colnames(resp))){

diff --git a/R/tam_itempartable.R b/R/tam_itempartable.R
@@ -1,29 +1,28 @@
 ## File Name: tam_itempartable.R
-## File Version: 9.12
+## File Version: 9.131
 
 
 
-##################################################
-# create table of item parameters
+#--- create table of item parameters
 tam_itempartable <- function( resp, maxK, AXsi, B, ndim,
-            resp.ind, rprobs, n.ik, pi.k, order=FALSE ){
+            resp.ind, rprobs=NULL, n.ik=NULL, pi.k=NULL, order=FALSE )
+{
 
     if ( is.null(dimnames(B)[[1]] ) ){
         dimnames(B)[[1]] <- colnames(resp)
     }
 
     item1 <- data.frame( "item"=dimnames(B)[[1]] )
-    item1$N <- colSums(resp.ind )
-    item1$M <- colSums( resp.ind * resp, na.rm=TRUE) / colSums( resp.ind )
-    maxKi <- rowSums( 1 - is.na( AXsi ) ) - 1
+    item1$N <- colSums(resp.ind)
+    item1$M <- colSums(resp.ind * resp, na.rm=TRUE) / colSums( resp.ind )
+    maxKi <- rowSums( 1 - is.na(AXsi) ) - 1
     I <- nrow(item1)
     item1$xsi.item <- - AXsi[ cbind(1:I, maxKi+1) ] / maxKi
 
-    #****
-    # Item fit
+    #--- Item fit
     # probs ... [ classes, items, categories ]
-    probs <- aperm( rprobs, perm=c(3,1,2))
-    pi.k <- matrix( pi.k, ncol=1 )
+    # probs <- aperm( rprobs, perm=c(3,1,2))
+    # pi.k <- matrix( pi.k, ncol=1 )
     b0 <- sum( B[, 1, ], na.rm=TRUE )
     a0 <- 0
     if ( b0 + a0 > 0 ){

diff --git a/R/tam_jml_compute_Axsi.R b/R/tam_jml_compute_Axsi.R
@@ -0,0 +1,16 @@
+## File Name: tam_jml_compute_Axsi.R
+## File Version: 0.03
+
+tam_jml_compute_Axsi <- function(A, xsi, resp)
+{
+    dim_A <- dim(A)
+    K <- dim_A[2]
+    I <- dim_A[1]
+    AXsi <- matrix(NA, nrow=I, ncol=K)
+    colnames(AXsi) <- paste0("Cat",0:(K-1) )
+    rownames(AXsi) <- colnames(resp)
+    for (kk in 1:K){
+        AXsi[,kk] <- A[,kk,] %*% xsi
+    }
+    return(AXsi)
+}
diff --git a/R/tam_jml_itempartable.R b/R/tam_jml_itempartable.R
@@ -0,0 +1,12 @@
+## File Name: tam_jml_itempartable.R
+## File Version: 0.02
+
+
+tam_jml_itempartable <- function(resp, maxK, AXsi, B, resp.ind)
+{
+    ndim <- dim(B)[3]
+    item <- tam_itempartable( resp=resp, maxK=maxK, AXsi=AXsi,
+                    B=B, ndim=ndim, resp.ind=resp.ind, rprobs=NULL,
+                    n.ik=NULL, pi.k=NULL, order=FALSE )
+    return(item)
+}
diff --git a/R/tam_jml_proc_abilities.R b/R/tam_jml_proc_abilities.R
@@ -0,0 +1,16 @@
+## File Name: tam_jml_proc_abilities.R
+## File Version: 0.02
+
+
+tam_jml_proc_abilities <- function(theta, pweights, B)
+{
+    ndim <- dim(B)[3]
+    if (is.vector(theta)){
+        M <- weighted_mean(x=theta, w=pweights)
+        SD <- weighted_sd(x=theta, w=pweights)
+    } else {
+        M <- NULL
+        SD <- NULL
+    }
+    res <- list(M=M, SD=SD, ndim=ndim)
+}
diff --git a/R/tam_jml_version1.R b/R/tam_jml_version1.R
@@ -1,22 +1,13 @@
 ## File Name: tam_jml_version1.R
-## File Version: 9.34
+## File Version: 9.356
 
 tam_jml_version1 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
                      bias=TRUE, xsi.fixed=NULL,  xsi.inits=NULL,
                      theta.fixed=NULL,
                      A=NULL, B=NULL, Q=NULL, ndim=1,
-                     pweights=NULL, control=list()
-                     # control can be specified by the user
-){
-
-  #------------------------------------
-  # INPUT:
-  # control:
-  #      control=list( nodes=seq(-6,6,len=15),
-  #                              convD=.001,conv=.0001, convM=.0001, Msteps=30,
-  #                   maxiter=1000, progress=TRUE)
-  # progress ... if TRUE, then display progress
-  #-------------------------------------
+                     pweights=NULL, control=list() )
+{
+
 
   maxiter <- conv <- progress <- tamobj <- convM <- Msteps <- NULL
   R <- NULL
@@ -107,14 +98,6 @@ tam_jml_version1 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     indexIP.list[[kk]] <- which( indexIP[,kk] > 0 )
   }
 
-  # These sufficient statistics must be changed
-  # to make it more general
-  # First extension:  pweights and dependent on A; needs to be further extended (e.g., different number of categories)
-  # Second extension: multiple category option       -> resp \in 0:maxKi (see method definition calc_posterior_TK)
-  #                                                  -> length(ItemScore)=np (see diff computation in M Step)
-  #                   multiple category option Bugfix
-  #                                                  -> dim(cResp)=(nstud, nitems*maxK)
-  #                                                  -> adapt dim(A) to dim(cResp) for sufficient statistic (cf. print.designMatrices)
 
   col.index <- rep( 1:nitems, each=maxK )
   cResp <- resp[, col.index  ]*resp.ind[, col.index ]
@@ -127,7 +110,6 @@ tam_jml_version1 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
   cB[is.na(cB)] <- 0
 
   # Item sufficient statistics
-  # ItemScore <- (cResp %*% cA) %t*% pweights
   ItemScore <- crossprod(cResp %*% cA, pweights )
 
   # Computer possible maximum parameter score for each person
@@ -193,15 +175,13 @@ tam_jml_version1 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     }
     olddeviance <- deviance
 
-    #**********************
-    #update theta, ability estimates
-
+    #-- update theta (ability estimates)
     jmlAbility <- tam_jml_wle( tamobj, resp, resp.ind, A, B, nstud, nitems, maxK, convM,
-                 PersonScores, theta, xsi, Msteps, WLE=FALSE,
-                 theta.fixed=theta.fixed)
+                    PersonScores, theta, xsi, Msteps, WLE=FALSE,
+                    theta.fixed=theta.fixed)
 
     theta <- jmlAbility$theta
-    if (is.null( xsi.fixed)){
+    if (is.null(xsi.fixed)){
         theta <- theta - mean(theta)
     }
     meanChangeWLE <- jmlAbility$meanChangeWLE
@@ -299,7 +279,7 @@ tam_jml_version1 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
                "PersonMax"=PersonMaxB, "ItemMax"=ItemMax,
                "deviance"=deviance, "deviance.history"=deviance.history,
                "resp"=resp, "resp.ind"=resp.ind, "group"=group,
-               "pweights"=pweights, "A"=A, "B"=B,
+               "pweights"=pweights, "A"=A, "B"=B, AXsi=AXsi, 
                "nitems"=nitems, "maxK"=maxK,
                "nstud"=nstud, "resp.ind.list"=resp.ind.list,
                "xsi.fixed"=xsi.fixed, "deviance"=deviance,

diff --git a/R/tam_jml_version2.R b/R/tam_jml_version2.R
@@ -1,9 +1,10 @@
 ## File Name: tam_jml_version2.R
-## File Version: 9.503
+## File Version: 9.514
 
 tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
             bias=TRUE, xsi.fixed=NULL,  xsi.inits=NULL, A=NULL, B=NULL, Q=NULL,
-            ndim=1, pweights=NULL, control=list())
+            ndim=1, pweights=NULL, control=list(), constraint="cases",
+            damp=.1)
 {
 
     maxiter <- conv <- progress <- tamobj <- convM <- Msteps <- NULL
@@ -34,7 +35,7 @@ tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
 
     #-- create design matrices
     design <- designMatrices( modeltype="PCM", maxKi=NULL, resp=resp,
-                        A=A, B=B, Q=Q, R=R, ndim=ndim )
+                        A=A, B=B, Q=Q, R=R, ndim=ndim, constraint=constraint )
     A <- design$A
     A.0 <- A
     A.0[ is.na(A.0) ] <- 0
@@ -43,7 +44,6 @@ tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     B.0[ is.na(B.0) ] <- 0
     cA <- design$flatA
     cA[is.na(cA)] <- 0
-
     # number of parameters
     np <- dim(A)[[3]]
     errorP <- rep(0,np)
@@ -151,6 +151,9 @@ tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     #Initialise theta (WLE) values for all students
     theta <- log(PersonScores/(PersonMaxB-PersonScores)) #log of odds ratio of raw score
 
+    # center theta?
+    center_theta <- is.null(xsi.fixed) & (constraint=="cases")
+
     deviance <- 0
     deviance.history <- matrix( 0, nrow=maxiter, ncol=2)
     colnames(deviance.history) <- c("iter", "deviance")
@@ -164,6 +167,7 @@ tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
     #-------------- start iterations -------------------------------
     while ( (( maxthetachange > conv) | (maxChangeP > conv))  & (iter < maxiter) ) {
 
+        xsi_old <- xsi
         iter <- iter + 1
         if (progress){
             cat(disp)
@@ -178,12 +182,13 @@ tam_jml_version2 <- function( resp, group=NULL, adj=.3, disattenuate=FALSE,
                             A=A, B=B, nstud=nrow(rp3.sel), nitems=nitems, maxK=maxK,
                             convM=convM, PersonScores=PersonScores[ rp3.sel$caseid ],
                             theta=theta[ rp3.sel$caseid,, drop=FALSE],
-                            xsi=xsi, Msteps=Msteps, WLE=FALSE)
+                            xsi=xsi, Msteps=Msteps, WLE=FALSE, damp=damp)
         theta <- jmlAbility$theta
         theta <- theta[ rp3$theta.index,, drop=FALSE]
-
-        if (is.null( xsi.fixed)){
+        if (center_theta){
             theta <- theta - mean(theta)
+        } else {
+            damp <- 1 - (1-damp)*.99
         }
         meanChangeWLE <- jmlAbility$meanChangeWLE
         maxthetachange <- max( abs( theta - theta_old ) )