3.6-8

alexanderrobitzsch · Jun 27, 2020 · 9496f03 · 9496f03
1 parent 281d926
commit 9496f03
Show file tree

Hide file tree

Showing 93 changed files with 2,323 additions and 2,085 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,10 +1,12 @@
 Package: TAM
 Type: Package
 Title: Test Analysis Modules
-Version: 3.6-2
-Date: 2020-05-06 12:34:35
+Version: 3.6-8
+Date: 2020-06-27 21:47:07
 Author:
-    Alexander Robitzsch [aut, cre], Thomas Kiefer [aut], Margaret Wu [aut]
+    Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>),
+    Thomas Kiefer [aut], 
+    Margaret Wu [aut]
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 
     Includes marginal maximum likelihood estimation and joint maximum

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

diff --git a/R/tam.fa.R b/R/tam.fa.R
@@ -1,10 +1,10 @@
 ## File Name: tam.fa.R
-## File Version: 9.256
+## File Version: 9.258
 
 
 #---- Exploratory Factor Analysis and Bifactor Models
 tam.fa <- function( resp, irtmodel, dims=NULL, nfactors=NULL,
-            pid=NULL,pweights=NULL, verbose=TRUE, control=list() )
+            pid=NULL,pweights=NULL, verbose=TRUE, control=list(), ... )
 {
     require_namespace_msg("GPArotation")
     require_namespace_msg("psych")
@@ -91,11 +91,11 @@ tam.fa <- function( resp, irtmodel, dims=NULL, nfactors=NULL,
     if ( irtmodel %in% c("bifactor2","efa") ){
         res <- tam.mml.2pl( resp=resp, Q=Q, irtmodel=irtmodel2,
                     variance.fixed=variance.fixed, pid=pid,
-                    pweights=pweights, control=con )
+                    pweights=pweights, control=con, ... )
     }
     if ( irtmodel=="bifactor1"){
         res <- tam.mml( resp=resp, Q=Q, variance.fixed=variance.fixed, pid=pid,
-                    pweights=pweights, control=con )
+                    pweights=pweights, control=con, ... )
     }
     #****
     # calculate standardized loadings

diff --git a/R/tam.linking.R b/R/tam.linking.R
@@ -1,5 +1,5 @@
 ## File Name: tam.linking.R
-## File Version: 0.345
+## File Version: 0.347
 
 tam.linking <- function( tamobj_list, type="Hae", method="joint",
     pow_rob_hae=1, eps_rob_hae=1e-4, theta=NULL, wgt=NULL, wgt_sd=2, fix.slope=FALSE,

diff --git a/R/tam.mml.3pl.R b/R/tam.mml.3pl.R
@@ -1,5 +1,5 @@
 ## File Name: tam.mml.3pl.R
-## File Version: 9.857
+## File Version: 9.872
 
 tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
             formulaY=NULL, dataY=NULL,
@@ -51,9 +51,10 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
     resp <- as.matrix(resp)
     resp0 <- resp <- add.colnames.resp(resp)
 
-    #********************
-    # create E design matrix from different input matrices
-    res0 <- tam_mml_3pl_create_E( resp, E, Q, gammaslope.des, Q.fixed )
+    #--- create E design matrix from different input matrices
+    E_null <- is.null(E)
+    res0 <- tam_mml_3pl_create_E( resp=resp, E=E, Q=Q,
+                    gammaslope.des=gammaslope.des, Q.fixed=Q.fixed )
     E <- res0$E
     if ( is.null(gammaslope.fixed ) ){
         gammaslope.fixed <- res0$gammaslope.fixed
@@ -86,6 +87,7 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
     }
     B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E )
 
+
     #***********************
     if ( is.null(A)){ printxsi <- FALSE  } else { printxsi <- TRUE }
 
@@ -170,8 +172,8 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
     if ( snodes > 0 ){ nnodes <- snodes }
 
     #--- print information about nodes
-    res <- tam_mml_progress_proc_nodes( progress=progress, snodes=snodes, nnodes=nnodes,
-                    skillspace=skillspace, QMC=QMC )
+    res <- tam_mml_progress_proc_nodes( progress=progress, snodes=snodes,
+                    nnodes=nnodes, skillspace=skillspace, QMC=QMC )
 
     # maximum no. of categories per item. Assuming dichotomous
     maxK <- max( resp, na.rm=TRUE ) + 1
@@ -355,6 +357,7 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
       gammaslope <- .mml.3pl.gammaslope.center( gammaslope, gammaslope.center.index,
                          gammaslope.center.value  )
 
+
     #******
     # prior distribution guessing parameter
     if ( ! is.null(guess.prior) ){
@@ -363,8 +366,8 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
         guess[ i1 ] <- guess.mean[i1]
         guess.prior[ guess.prior==0 ] <- 1E-3
     }
-    #******
-    # prior distribution slope parameter
+
+    #---- prior distribution slope parameter
     if ( ( ! is.null(gammaslope.prior) ) & ( ! init.gammaslope) ){
         i1 <- which( gammaslope.prior[,2] < 10 )
         gammaslope[ i1 ] <- gammaslope.prior[i1,1]
@@ -396,9 +399,9 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
     }
     AXsi[iIndex,] <- AXsi.tmp[,,1]
 
-    #****
-    # compute B matrix
-    B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E )
+    #--- compute B matrix
+    B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E, B=B,
+                    skip_B=FALSE)
 
     ##**SE
     se.xsi <- 0*xsi
@@ -560,7 +563,8 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
         if ( skillspace !="normal" ){
             res <- tam_mml_3pl_skillspace( Ngroup=Ngroup, pi.k=pi.k,
                         delta.designmatrix=delta.designmatrix, G=G, delta=delta,
-                        delta.fixed=delta.fixed, hwt=hwt, resp.ind=resp.ind, pweightsM=pweightsM,
+                        delta.fixed=delta.fixed, hwt=hwt, resp.ind=resp.ind,
+                        pweightsM=pweightsM,
                         pweights=pweights, group1.list=group1.list,
                         delta_acceleration=delta_acceleration, iter=iter )
             pi.k <- res$pi.k
@@ -575,8 +579,8 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
 
         #******
         # generate input for fixed parameters
-    xsi.fixed.estimated <- tam_generate_xsi_fixed_estimated( xsi=xsi, A=A )
-    B.fixed.estimated <- tam_generate_B_fixed_estimated(B=B)
+        xsi.fixed.estimated <- tam_generate_xsi_fixed_estimated( xsi=xsi, A=A )
+        B.fixed.estimated <- tam_generate_B_fixed_estimated(B=B)
 
         ######################################
         # calculation of expected counts
@@ -604,22 +608,26 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
 
 # cat("\nM steps intercepts") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
 
-      ###############################################
-      # M-step item slopes
+
+        ###############################################
+        # M-step item slopes
         if ( est.some.slopes){
             oldgamma <- gammaslope
-            res <- tam_mml_3pl_mstep_item_slopes( max.increment=max.increment, np=np,
-                        Msteps=Msteps, nitems=nitems, A=A, AXsi=AXsi, B=B, xsi=xsi, guess=guess,
-                        theta=theta, nnodes=nnodes, maxK=maxK, progress=progress, ItemScore=ItemScore,
-                        fac.oldxsi=fac.oldxsi, rprobs=rprobs, xsi.fixed=xsi.fixed, convM=convM,
-                        rprobs0=rprobs0, n.ik=n.ik, N.ik=N.ik, gammaslope=gammaslope, E=E,
-                        FdesM=FdesM, dimFdes=dimFdes, gammaslope.fixed=gammaslope.fixed,
-                        gammaslope.prior=gammaslope.prior, maxgamma=maxgamma, Edes=Edes,
-                        gammaslope.constr.V=gammaslope.constr.V, V1=V1, e2=e2,
-                        gammaslope.center.index=gammaslope.center.index,
+            res <- tam_mml_3pl_mstep_item_slopes( max.increment=max.increment,
+                        np=np, Msteps=Msteps, nitems=nitems, A=A, AXsi=AXsi, B=B,
+                        xsi=xsi, guess=guess, theta=theta, nnodes=nnodes,
+                        maxK=maxK, progress=progress, ItemScore=ItemScore,
+                        fac.oldxsi=fac.oldxsi, rprobs=rprobs, xsi.fixed=xsi.fixed,
+                        convM=convM, rprobs0=rprobs0, n.ik=n.ik, N.ik=N.ik,
+                        gammaslope=gammaslope, E=E, FdesM=FdesM, dimFdes=dimFdes,
+                        gammaslope.fixed=gammaslope.fixed,
+                        gammaslope.prior=gammaslope.prior, maxgamma=maxgamma,
+                        Edes=Edes, gammaslope.constr.V=gammaslope.constr.V, V1=V1,
+                        e2=e2, gammaslope.center.index=gammaslope.center.index,
                         gammaslope.center.value=gammaslope.center.value,
                         userfct.gammaslope=userfct.gammaslope,
-                        gammaslope_acceleration=gammaslope_acceleration, V=V )
+                        gammaslope_acceleration=gammaslope_acceleration, V=V,
+                        skip_B=FALSE)
             gammaslope <- res$gammaslope
             se.gammaslope <- res$se.gammaslope
             gammaslope.change <- res$gammachange
@@ -628,8 +636,8 @@ tam.mml.3pl <- function( resp, Y=NULL, group=NULL,
             B <- res$B
         }
 
-# cat("\nM steps slopes") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
 
+# cat("\nM steps slopes") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
 
         #--- guessing parameter estimation
         if ( est.some.guess ){

diff --git a/R/tam.np.R b/R/tam.np.R
@@ -1,5 +1,5 @@
 ## File Name: tam.np.R
-## File Version: 0.413
+## File Version: 0.421
 
 
 tam.np <- function( dat, probs_init=NULL, pweights=NULL, lambda=NULL, control=list(),
@@ -103,9 +103,9 @@ tam.np <- function( dat, probs_init=NULL, pweights=NULL, lambda=NULL, control=li
                     lambda=lambda, iter=iter, dev=dev, n_basis=n_basis,
                     penalty_type=penalty_type, target_fct=target_fct,
                     index_target=index_target, index_basis=index_basis, model=model)
+            probs <- res$probs
         }
         spline_optim <- res$spline_optim
-        probs <- res$probs
         pars <- res$pars
         n_reg <- res$n_reg
         par_reg_penalty <- res$par_reg_penalty
@@ -123,11 +123,12 @@ tam.np <- function( dat, probs_init=NULL, pweights=NULL, lambda=NULL, control=li
         #-- print progress
         res <- tam_mml_progress_em(progress=progress, deviance=dev,
                     deviance_change=deviance_change, iter=iter,
-                    rel_deviance_change=rel_deviance_change, is_mml_3pl=FALSE, xsi_change=0,
-                    beta_change=0, variance_change=0, B_change=0, skillspace='np', delta_change=0,
-                    digits_pars=6, devch=devch,    penalty_xsi=0, is_np=TRUE, is_latreg=TRUE,
-                    np_change=np_change, par_reg_penalty=par_reg_penalty, n_reg=n_reg,
-                    AIC=AIC, n_est=n_est, n_reg_max=n_reg_max)
+                    rel_deviance_change=rel_deviance_change, is_mml_3pl=FALSE,
+                    xsi_change=0, beta_change=0, variance_change=0, B_change=0,
+                    skillspace='np', delta_change=0, digits_pars=6, devch=devch,
+                    penalty_xsi=0, is_np=TRUE, is_latreg=TRUE, np_change=np_change,
+                    par_reg_penalty=par_reg_penalty, n_reg=n_reg, AIC=AIC, 
+                    n_est=n_est, n_reg_max=n_reg_max)
         #-- convergence
         iter <- iter + 1
         if (iter > maxiter){
@@ -168,13 +169,14 @@ tam.np <- function( dat, probs_init=NULL, pweights=NULL, lambda=NULL, control=li
     #--- output
     s2 <- Sys.time()
     time <- c(s1, s2)
-    res <- list( CALL=CALL, dat=dat, dat2=dat2, dat_resp=dat_resp, n.ik=n.ik, N.ik=N.ik,
-                item=item, rprobs=probs, pi.k=pi.k, nodes=nodes, pweights=pweights, like=f.yi.qk,
-                hwt=f.qk.yi, iter=iter, loglike=loglike, AIC=AIC, converged=converged,
-                iter=iter, time=time, dev=dev, theta=theta, G=1, pars=pars, n_est=n_est,
-                n_reg=n_reg, regularized=regularized, basis_type=basis_type,
-                n_basis=n_basis, desmat=desmat, ic=ic, pid=NULL,
-                orthonormalize=orthonormalize, penalty_type=penalty_type,
+    res <- list( CALL=CALL, dat=dat, dat2=dat2, dat_resp=dat_resp, n.ik=n.ik, 
+                N.ik=N.ik, item=item, rprobs=probs, pi.k=pi.k, nodes=nodes, 
+                pweights=pweights, like=f.yi.qk, hwt=f.qk.yi, iter=iter, 
+                loglike=loglike, AIC=AIC, converged=converged,
+                iter=iter, time=time, dev=dev, theta=theta, G=1, pars=pars, 
+                n_est=n_est, n_reg=n_reg, regularized=regularized, 
+                basis_type=basis_type, n_basis=n_basis, desmat=desmat, ic=ic, 
+                pid=NULL, orthonormalize=orthonormalize, penalty_type=penalty_type,
                 pen_val=pen_val, use_basis=use_basis, model=model, sigma=sigma,
                 ll_individual=ll_individual, control=control)
     class(res) <- "tam.np"

diff --git a/R/tam_linking_irf_discrepancy.R b/R/tam_linking_irf_discrepancy.R
@@ -1,10 +1,10 @@
 ## File Name: tam_linking_irf_discrepancy.R
-## File Version: 0.054
+## File Version: 0.055
 
 tam_linking_irf_discrepancy <- function(probs1, probs2, wgt, type,
     pow_rob_hae=1, eps_rob_hae=1e-4)
 {
-    K <- dim(probs1)[3]
+    K <- min( dim(probs1)[3], dim(probs2)[3])
     crit <- 0
     #-- define Haebara criterion function
     if (type %in% c("Hae","RobHae") ){

diff --git a/R/tam_linking_joint.R b/R/tam_linking_joint.R
@@ -1,5 +1,5 @@
 ## File Name: tam_linking_joint.R
-## File Version: 0.078
+## File Version: 0.082
 
 tam_linking_joint <- function(NM, parameters_list, linking_args, verbose=TRUE)
 {
@@ -44,7 +44,6 @@ tam_linking_joint <- function(NM, parameters_list, linking_args, verbose=TRUE)
     if (!is.null(par_init)){
         par <- par_init
     }
-
     linking_criterion_multiple_studies <- function(x){
         bvec <- c(0, x[1:(NM-1)])
         avec <- c(1, x[NM-1 + 1:(NM-1)])

diff --git a/R/tam_linking_joint_calc_probs.R b/R/tam_linking_joint_calc_probs.R
@@ -1,5 +1,5 @@
 ## File Name: tam_linking_joint_calc_probs.R
-## File Version: 0.06
+## File Version: 0.09
 
 tam_linking_joint_calc_probs <- function(a, b, parameters_list_mm, theta)
 {

diff --git a/R/tam_mml_3pl_computeB.R b/R/tam_mml_3pl_computeB.R
@@ -1,12 +1,14 @@
 ## File Name: tam_mml_3pl_computeB.R
-## File Version: 0.07
+## File Version: 0.11
 
 
 # function for computation of item loadings
-tam_mml_3pl_computeB <- function( Edes, gammaslope, E )
+tam_mml_3pl_computeB <- function( Edes, gammaslope, E, skip_B=FALSE, B=NULL )
 {
-    B <- tam_rcpp_mml_3pl_compute_B( Edes=Edes, gammaslope=gammaslope,
+    if (! skip_B){
+        B <- tam_rcpp_mml_3pl_compute_B( Edes=Edes, gammaslope=gammaslope,
                     dimE=dim(E) )$B
-    B <- array( B, dim(E)[1:3] )
+        B <- array( B, dim(E)[1:3] )
+    }
     return(B)
 }
diff --git a/R/tam_mml_3pl_create_E.R b/R/tam_mml_3pl_create_E.R
@@ -1,11 +1,9 @@
 ## File Name: tam_mml_3pl_create_E.R
-## File Version: 0.06
+## File Version: 0.07
 
 
-####################################################
-# create E matrix
-tam_mml_3pl_create_E <- function( resp, E, Q, gammaslope.des,
-        Q.fixed=NULL )
+#--- create E matrix
+tam_mml_3pl_create_E <- function( resp, E, Q, gammaslope.des, Q.fixed=NULL )
 {
     Qdes <- NULL
     gammaslope.fixed <- NULL

diff --git a/R/tam_mml_3pl_mstep_item_slopes.R b/R/tam_mml_3pl_mstep_item_slopes.R
@@ -1,15 +1,16 @@
 ## File Name: tam_mml_3pl_mstep_item_slopes.R
-## File Version: 9.63
+## File Version: 9.645
 
-########################################
-# tam.mml.3pl estimate item slopes
+
+#--- tam.mml.3pl estimate item slopes
 tam_mml_3pl_mstep_item_slopes <- function( max.increment, np,
             Msteps, nitems, A, AXsi, B, xsi, guess, theta, nnodes, maxK,
             progress,ItemScore, fac.oldxsi, rprobs, xsi.fixed, convM, rprobs0,
             n.ik, N.ik, gammaslope, E, FdesM, dimFdes,
             gammaslope.fixed, gammaslope.prior, maxgamma=9.99, Edes,
             gammaslope.constr.V, V1, e2, gammaslope.center.index,
-            gammaslope.center.value, userfct.gammaslope, gammaslope_acceleration, V )
+            gammaslope.center.value, userfct.gammaslope, gammaslope_acceleration, V,
+            skip_B=FALSE)
 {
     if (progress){
         cat("\nM Step Slopes       |")
@@ -26,7 +27,8 @@ tam_mml_3pl_mstep_item_slopes <- function( max.increment, np,
 
     while( ( iter <=msteps ) & ( parchange > convM)  ){
         Xlambda0 <- gammaslope <- Xlambda
-        B <- tam_mml_3pl_computeB( Edes, gammaslope, E )
+        B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E,
+                        B=B, skip_B=skip_B)
 
         # calculate probabilities
         res <- tam_mml_3pl_calc_prob( iIndex=1:nitems, A=A, AXsi=AXsi, B=B, xsi=xsi,
@@ -128,8 +130,9 @@ tam_mml_3pl_mstep_item_slopes <- function( max.increment, np,
     gammaslope <- fac.oldxsi * gammaslope0 + ( 1 - fac.oldxsi)*gammaslope
 
     #--- gammaslope acceleration
-    if ( gammaslope_acceleration$acceleration !="none" ){
-        gammaslope_acceleration <- tam_accelerate_parameters( xsi_acceleration=gammaslope_acceleration,
+    if ( gammaslope_acceleration$acceleration!="none" ){
+        gammaslope_acceleration <- tam_accelerate_parameters(
+                        xsi_acceleration=gammaslope_acceleration,
                         xsi=gammaslope, iter=iter, itermin=3)
         gammaslope <- gammaslope_acceleration$parm
     }
@@ -138,15 +141,16 @@ tam_mml_3pl_mstep_item_slopes <- function( max.increment, np,
     gammaslope_change <- tam_parameter_change( gammaslope, gammaslope0)
 
     #--- recompute B
-    B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E )
+    B <- tam_mml_3pl_computeB( Edes=Edes, gammaslope=gammaslope, E=E,
+                    skip_B=skip_B, B=B )
 
     #---- OUTPUT
-    res <- list("gammaslope"=Xlambda, "se.gammaslope"=se.Xlambda,
-                 "max.increment.b"=max.increment,
-                 "gammachange"=max( abs( Xlambda00 - Xlambda) ),
-                 gammaslope_change=gammaslope_change,
-                 gammaslope_acceleration=gammaslope_acceleration, B=B
-                 )
+    res <- list(gammaslope=Xlambda, se.gammaslope=se.Xlambda,
+                max.increment.b=max.increment,
+                gammachange=max( abs( Xlambda00 - Xlambda) ),
+                gammaslope_change=gammaslope_change,
+                gammaslope_acceleration=gammaslope_acceleration, B=B
+                )
     return(res)
 }
 #----------------------------------------------------------