4.2-93

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 4.2-89
-Date: 2024-11-13 10:10:05
+Version: 4.2-93
+Date: 2024-11-28 08:56:47
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 

NAMESPACE

@@ -251,6 +251,7 @@ export(linking.haberman)
 export(linking.haberman.lq)
 export(linking.haebara)
 export(linking.robust)
+export(locpolycor)
 export(lq_fit)
 export(lq_fit_estimate_power)
 export(lsdm)

R/RcppExports.R

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

R/brm.sim.R

@@ -1,5 +1,5 @@
 ## File Name: brm.sim.R
-## File Version: 0.123
+## File Version: 0.124
 
 
 #-- brm.sim
@@ -8,12 +8,11 @@ brm.sim <- function( theta, delta, tau, K=NULL)
     I <- length(delta)
     N <- length(theta)
     dat <- matrix( 0, nrow=N, ncol=I )
-    colnames(dat) <- paste0( 'I', 1:9 )
+    colnames(dat) <- paste0( 'I', 1L:9 )
     if ( ! is.null(K) ){
         br <- seq( 0, 1, len=K+1 )
     }
     for (ii in 1L:I){
-        # ii <- 1
         m1 <- exp( ( theta - delta[ii] + tau[ii] ) / 2 )
         n1 <- exp( ( - theta + delta[ii] + tau[ii] ) / 2 )
         dat[,ii] <- stats::rbeta( N, shape1=m1, shape2=n1 )

R/locpolycor.R

@@ -0,0 +1,124 @@
+## File Name: locpolycor.R
+## File Version: 0.251
+
+
+locpolycor <- function(y, data.mod, moderator.grid, h=1.1,
+                    model_thresh, model_polycor, kernel="gaussian",
+                    eps=1e-10)
+{
+    #- compute weights
+    llw <- lsem_local_weights(data.mod=x, moderator.grid=moderator.grid, h=h,
+                kernel=kernel)
+    weights_grid <- llw$weights
+
+    #- estimate thresholds
+    y <- as.matrix(y)
+    I <- ncol(y)
+    G <- length(moderator.grid)
+
+    #- item-wise optimization
+    thresh_list <- as.list(1L:I)
+    K_list <- as.list(1L:I)
+    thresh_ind_list <- as.list(1L:I)
+    items <- colnames(y)
+    names(thresh_list) <- items
+    names(thresh_ind_list) <- items
+    names(K_list) <- items
+    Kmax <- max(y, na.rm=TRUE)
+    thresh_stat <- list()
+    thresh0 <- matrix(NA, nrow=Kmax, ncol=I)
+    rownames(thresh0) <- paste0('t',1L:Kmax)
+    colnames(thresh0) <- items
+    for (gg in 1L:G){
+        thresh_stat[[gg]] <- thresh0
+    }
+
+    for (ii in 1L:I){
+        cat(paste0( '-- compute thresholds for item ', ii, '\n') )
+        utils::flush.console()
+        y1 <- y[,ii]
+        ind <- which( ! is.na(y1) )
+        K <- max(y1)
+        thresh_item <- matrix(NA, nrow=K, ncol=G)
+        thresh_ind_cbind <- NULL
+        data_mod <- data.mod[ind]
+        res_item <- list()
+        par_init <- NULL
+        for (gg in 1L:G){
+            x0 <- moderator.grid[gg]
+            w <- weights_grid[ind, gg]
+            res <- locpolycor_est_thresh_item(y=y1, data.mod=data_mod, x0=x0,
+                        w=w, model=model_thresh, par_init=par_init, eps=eps)
+            par_init <- res$res_optim$par
+            thresh_ind_cbind <- cbind( thresh_ind_cbind, res$thresh_ind)
+            thresh_item[,gg] <- res$thresh
+            res_item[[gg]] <- res
+            thresh_stat[[gg]][1L:K,ii] <- res$thresh
+        }
+
+        thresh_list[[ii]] <- thresh_item
+
+        K_list[[ii]] <- K
+
+        #* individual predictions of thresholds
+        thresh_ind <- matrix(NA, nrow=N, ncol=K)
+        distmat <- abs( outer( data_mod, moderator.grid, '-' ) )
+        smallest <- rowKSmallest2.sirt(matr=distmat, K=2)
+        smallest <- smallest$smallind[,1L:2]
+
+        # define predictions
+        thresh_ind <- res_item[[1]]$thresh_ind
+        ind <- which(data_mod > moderator.grid[G])
+        if (length(ind)>0){
+            thresh_ind[ind, ] <- ( res_item[[G]]$thresh_ind )[ ind,, drop=FALSE ]
+        }
+
+        ind <- which( ( data_mod <=moderator.grid[G]) & (data_mod >=moderator.grid[1]))
+        sm1 <- smallest[,1]
+        sm2 <- smallest[,2]
+        mg1 <- moderator.grid[sm1]
+        mg2 <- moderator.grid[sm2]
+        t3 <- t2 <- t1 <- matrix(NA, nrow=N, ncol=K)
+        for (kk in 1L:K){
+            mat1 <- cbind( 1L:N, kk+K*(sm1-1) )
+            t1[,kk] <- thresh_ind_cbind[ mat1 ]
+            mat2 <- cbind( 1L:N, kk+K*(sm2-1) )
+            t2[,kk] <- thresh_ind_cbind[ mat2 ]
+            # linear interpolation of individual thresholds
+            t3[,kk] <- t1[,kk] + ( t2[,kk]-t1[,kk]) * (data_mod-mg1) / (mg2-mg1)
+        }
+        thresh_ind[ind,] <- t3[ind,]
+        thresh_ind_list[[ii]] <- thresh_ind
+
+    }  # end ii
+
+    #-- local polychoric correlations
+    polycor0 <- diag(1,I)
+    rownames(polycor0) <- colnames(polycor0) <- items
+    polycor_stat <- list()
+    for (gg in 1L:G){
+        polycor_stat[[gg]] <- polycor0
+    }
+    for (ii in 1L:(I-1) ){
+        for (jj in (ii+1):I){
+            cat(paste0( '** compute polychoric correlation for item pair (',
+                                ii, ',', jj, ')\n') )
+            res <- locpolycor_est_polycor_itempair(y=y, ii=ii, jj=jj, data.mod=data.mod,
+                            moderator.grid=moderator.grid, weights_grid=weights_grid,
+                            model=model_polycor, thresh_ind_list=thresh_ind_list,
+                            x0=x0, eps=eps)
+            polycor1 <- res$polycor
+            for (gg in 1L:G){
+                polycor_stat[[gg]][ii,jj] <- polycor1[[gg]]
+                polycor_stat[[gg]][jj,ii] <- polycor1[[gg]]
+            }
+        }
+    }
+
+    #-- output
+    res <- list(thresh_list=thresh_list, thresh_stat=thresh_stat,
+                    polycor_stat=polycor_stat, thresh_ind_list=thresh_ind_list,
+                    K_list=K_list, I=I, moderator.grid=moderator.grid,
+                    weights_grid=weights_grid)
+    return(res)
+}

R/locpolycor_est_polycor_itempair.R

@@ -0,0 +1,53 @@
+## File Name: locpolycor_est_polycor_itempair.R
+## File Version: 0.089
+
+
+locpolycor_est_polycor_itempair <- function(y, ii, jj, data.mod, moderator.grid,
+        weights_grid, thresh_ind_list, x0, model="const",
+        lower=-.99, upper=.99, EE=99.99, par_init=0.2,
+        eps=1e-10, optimizer="optim")
+{
+    y1 <- y[,ii]
+    y2 <- y[,jj]
+    data_mod <- data.mod
+    G <- length(moderator.grid)
+    par0 <- par_init
+    if (model=='lin'){
+        par0 <- c(par0, 0)
+    }
+    K1 <- max(y1)
+    K2 <- max(y2)
+    N <- length(y1)
+    x1 <- data_mod
+
+    W <- polycor <- rep(NA,G)
+
+    for (gg in 1L:G){
+        x0 <- moderator.grid[gg]
+        w <- weights_grid[,gg]
+        thresh_ii <- thresh_ind_list[[ii]]
+        thresh_jj <- thresh_ind_list[[jj]]
+        thresh_ii <- cbind(-EE, thresh_ii, EE)
+        thresh_jj <- cbind(-EE, thresh_jj, EE)
+        thresh_ii1 <- thresh_ii[ cbind(1L:N, y1+2) ]
+        thresh_ii2 <- thresh_ii[ cbind(1L:N, y1+1) ]
+        thresh_jj1 <- thresh_jj[ cbind(1L:N, y2+2) ]
+        thresh_jj2 <- thresh_jj[ cbind(1L:N, y2+1) ]
+
+        res_optim <- sirt_optimizer(optimizer=optimizer, par=par0,
+                            fn=locpolycor_est_polycor_opt_fun,
+                            w=w, thresh_ii1=thresh_ii1,
+                            thresh_ii2=thresh_ii2, thresh_jj1=thresh_jj1,
+                            thresh_jj2=thresh_jj2, eps=eps,
+                            x1=data.mod, x0=x0, model=model, method='L-BFGS-B',
+                            lower=lower, upper=upper, hessian=FALSE,
+                            package='pbivnorm')
+        par0 <- res_optim$par
+        polycor[gg] <- par0[1]
+        W[gg] <- sum(w)
+    }
+
+    #-- output
+    res <- list( polycor=polycor, ii=ii, jj=jj, W=W, N=N)
+    return(res)
+}

R/locpolycor_est_polycor_opt_fun.R

@@ -0,0 +1,33 @@
+## File Name: locpolycor_est_polycor_opt_fun.R
+## File Version: 0.084
+
+
+locpolycor_est_polycor_opt_fun <- function(x, w, x1, x0,
+        thresh_ii1, thresh_ii2, thresh_jj1, thresh_jj2,
+        model="const", eps=1e-10,
+        package="pbivnorm", maxcor=0.99)
+{
+    requireNamespace(package)
+    if (package=='pbivnorm'){
+        pbvnorm_fun <- pbivnorm::pbivnorm
+    }
+    if (package=='pbv'){
+        pbvnorm_fun <- pbv::pbvnorm
+    }
+    if (model=='const'){
+        rho1 <- x
+    }
+    if (model=='lin'){
+        rho1 <- x[1] + x[2] * (x1-x0)
+        rho1 <- pmin( maxcor, pmax(-maxcor, rho1) )
+    }
+
+    d11 <- pbvnorm_fun(x=thresh_ii1, y=thresh_jj1, rho=rho1)
+    d10 <- pbvnorm_fun(x=thresh_ii1, y=thresh_jj2, rho=rho1)
+    d01 <- pbvnorm_fun(x=thresh_ii2, y=thresh_jj1, rho=rho1)
+    d00 <- pbvnorm_fun(x=thresh_ii2, y=thresh_jj2, rho=rho1)
+    val <- d11 - d10 - d01 + d00
+    val <- ifelse(val<eps, eps, val)
+    val <- - sum( w*log(val+eps) )
+    return(val)
+}

R/locpolycor_est_thresh_grad_fun.R

@@ -0,0 +1,49 @@
+## File Name: locpolycor_est_thresh_grad_fun.R
+## File Version: 0.141
+
+
+locpolycor_est_thresh_grad_fun <- function(x, y, x1, w, x0, model, K, eps=1e-10)
+{
+    if (model=='const'){
+        par1 <- c(-Inf,x,Inf)
+        pred1 <- par1[y+2]
+        pred2 <- par1[y+1]
+    }
+    if (model=='lin'){
+        par1 <- c(-Inf,x[1L:K],Inf)
+        par2 <- c(-Inf,x[K+(1L:K)],Inf)
+        pred1 <- par1[y+2] + par2[y+2]*(x1-x0)
+        pred1 <- ifelse(y==K, par1[y+2], pred1)
+        pred2 <- par1[y+1] + par2[y+1]*(x1-x0)
+        pred2 <- ifelse(y==0, par2[y+1], pred2)
+    }
+    arg <- stats::pnorm(pred1) - stats::pnorm(pred2)
+    dpred1 <- stats::dnorm(pred1)
+    dpred2 <- stats::dnorm(pred2)
+    if (model=='lin'){
+        lpred1 <- stats::dnorm(pred1)*(x1-x0)
+        lpred2 <- stats::dnorm(pred2)*(x1-x0)
+    }
+
+    arg <- ifelse( arg < eps, eps, arg )
+    NP <- length(x)
+    grad <- rep(NA, NP)
+    for (pp in 1L:NP){
+        # P(Y=0)=Phi(tau1) - 0
+        # P(Y=1)=P(Y<=1) - P(Y=0)=Phi(tau2) - Phi(tau1), ...
+        arg1 <- dpred1
+        arg2 <- dpred2
+        uu <- pp
+        if (model=='lin' & (pp>(NP/2) )){
+            arg1 <- lpred1
+            arg2 <- lpred2
+            uu <- pp - NP/2
+        }
+        # d log(arg) / d par=d par / arg
+        a1 <- ifelse(y==uu-1, arg1, 0)
+        a2 <- ifelse(y==uu, -arg2, 0)
+        ad <- a1+a2
+        grad[pp] <- - sum(w*ad/arg)
+    }
+    return(grad)
+}

R/locpolycor_est_thresh_item.R

@@ -0,0 +1,55 @@
+## File Name: locpolycor_est_thresh_item.R
+## File Version: 0.189
+
+
+locpolycor_est_thresh_item <- function(y, data.mod, x0, w, model, par_init=NULL,
+        eps=1e-10, optimizer="optim", use_deriv=TRUE)
+{
+    K <- max(y)
+    N <- length(y)
+    x1 <- data.mod
+    # starting values
+    par0 <- stats::qnorm( seq(1/(K+1), K/(K+1), length=K) )
+    if (model=='lin'){
+        par1 <- 0*par0
+        names(par1) <- paste0('lin', 1L:K)
+        par0 <- c(par0, par1)
+    }
+    if (! is.null(par_init) ){
+        par0 <- par_init
+    }
+    #- optimization
+    if (use_deriv){
+        grad <- locpolycor_est_thresh_grad_fun
+    } else {
+        grad <- NULL
+    }
+
+    res_optim <- sirt_optimizer(optimizer=optimizer, par=par0,
+                        fn=locpolycor_est_thresh_opt_fun,
+                        grad=locpolycor_est_thresh_grad_fun,
+                        y=y, x1=x1, x0=x0, w=w, model=model, method='L-BFGS-B',
+                        K=K, eps=eps, hessian=FALSE )
+
+    #- local thresholds across all persons at x0
+    thresh <- res_optim$par
+    if (model=='lin'){
+        thresh <- thresh[1L:K]
+    }
+    #- individual threshold predictions
+    thresh_ind <- matrix(NA, nrow=N, ncol=K)
+    colnames(thresh_ind) <- paste0('Cat',1L:K)
+    for (kk in 1L:K){
+        if (model=='const'){
+            thresh_ind[,kk] <- thresh[kk]
+        }
+        if (model=='lin'){
+            thresh_ind[,kk] <- thresh[kk] + res_optim$par[K+kk] * (x1-x0)
+        }
+    }
+
+    #*** output
+    res <- list(thresh=thresh, thresh_ind=thresh_ind, x0=x0, res_optim=res_optim, N=N,
+                        W=sum(w), K=K, model=model)
+    return(res)
+}

R/locpolycor_est_thresh_opt_fun.R

@@ -0,0 +1,25 @@
+## File Name: locpolycor_est_thresh_opt_fun.R
+## File Version: 0.153
+
+
+
+locpolycor_est_thresh_opt_fun <- function(x, y, x1, w, x0, model, K, eps=1e-10)
+{
+    if (model=='const'){
+        par1 <- c(-Inf,x,Inf)
+        pred1 <- par1[y+2]
+        pred2 <- par1[y+1]
+    }
+    if (model=='lin'){
+        par1 <- c(-Inf,x[1L:K],Inf)
+        par2 <- c(-Inf,x[K+(1L:K)],Inf)
+        pred1 <- par1[y+2] + par2[y+2]*(x1-x0)
+        pred1 <- ifelse(y==K, par1[y+2], pred1)
+        pred2 <- par1[y+1] + par2[y+1]*(x1-x0)
+        pred2 <- ifelse(y==0, par2[y+1], pred2)
+    }
+    arg <- stats::pnorm(pred1) - stats::pnorm(pred2)
+    arg <- ifelse( arg < eps, eps, arg )
+    val <- - sum(w*log(arg))
+    return(val)
+}