3.5-36

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.5-27
-Date: 2019-05-12 13:58:32
+Version: 3.5-36
+Date: 2019-05-15 16:30:21
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

R/gom.em.R

@@ -1,13 +1,13 @@
 ## File Name: gom.em.R
-## File Version: 5.296
+## File Version: 5.363
 
 
 #-- gom EM algorithm
 gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
     theta0.k=seq(-5,5,len=15), xsi0.k=exp(seq(-6,3,len=15)),
-    max.increment=.3, numdiff.parm=.001, maxdevchange=1e-6,
+    max.increment=.3, numdiff.parm=1e-4, maxdevchange=1e-6,
     globconv=1e-4, maxiter=1000, msteps=4, mstepconv=.001,
-    theta_adjust=TRUE, progress=TRUE )
+    theta_adjust=TRUE, lambda.inits=NULL, pi.k.inits=NULL, progress=TRUE )
 {
     CALL <- match.call()
     s1 <- Sys.time()
@@ -52,8 +52,12 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
     if (model=="GOM"){
         theta.k <- gom_em_calc_theta(K=K, problevels=problevels)
         TP <- nrow(theta.k)
-        pi.k <- as.vector(rep(1/TP, TP ))
-        lambda <- gom_em_inits_lambda(I=I, K=K)
+        if (is.null(pi.k.inits)){
+            pi.k <- as.vector(rep(1/TP, TP ))
+        } else {
+            pi.k <- pi.k.inits
+        }
+        lambda <- gom_em_inits_lambda(I=I, K=K, lambda.inits=lambda.inits)
         theta0.k <- NULL
     }
 
@@ -69,10 +73,15 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
         TP <- nrow(theta_grid)
         mu <- rep(0,K-1)
         sigma <- rep(1,K-1)
-        pi.k <- rep(1/TP, TP)
+        if (is.null(pi.k.inits)){
+            pi.k <- as.vector(rep(1/TP, TP ))
+        } else {
+            pi.k <- pi.k.inits
+        }
         theta.k <- gom_em_normal_to_membership_scores(theta_grid=theta_grid, K=K, TP=TP)
-        lambda <- gom_em_inits_lambda(I=I, K=K)
+        lambda <- gom_em_inits_lambda(I=I, K=K, lambda.inits=lambda.inits)
         theta0.k <- NULL
+        fac1 <- 1
     }
 
     # inits
@@ -84,7 +93,7 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
     disp <- "...........................................................\n"
 
     #---------- start EM algorithm --------------------
-    while( ( ( maxdevchange < devchange ) | (globconv < conv) ) & ( iter < maxiter ) ){
+    while( ( ( maxdevchange < devchange ) & (globconv < conv) ) & ( iter < maxiter ) ){
         if (progress){
             cat(disp)
             cat("Iteration", iter+1, "   ", paste( Sys.time() ), "\n" )
@@ -94,12 +103,14 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
         pi.k0 <- pi.k
         lambda0 <- lambda
         b0 <- b
+# a0 <- Sys.time()
         # calculate probabilities
         res <- gom_em_calc_probs( lambda=lambda, theta.k=theta.k, b=b,
                         theta0.k=theta0.k )
         probs <- res$probs
         probs <- problong2probarray( probres=probs, I=I, TP=TP )
 
+
         # calculate counts
         probsM <- matrix( aperm( probs, c(2,1,3) ), nrow=I*2, ncol=TP )
         res1 <- calcpost( dat2=dat2, dat2resp=dat2.resp, probs=probsM,
@@ -116,7 +127,8 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
             mu <- res$center
             Sigma <- res$cov
             #- readjust theta_grid
-            if (theta_adjust){
+            if (theta_adjust ){
+                fac1_old <- fac1
                 mth <- max(theta0)
                 sig_max <- max( sqrt(diag(Sigma)))
                 fac1 <- sig_max * ( 6 / mth    )
@@ -138,7 +150,7 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
         }
         if (model=="GOMRasch"){
             res <- gom_em_est_b( lambda=lambda, I=I, K=K, n.ik=n.ik, b=b, theta0.k=theta0.k,
-                        numdiff.parm=.001, max.increment=max.increment, theta.k=theta.k,
+                        numdiff.parm=numdiff.parm, max.increment=max.increment, theta.k=theta.k,
                         msteps=msteps, mstepconv=mstepconv, eps=.001, progress=progress )
             b <- res$b
             se.b <- res$se.b
@@ -169,6 +181,46 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
                     paste( round(max(abs(b-b0)),6), collapse=" " ), "\n", sep=""))
             }
         }
+    } #--- end EM algorithm
+
+    #--- Newton-Raphson steps: include it later
+
+    if (FALSE){
+        lambda_logit <- as.vector( stats::qlogis(lambda) )
+        pi_k_logit <- sirt_logit_probs(y=pi.k)
+        ind_lambda <- seq(1,I*K)
+        ind_pi <- max(ind_lambda) + seq(1,TP-1)
+        x0 <- c(lambda_logit, pi_k_logit)
+        #- define optimization function
+        gom_em_loglike <- function(x, ...){
+            res <- gom_em_loglike_parameter_conversion(x=x, ind_lambda=ind_lambda,
+                        ind_pi=ind_pi, I=I, K=K)
+            lambda <- res$lambda
+            pi.k <- res$pi.k
+            res <- gom_em_calc_probs( lambda=lambda, theta.k=theta.k, b=NULL,
+                                        theta0.k=theta0.k )
+            probs <- res$probs
+            probs <- problong2probarray( probres=probs, I=I, TP=TP )
+            # calculate counts
+            probsM <- matrix( aperm( probs, c(2,1,3) ), nrow=I*2, ncol=TP )
+            res1 <- calcpost( dat2=dat2, dat2resp=dat2.resp, probs=probsM,
+                                        dat2ind=dat2.ind, pik=pi.k, K=1 )
+            f.yi.qk <- res1$fyiqk
+            pik_m <- sirt_matrix2(pi.k, nrow=nrow(dat2))
+            ll <- -2*sum( log( rowSums( f.yi.qk * pik_m ) ) )
+            return(ll)
+        }
+
+        #-- optimization
+        NP <- length(x0)
+        bound <- 8
+        lower <- rep(-bound,NP)
+        upper <- rep(bound, NP)
+        res <- stats::nlminb(start=x0, objective=gom_em_loglike, control=list(trace=1),
+                    lower=lower, upper=upper)
+        par <- res$par
+        lambda <- gom_em_extract_lambda_matrix(lambda_logit=par[ind_lambda], I=I, K=K)
+        pi.k <- sirt_logit_to_probs(y=par[ind_pi])
     }
 
     #--------------- arrange OUTPUT
@@ -282,3 +334,6 @@ gom.em <- function( dat, K=NULL, problevels=NULL, model="GOM",
     class(res) <- "gom"
     return(res)
 }
+
+
+# a1 <- Sys.time(); cat("probs \n"); print(a1-a0); a0 <- a1

R/gom_em_calc_probs.R

@@ -1,5 +1,5 @@
 ## File Name: gom_em_calc_probs.R
-## File Version: 0.12
+## File Version: 0.13
 
 
 
@@ -14,7 +14,7 @@ gom_em_calc_probs <- function( lambda, theta.k, b=NULL, theta0.k=NULL )
     probsL <- array( 0, dim=c( nrow(lambda), 2, nrow(theta.k) ) )
     probsL[,2,] <- probs
     probsL[,1,] <- 1-probs
-    res <- list("probs"=probs, "probsL"=probsL)
+    res <- list(probs=probs, probsL=probsL)
     return(res)
 }
 

R/gom_em_est_lambda.R

@@ -1,10 +1,10 @@
 ## File Name: gom_em_est_lambda.R
-## File Version: 0.14
+## File Version: 0.196
 
 
 # estimation of lambda parameters
 gom_em_est_lambda <- function( lambda, I, K, n.ik, numdiff.parm=.001,
-        max.increment=1,theta.k, msteps, mstepconv, eps=.001, progress=TRUE)
+        max.increment=1, theta.k, msteps, mstepconv, eps=2*numdiff.parm, progress=TRUE)
 {
     h <- numdiff.parm
     diffindex <- 1:I
@@ -16,18 +16,32 @@ gom_em_est_lambda <- function( lambda, I, K, n.ik, numdiff.parm=.001,
     }
     it <- 0
     conv1 <- 1000
+
     while( ( it < msteps ) & ( conv1 > mstepconv ) ){
         lambda0 <- lambda
+        # nu <- stats::qlogis(lambda)
+        # nu0 <- nu
         for (kk in 1:K){
             Q1 <- Q0
             Q1[,kk] <- 1
-            pjk <- gom_em_calc_probs( lambda=lambda, theta.k=theta.k )$probsL
-            pjk1 <- gom_em_calc_probs( lambda=lambda+h*Q1, theta.k=theta.k )$probsL
-            pjk2 <- gom_em_calc_probs( lambda=lambda-h*Q1, theta.k=theta.k )$probsL
+            # parm_temp <- lambda0
+            parm_temp <- lambda
+            # lambda00 <- stats::plogis(nu0)
+            # lambda01 <- stats::plogis(nu0+h*Q1)
+            # lambda02 <- stats::plogis(nu0-h*Q1)
+            lambda00 <- parm_temp
+            lambda01 <- parm_temp + h*Q1
+            lambda02 <- parm_temp - h*Q1
+
+            pjk <- gom_em_calc_probs( lambda=lambda00, theta.k=theta.k )$probsL
+            pjk1 <- gom_em_calc_probs( lambda=lambda01, theta.k=theta.k )$probsL
+            pjk2 <- gom_em_calc_probs( lambda=lambda02, theta.k=theta.k )$probsL
             res <- gom_em_numdiff_index( pjk, pjk1, pjk2, an.ik, diffindex,
                         max.increment=max.increment, numdiff.parm )
             increment <- Q1*matrix( res$increment, nrow=I, ncol=K)
             lambda <- lambda + increment
+            # nu <- nu + increment
+            # lambda <- stats::plogis(nu)
             lambda[ lambda[,kk] < eps, kk ] <- eps
             lambda[ lambda[,kk] > 1-eps, kk ] <- 1 - eps
             se.lambda[,kk] <- sqrt(-1/res$d2)

R/gom_em_extract_lambda_matrix.R

@@ -0,0 +1,8 @@
+## File Name: gom_em_extract_lambda_matrix.R
+## File Version: 0.01
+
+gom_em_extract_lambda_matrix <- function(lambda_logit, I, K)
+{
+    lambda <- matrix( stats::plogis(lambda_logit), nrow=I, ncol=K)
+    return(lambda)
+}

R/gome_em_ic.R → R/gom_em_ic.R

@@ -1,5 +1,5 @@
-## File Name: gome_em_ic.R
-## File Version: 0.02
+## File Name: gom_em_ic.R
+## File Version: 0.03
 
 gom_em_ic <- function(dev, dat2, I, K, TP, model)
 {

R/gom_em_inits_lambda.R

@@ -1,8 +1,12 @@
 ## File Name: gom_em_inits_lambda.R
-## File Version: 0.01
+## File Version: 0.02
 
-gom_em_inits_lambda <- function(I, K)
+gom_em_inits_lambda <- function(I, K, lambda.inits=NULL)
 {
-    lambda <- matrix( .75*seq( 1/(2*K), 1, 1/K), I, K, byrow=TRUE )
+    if (is.null(lambda.inits)){
+        lambda <- matrix( .75*seq( 1/(2*K), 1, 1/K), I, K, byrow=TRUE )
+    } else {
+        lambda <- lambda.inits
+    }
     return(lambda)
 }

R/gom_em_loglike_calc_probs.R

@@ -0,0 +1,15 @@
+## File Name: gom_em_loglike_calc_probs.R
+## File Version: 0.03
+
+gom_em_loglike_calc_probs <- function(x, ind_pi, ind_lambda, I, K, theta.k,
+    theta0.k)
+{
+    pi_k_logit <- x[ind_pi]
+    lambda_logit <- x[ind_lambda]
+    pi.k <- sirt_logit_to_probs(y=pi_k_logit)
+    lambda <- gom_em_extract_lambda_matrix(lambda_logit=lambda_logit, I=I, K=K)
+    res <- gom_em_calc_probs( lambda=lambda, theta.k=theta.k, b=NULL,
+                                theta0.k=theta0.k )
+    probs <- res$probs
+    return(probs)
+}

R/gom_em_loglike_parameter_conversion.R

@@ -0,0 +1,14 @@
+## File Name: gom_em_loglike_parameter_conversion.R
+## File Version: 0.02
+
+
+gom_em_loglike_parameter_conversion <- function(x, ind_lambda, ind_pi, I, K)
+{
+    pi_k_logit <- x[ind_pi]
+    lambda_logit <- x[ind_lambda]
+    pi.k <- sirt_logit_to_probs(y=pi_k_logit)
+    lambda <- gom_em_extract_lambda_matrix(lambda_logit=lambda_logit, I=I, K=K)
+    #--- output
+    res <- list(lambda=lambda, pi.k=pi.k)
+    return(res)
+}

R/gom_em_numdiff_index.R

@@ -1,11 +1,11 @@
 ## File Name: gom_em_numdiff_index.R
-## File Version: 0.13
+## File Version: 0.14
 
 
 #--- general function for numerical differentiation
 #--- diffindex aggregates across super items
 gom_em_numdiff_index <- function( pjk, pjk1, pjk2, an.ik, diffindex,
-        max.increment, numdiff.parm, eps=1e-20 )
+        max.increment, numdiff.parm, eps=1e-5 )
 {
     h <- numdiff.parm
     ll0 <- rowSums( an.ik * log(pjk+eps) )

R/lsem.bootstrap.R

@@ -1,5 +1,5 @@
 ## File Name: lsem.bootstrap.R
-## File Version: 0.295
+## File Version: 0.300
 
 
 lsem.bootstrap <- function(object, R=100, verbose=TRUE)
@@ -48,7 +48,7 @@ lsem.bootstrap <- function(object, R=100, verbose=TRUE)
 
     #- include new objects in output
     object$parameters_boot <- parameters_boot
-    object$parameters <- object$parameters
+    object$parameters <- parameters
     object$R <- R
     object$class_boot <- TRUE
     object$fitstats_joint <- fitstats_joint

R/rasch_mml2_difference_quotient.R

@@ -1,7 +1,7 @@
 ## File Name: rasch_mml2_difference_quotient.R
-## File Version: 0.11
+## File Version: 0.12
 
-rasch_mml2_difference_quotient <- function(ll0, ll1, ll2, h, eps=1e-10)
+rasch_mml2_difference_quotient <- function(ll0, ll1, ll2, h, eps=1e-6)
 {
     # first order derivative
     # f(x+h) - f(x-h)=2*f'(x)*h

R/sirt_add_increment.R

@@ -0,0 +1,9 @@
+## File Name: sirt_add_increment.R
+## File Version: 0.01
+
+sirt_add_increment <- function(x, pos, value)
+{
+    y <- x
+    y[pos] <- x[pos] + value
+    return(y)
+}

R/sirt_logit_probs.R

@@ -0,0 +1,9 @@
+## File Name: sirt_logit_probs.R
+## File Version: 0.01
+
+sirt_logit_probs <- function(y)
+{
+    K <- length(y)
+    y_logit <- log( y[-K] / y[K] )
+    return(y_logit)
+}

R/sirt_logit_to_probs.R

@@ -0,0 +1,13 @@
+## File Name: sirt_logit_to_probs.R
+## File Version: 0.01
+
+
+sirt_logit_to_probs <- function(y)
+{
+    K1 <- length(y)
+    x <- rep(0,K1+1)
+    x[1:K1] <- y
+    x <- exp(x)
+    x <- x / sum(x)
+    return(x)
+}

README.md

@@ -22,9 +22,9 @@ The CRAN version can be installed from within R using:
 utils::install.packages("sirt")
 ```
 
-#### GitHub version `sirt` 3.5-27 (2019-05-12)
+#### GitHub version `sirt` 3.5-36 (2019-05-15)
 
-[![](https://img.shields.io/badge/github%20version-3.5--27-orange.svg)](https://github.com/alexanderrobitzsch/sirt)  
+[![](https://img.shields.io/badge/github%20version-3.5--36-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: