rasch.copula2: Multidimensional IRT Copula Model

Description

This function handles local dependence by specifying copulas for residuals in multidimensional item response models for dichotomous item responses (Braeken, 2011; Braeken, Tuerlinckx & de Boeck, 2007; Schroeders, Robitzsch & Schipolowski, 2014). Estimation is allowed for item difficulties, item slopes and a generalized logistic link function (Stukel, 1988).

The function rasch.copula3 allows the estimation of multidimensional models while rasch.copula2 only handles unidimensional models.

Usage

rasch.copula2(dat, itemcluster, copula.type="bound.mixt",
    progress=TRUE, mmliter=1000, delta=NULL,
    theta.k=seq(-4, 4, len=21), alpha1=0, alpha2=0,
    numdiff.parm=1e-06,  est.b=seq(1, ncol(dat)),
    est.a=rep(1, ncol(dat)), est.delta=NULL, b.init=NULL, a.init=NULL,
    est.alpha=FALSE, glob.conv=0.0001, alpha.conv=1e-04, conv1=0.001,
    dev.crit=.2, increment.factor=1.01)
rasch.copula3(dat, itemcluster, dims=NULL, copula.type="bound.mixt",
    progress=TRUE, mmliter=1000, delta=NULL,
    theta.k=seq(-4, 4, len=21), alpha1=0, alpha2=0,
    numdiff.parm=1e-06,  est.b=seq(1, ncol(dat)),
    est.a=rep(1, ncol(dat)), est.delta=NULL, b.init=NULL, a.init=NULL,
    est.alpha=FALSE, glob.conv=0.0001, alpha.conv=1e-04, conv1=0.001,
    dev.crit=.2, rho.init=.5, increment.factor=1.01)
# S3 method for rasch.copula2
summary(object,...)
# S3 method for rasch.copula3
summary(object,...)
# S3 method for rasch.copula2
anova(object,...)
# S3 method for rasch.copula3
anova(object,...)
# S3 method for rasch.copula2
logLik(object,...)
# S3 method for rasch.copula3
logLik(object,...)
# S3 method for rasch.copula2
IRT.likelihood(object,...)
# S3 method for rasch.copula3
IRT.likelihood(object,...)
# S3 method for rasch.copula2
IRT.posterior(object,...)
# S3 method for rasch.copula3
IRT.posterior(object,...)

Arguments

dat

An $N \times I$ data frame. Cases with only missing responses are removed from the analysis.

itemcluster

An integer vector of length $I$ (number of items). Items with the same integers define a joint item cluster of (positively) locally dependent items. Values of zero indicate that the corresponding item is not included in any item cluster of dependent responses.

dims

A vector indicating to which dimension an item is allocated. The default is that all items load on the first dimension.

copula.type

A character or a vector containing one of the following copula types: bound.mixt (boundary mixture copula), cook.johnson (Cook-Johnson copula) or frank (Frank copula) (see Braeken, 2011). The vector copula.type must match the number of different itemclusters. For every itemcluster, a different copula type may be specified (see Examples).

progress

Print progress? Default is TRUE.

mmliter

Maximum number of iterations.

delta

An optional vector of starting values for the dependency parameter delta.

theta.k

Discretized trait distribution

alpha1

alpha1 parameter in the generalized logistic item response model (Stukel, 1988). The default is 0 which leads together with alpha2=0 to the logistic link function.

alpha2

alpha2 parameter in the generalized logistic item response model

numdiff.parm

Parameter for numerical differentiation

est.b

Integer vector of item difficulties to be estimated

est.a

Integer vector of item discriminations to be estimated

est.delta

Integer vector of length length(itemcluster). Nonzero integers correspond to delta parameters which are estimated. Equal integers indicate parameter equality constraints.

b.init

Initial $b$ parameters

a.init

Initial $a$ parameters

est.alpha

Should both alpha parameters be estimated? Default is FALSE.

glob.conv

Convergence criterion for all parameters

alpha.conv

Maximal change in alpha parameters for convergence

conv1

Maximal change in item parameters for convergence

dev.crit

Maximal change in the deviance. Default is .2.

rho.init

Initial value for off-diagonal elements in correlation matrix

increment.factor

A numeric value larger than one which controls the size of increments in iterations. To stabilize convergence, choose values 1.05 or 1.1 in some situations.

object

Object of class rasch.copula2 or rasch.copula3

…

Further arguments to be passed

Value

A list with following entries

N.itemclusters

Number of item clusters

item

Estimated item parameters

iter

Number of iterations

dev

Deviance

delta

Estimated dependency parameters $δ$

Estimated item difficulties

Estimated item slopes

Mean

sigma

Standard deviation

alpha1

Parameter $α_{1}$ in the generalized item response model

alpha2

Parameter $α_{2}$ in the generalized item response model

Information criteria

theta.k

Discretized ability distribution

pi.k

Fixed $θ$ distribution

deviance

Deviance

pattern

Item response patterns with frequencies and posterior distribution

person

Data frame with person parameters

datalist

List of generated data frames during estimation

EAP.rel

Reliability of the EAP

copula.type

Type of copula

summary.delta

Summary for estimated $δ$ parameters

f.qk.yi

Individual posterior

f.yi.qk

Individual likelihood

…

References

Braeken, J. (2011). A boundary mixture approach to violations of conditional independence. Psychometrika, 76, 57-76.

Braeken, J., & Tuerlinckx, F. (2009). Investigating latent constructs with item response models: A MATLAB IRTm toolbox. Behavior Research Methods, 41, 1127-1137.

Braeken, J., Tuerlinckx, F., & De Boeck, P. (2007). Copulas for residual dependencies. Psychometrika, 72, 393-411.

Schroeders, U., Robitzsch, A., & Schipolowski, S. (2014). A comparison of different psychometric approaches to modeling testlet structures: An example with C-tests. Journal of Educational Measurement, 51(4), 400-418.

Stukel, T. A. (1988). Generalized logistic models. Journal of the American Statistical Association, 83, 426-431.

Examples

Run this code

# NOT RUN {
#############################################################################
# EXAMPLE 1: Reading Data
#############################################################################

data(data.read)
dat <- data.read

# define item clusters
itemcluster <- rep( 1:3, each=4 )

# estimate Copula model
mod1 <- sirt::rasch.copula2( dat=dat, itemcluster=itemcluster)

# estimate Rasch model
mod2 <- sirt::rasch.copula2( dat=dat, itemcluster=itemcluster,
        delta=rep(0,3), est.delta=rep(0,3 ) )
summary(mod1)
summary(mod2)

# }
# NOT RUN {
#############################################################################
# EXAMPLE 2: 11 items nested within 2 item clusters (testlets)
#    with 2 resp. 3 dependent and 6 independent items
#############################################################################

set.seed(5698)
I <- 11                             # number of items
n <- 3000                           # number of persons
b <- seq(-2,2, len=I)               # item difficulties
theta <- stats::rnorm( n, sd=1 ) # person abilities
# define item clusters
itemcluster <- rep(0,I)
itemcluster[ c(3,5 )] <- 1
itemcluster[c(2,4,9)] <- 2
# residual correlations
rho <- c( .7, .5 )

# simulate data
dat <- sirt::sim.rasch.dep( theta, b, itemcluster, rho )
colnames(dat) <- paste("I", seq(1,ncol(dat)), sep="")

# estimate Rasch copula model
mod1 <- sirt::rasch.copula2( dat, itemcluster=itemcluster )
summary(mod1)

# both item clusters have Cook-Johnson copula as dependency
mod1c <- sirt::rasch.copula2( dat, itemcluster=itemcluster,
            copula.type="cook.johnson")
summary(mod1c)

# first item boundary mixture and second item Cook-Johnson copula
mod1d <- sirt::rasch.copula2( dat, itemcluster=itemcluster,
            copula.type=c( "bound.mixt", "cook.johnson" ) )
summary(mod1d)

# compare result with Rasch model estimation in rasch.copula2
# delta must be set to zero
mod2 <- sirt::rasch.copula2( dat, itemcluster=itemcluster, delta=c(0,0),
            est.delta=c(0,0) )
summary(mod2)

#############################################################################
# EXAMPLE 3: 12 items nested within 3 item clusters (testlets)
#   Cluster 1 -> Items 1-4; Cluster 2 -> Items 6-9;  Cluster 3 -> Items 10-12
#############################################################################

set.seed(967)
I <- 12                             # number of items
n <- 450                            # number of persons
b <- seq(-2,2, len=I)               # item difficulties
b <- sample(b)                      # sample item difficulties
theta <- stats::rnorm( n, sd=1 ) # person abilities
# itemcluster
itemcluster <- rep(0,I)
itemcluster[ 1:4 ] <- 1
itemcluster[ 6:9 ] <- 2
itemcluster[ 10:12 ] <- 3
# residual correlations
rho <- c( .35, .25, .30 )

# simulate data
dat <- sirt::sim.rasch.dep( theta, b, itemcluster, rho )
colnames(dat) <- paste("I", seq(1,ncol(dat)), sep="")

# estimate Rasch copula model
mod1 <- sirt::rasch.copula2( dat, itemcluster=itemcluster )
summary(mod1)

# person parameter estimation assuming the Rasch copula model
pmod1 <- sirt::person.parameter.rasch.copula(raschcopula.object=mod1 )

# Rasch model estimation
mod2 <- sirt::rasch.copula2( dat, itemcluster=itemcluster,
             delta=rep(0,3), est.delta=rep(0,3) )
summary(mod1)
summary(mod2)

#############################################################################
# EXAMPLE 4: Two-dimensional copula model
#############################################################################

set.seed(5698)
I <- 9
n <- 1500                           # number of persons
b <- seq(-2,2, len=I)               # item difficulties
theta0 <- stats::rnorm( n, sd=sqrt( .6 ) )

#*** Dimension 1
theta <- theta0 + stats::rnorm( n, sd=sqrt( .4 ) )   # person abilities
# itemcluster
itemcluster <- rep(0,I)
itemcluster[ c(3,5 )] <- 1
itemcluster[c(2,4,9)] <- 2
itemcluster1 <- itemcluster
# residual correlations
rho <- c( .7, .5 )
# simulate data
dat <- sirt::sim.rasch.dep( theta, b, itemcluster, rho )
colnames(dat) <- paste("A", seq(1,ncol(dat)), sep="")
dat1 <- dat
# estimate model of dimension 1
mod0a <- sirt::rasch.copula2( dat1, itemcluster=itemcluster1)
summary(mod0a)

#*** Dimension 2
theta <- theta0 + stats::rnorm( n, sd=sqrt( .8 ) )        # person abilities
# itemcluster
itemcluster <- rep(0,I)
itemcluster[ c(3,7,8 )] <- 1
itemcluster[c(4,6)] <- 2
itemcluster2 <- itemcluster
# residual correlations
rho <- c( .2, .4 )
# simulate data
dat <- sirt::sim.rasch.dep( theta, b, itemcluster, rho )
colnames(dat) <- paste("B", seq(1,ncol(dat)), sep="")
dat2 <- dat
# estimate model of dimension 2
mod0b <- sirt::rasch.copula2( dat2, itemcluster=itemcluster2)
summary(mod0b)

# both dimensions
dat <- cbind( dat1, dat2 )
itemcluster2 <- ifelse( itemcluster2 > 0, itemcluster2 +2, 0 )
itemcluster <- c( itemcluster1, itemcluster2 )
dims <- rep( 1:2, each=I)

# estimate two-dimensional copula model
mod1 <- sirt::rasch.copula3( dat, itemcluster=itemcluster, dims=dims, est.a=dims,
            theta.k=seq(-5,5,len=15) )
summary(mod1)

#############################################################################
# EXAMPLE 5: Subset of data Example 2
#############################################################################

set.seed(5698)
I <- 11                             # number of items
n <- 3000                           # number of persons
b <- seq(-2,2, len=I)               # item difficulties
theta <- stats::rnorm( n, sd=1.3 )  # person abilities
# define item clusters
itemcluster <- rep(0,I)
itemcluster[ c(3,5)] <- 1
itemcluster[c(2,4,9)] <- 2
# residual correlations
rho <- c( .7, .5 )
# simulate data
dat <- sirt::sim.rasch.dep( theta, b, itemcluster, rho )
colnames(dat) <- paste("I", seq(1,ncol(dat)), sep="")

# select subdataset with only one dependent item cluster
item.sel <- scan( what="character", nlines=1 )
    I1 I6 I7 I8 I10 I11 I3 I5
dat1 <- dat[,item.sel]

#******************
#*** Model 1a: estimate Copula model in sirt
itemcluster <- rep(0,8)
itemcluster[c(7,8)] <- 1
mod1a <- sirt::rasch.copula2( dat3, itemcluster=itemcluster )
summary(mod1a)

#******************
#*** Model 1b: estimate Copula model in mirt
library(mirt)
#*** redefine dataset for estimation in mirt
dat2 <- dat1[, itemcluster==0 ]
dat2 <- as.data.frame(dat2)
# combine items 3 and 5
dat2$C35 <- dat1[,"I3"] + 2*dat1[,"I5"]
table( dat2$C35, paste0( dat1[,"I3"],dat1[,"I5"]) )
#* define mirt model
mirtmodel <- mirt::mirt.model("
      F=1-7
      CONSTRAIN=(1-7,a1)
      " )
#-- Copula function with two dependent items
# define item category function for pseudo-items like C35
P.copula2 <- function(par,Theta, ncat){
     b1 <- par[1]
     b2 <- par[2]
     a1 <- par[3]
     ldelta <- par[4]
     P1 <- stats::plogis( a1*(Theta - b1 ) )
     P2 <- stats::plogis( a1*(Theta - b2 ) )
     Q1 <- 1-P1
     Q2 <- 1-P2
     # define vector-wise minimum function
     minf2 <- function( x1, x2 ){
         ifelse( x1 < x2, x1, x2 )
                                }
     # distribution under independence
     F00 <- Q1*Q2
     F10 <- Q1*Q2 + P1*Q2
     F01 <- Q1*Q2 + Q1*P2
     F11 <- 1+0*Q1
     F_ind <- c(F00,F10,F01,F11)
     # distribution under maximal dependence
     F00 <- minf2(Q1,Q2)
     F10 <- Q2              #=minf2(1,Q2)
     F01 <- Q1              #=minf2(Q1,1)
     F11 <- 1+0*Q1          #=minf2(1,1)
     F_dep <- c(F00,F10,F01,F11)
     # compute mixture distribution
     delta <- stats::plogis(ldelta)
     F_tot <- (1-delta)*F_ind + delta * F_dep
     # recalculate probabilities of mixture distribution
     L1 <- length(Q1)
     v1 <- 1:L1
     F00 <- F_tot[v1]
     F10 <- F_tot[v1+L1]
     F01 <- F_tot[v1+2*L1]
     F11 <- F_tot[v1+3*L1]
     P00 <- F00
     P10 <- F10 - F00
     P01 <- F01 - F00
     P11 <- 1 - F10 - F01 + F00
     prob_tot <- c( P00, P10, P01, P11 )
     return(prob_tot)
        }
# create item
copula2 <- mirt::createItem(name="copula2", par=c(b1=0, b2=0.2, a1=1, ldelta=0),
                est=c(TRUE,TRUE,TRUE,TRUE), P=P.copula2,
                lbound=c(-Inf,-Inf,0,-Inf), ubound=c(Inf,Inf,Inf,Inf) )
# define item types
itemtype <- c( rep("2PL",6), "copula2" )
customItems <- list("copula2"=copula2)
# parameter table
mod.pars <- mirt::mirt(dat2, 1, itemtype=itemtype,
                customItems=customItems, pars='values')
# estimate model
mod1b <- mirt::mirt(dat2, mirtmodel, itemtype=itemtype, customItems=customItems,
                verbose=TRUE, pars=mod.pars,
                technical=list(customTheta=as.matrix(seq(-4,4,len=21)) ) )
# estimated coefficients
cmod <- sirt::mirt.wrapper.coef(mod)$coef

# compare common item discrimination
round( c("sirt"=mod1a$item$a[1], "mirt"=cmod$a1[1] ), 4 )
  ##     sirt   mirt
  ##   1.2845 1.2862
# compare delta parameter
round( c("sirt"=mod1a$item$delta[7], "mirt"=stats::plogis( cmod$ldelta[7] ) ), 4 )
  ##     sirt   mirt
  ##   0.6298 0.6297
# compare thresholds a*b
dfr <- cbind( "sirt"=mod1a$item$thresh,
               "mirt"=c(- cmod$d[-7],cmod$b1[7]*cmod$a1[1], cmod$b2[7]*cmod$a1[1]))
round(dfr,4)
  ##           sirt    mirt
  ##   [1,] -1.9236 -1.9231
  ##   [2,] -0.0565 -0.0562
  ##   [3,]  0.3993  0.3996
  ##   [4,]  0.8058  0.8061
  ##   [5,]  1.5293  1.5295
  ##   [6,]  1.9569  1.9572
  ##   [7,] -1.1414 -1.1404
  ##   [8,] -0.4005 -0.3996
# }

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning