slca: Structured Latent Class Analysis (SLCA)

Description

This function implements a structured latent class model for polytomous item responses (Formann, 1992).

Usage

slca(data, group = NULL, weights = rep(1, nrow(data)), Xdes, 
  Xlambda.init = NULL, Xlambda.fixed = NULL, Xlambda.constr.V=NULL , 
  Xlambda.constr.c=NULL ,  delta.designmatrix = NULL, 
  delta.init = NULL, delta.fixed = NULL, delta.linkfct = "log" ,
  maxiter = 1000, conv = 10^(-5), globconv = 10^(-5), msteps = 4, convM = 5e-04, 
  decrease.increments = FALSE, oldfac=0, seed=NULL , progress=TRUE , ...)
  
## S3 method for class 'slca':
summary(object, file = NULL , \dots)

## S3 method for class 'slca':
print(x, \dots)    

## S3 method for class 'slca':
plot(x, group=1 , ... )

Arguments

data

Matrix of polytomous item responses

group

Optional vector of group identifiers. For plot.slca it is a single integer group identified.

weights

Optional vector of sample weights

Xdes

Design matrix for $x_{ijh}$ with $q_{ihjv}$ entries. Therefore, it must be an array with four dimensions referring to items ($i$), categories ($h$), latent classes ($j$) and $\lambda$ parameters ($v$).

Xlambda.init

Initial $\lambda_x$ parameters

Xlambda.fixed

Fixed $\lambda_x$ parameters. These must be provided by a matrix with two columns: 1st column -- Parameter index, 2nd column: Fixed value.

Xlambda.constr.V

A design matrix for linear restrictions of the form $V_x \lambda_x = c_x$ for the $\lambda_x$ parameter.

Xlambda.constr.c

A vector for the linear restriction $V_x \lambda_x = c_x$ of the $\lambda_x$ parameter.

delta.designmatrix

Design matrix for delta parameters $\delta$ parametrizing the latent class distribution by log-linear smoothing (Xu & von Davier, 2008)

delta.init

Initial $\delta$ parameters

delta.fixed

Fixed $\delta$ parameters. This must be a matrix with three columns: 1st column: Parameter index, 2nd column: Group index, 3rd column: Fixed value

delta.linkfct

Link function for skill space reduction. This can be the log-linear link (log) or the logistic link function (logit).

maxiter

Maximum number of iterations

conv

Convergence criterion for item parameters and distribution parameters

globconv

Global deviance convergence criterion

msteps

Maximum number of M steps in estimating $b$ and $a$ item parameters. The default is to use 4 M steps.

convM

Convergence criterion in M step

decrease.increments

Should in the M step the increments of $a$ and $b$ parameters decrease during iterations? The default is FALSE. If there is an increase in deviance during estimation, setting decrease.increments to TRUE

oldfac

Factor $f$ between 0 and 1 to control convergence behavior. If $x_t$ denotes the estimated parameter in iteration $t$, then the regularized estimate $x_t^{\ast}$ is obtained by $x_t^{\ast} = f x_{t-1} + (1-f) x_t$. Therefore, values of oldfa

seed

Simulation seed for initial parameters. The default of NULL corresponds to a random seed.

progress

An optional logical indicating whether the function should print the progress of iteration in the estimation process.

object

A required object of class slca

file

Optional file name for a file in which summary should be sinked.

A required object of class slca

...

Optional parameters to be passed to or from other methods will be ignored.

Value

An object of class slca. The list contains the following entries:
itemData frame with conditional item probabilities
devianceDeviance
icInformation criteria, number of estimated parameters
XlambdaEstimated $\lambda_x$ parameters
se.XlambdaStandard error of $\lambda_x$ parameters
pi.kTrait distribution
pjkItem response probabilities evaluated for all classes
n.ikAn array of expected counts $n_{cikg}$ of ability class $c$ at item $i$ at category $k$ in group $g$
GNumber of groups
INumber of items
NNumber of persons
deltaParameter estimates for skillspace representation
covdeltaCovariance matrix of parameter estimates for skillspace representation
MLE.classClassified skills for each student (MLE)
MAP.classClassified skills for each student (MAP)
dataOriginal data frame
group.statGroup statistics (sample sizes, group labels)
p.xi.ajIndividual likelihood
posteriorIndividual posterior distribution
K.itemMaximal category per item
timeInfo about computation time
skillspaceUsed skillspace parametrization
iterNumber of iterations
seed.usedUsed simulation seed
Xlambda.initUsed initial lambda parameters
delta.initUsed initial delta parameters
convergedLogical indicating whether convergence was achieved.

Details

The structured latent class model allows for general constraints of items $i$ in categories $h$ and classes $j$. The item response model is $$P( X_{i} = h | j ) = \frac{ \exp( x_{ihj} ) }{ \sum_l \exp( x_{ilj} ) }$$ with linear constraints on the class specific probabilities $$x_{ihj} = \sum_v q_{ihjv} \lambda_{xv}$$ Linear restrictions on the $\lambda_x$ parameter can be specified by a matrix equation $V_x \lambda_x = c_x$ (see Xlambda.constr.V and Xlambda.constr.c; Neuhaus, 1996). The latent class distribution can be smoothed by a log-linear link function (Xu & von Davier, 2008) or a logistic link function (Formann, 1992). For class $j$ in group $g$ employing a link function $h$, it holds that $$h [ P( j| g) ] \propto \sum_w r_{jw} \delta_{gw}$$ where group-specific distributions are allowed. The values $r_{jw}$ are specified in the design matrix delta.designmatrix. This model contains classical uni- and multidimensional latent trait models, latent class analysis, located latent class analysis, cognitive diagnostic models, the general diagnostic model and mixture item response models as special cases (see Formann & Kohlmann, 1998; Formann, 2007).

References

Formann, A. K. (1992). Linear logistic latent class analysis for polytomous data. Journal of the American Statistical Association, 87, 476-486. Formann, A. K. (2007). (Almost) Equivalence between conditional and mixture maximum likelihood estimates for some models of the Rasch type. In M. von Davier & C. H. Carstensen (Eds.), Multivariate and mixture distribution Rasch models (pp. 177-189). New York: Springer. Formann, A. K., & Kohlmann, T. (1998). Structural latent class models. Sociological Methods & Research, 26, 530-565. Neuhaus, W. (1996). Optimal estimation under linear constraints. Astin Bulletin, 26, 233-245. Xu, X., & von Davier, M. (2008). Fitting the structured general diagnostic model to NAEP data. ETS Research Report ETS RR-08-27. Princeton, ETS.

Examples

Run this code

#############################################################################
# EXAMPLE 1: data.Students | (Generalized) Partial Credit Model
#############################################################################

data(data.Students)
dat <- data.Students[ , c("mj1","mj2","mj3","mj4","sc1" , "sc2") ]
# define discretized ability
theta.k <- seq( -6 , 6 , len=21 )

#*** Model 1: Partial credit model

# define design matrix for lambda
I <- ncol(dat)
maxK <- 4
TP <- length(theta.k)
NXlam <- I*(maxK-1) + 1       # number of estimated parameters
       # last parameter is joint slope parameter
Xdes <- array( 0 , dim=c(I , maxK , TP ,  NXlam ) )
# Item1Cat1 , ... , Item1Cat3, Item2Cat1, ..., 
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[3]] <- paste0("Class" , 1:TP )
v2 <- unlist( sapply( 1:I , FUN = function(ii){ # ii
    paste0( paste0( colnames(dat)[ii]  , "_b"  ) , "Cat" , 1:(maxK-1) )
                } , simplify=FALSE) )                
dimnames(Xdes)[[4]] <- c( v2 , "a" )
# define theta design and item discriminations
for (ii in 1:I){
for (hh in 1:(maxK-1) ){
    Xdes[ii , hh + 1 , , NXlam ] <- hh * theta.k
                    }
            }
# item intercepts
for (ii in 1:I){
for (hh in 1:(maxK-1) ){
    # ii <- 1  # Item    # hh <- 1  # category
    Xdes[ii,hh+1,, ( ii - 1)*(maxK-1) + hh] <- 1
                        }
                    }
#****
# skill space designmatrix
TP <- length(theta.k)
w1 <- stats::dnorm(theta.k)
w1 <- w1 / sum(w1)  
delta.designmatrix <- matrix( 1 , nrow=TP , ncol=1 )
delta.designmatrix[,1] <- log(w1)

# initial lambda parameters
Xlambda.init <- c( stats::rnorm( dim(Xdes)[[4]] - 1 ) , 1 )
# fixed delta parameter
delta.fixed <- cbind( 1 , 1 ,1 )

# estimate model
mod1 <- slca( dat , Xdes=Xdes , delta.designmatrix=delta.designmatrix , 
            Xlambda.init=Xlambda.init , delta.fixed = delta.fixed, maxiter = 70 )
summary(mod1)
plot(mod1, cex.names=.7 )

#*** Model 2: Partial credit model with some parameter constraints
# fixed lambda parameters
Xlambda.fixed <- cbind( c(1,19) , c(3.2,1.52 ) )
# 1st parameter = 3.2
# 19th parameter = 1.52 (joint item slope)
mod2 <- slca( dat , Xdes=Xdes , delta.designmatrix=delta.designmatrix , 
            delta.init=delta.init, Xlambda.init=Xlambda.init, delta.fixed=delta.fixed,
            Xlambda.fixed=Xlambda.fixed , maxiter = 70 )

#*** Model 3: Partial credit model with non-normal distribution
Xlambda.fixed <- cbind(  c(1,19) , c(3.2,1) ) # fix item slope to one
delta.designmatrix <- cbind( 1 , theta.k , theta.k^2 , theta.k^3 )
mod3 <- slca( dat , Xdes=Xdes,  delta.designmatrix=delta.designmatrix ,          
            Xlambda.fixed=Xlambda.fixed ,  maxiter = 200 )
summary(mod3)  

# non-normal distribution with convergence regularizing factor oldfac
mod3a <- slca( dat , Xdes=Xdes,  delta.designmatrix=delta.designmatrix ,          
            Xlambda.fixed=Xlambda.fixed, maxiter = 500 , oldfac=.95 )
summary(mod3a)  

#*** Model 4: Generalized Partial Credit Model

# estimate generalized partial credit model without restrictions on trait
# distribution and item parameters to ensure better convergence behavior
# Note that two parameters are not identifiable and information criteria
# have to be adapted.

#---
# define design matrix for lambda
I <- ncol(dat)
maxK <- 4
TP <- length(theta.k)
NXlam <- I*(maxK-1) + I       # number of estimated parameters
Xdes <- array( 0 , dim=c(I , maxK , TP ,  NXlam ) )
# Item1Cat1 , ... , Item1Cat3, Item2Cat1, ..., 
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[3]] <- paste0("Class" , 1:TP )
v2 <- unlist( sapply( 1:I , FUN = function(ii){ # ii
    paste0( paste0( colnames(dat)[ii]  , "_b"  ) , "Cat" , 1:(maxK-1) )
                } , simplify=FALSE) )                
dimnames(Xdes)[[4]] <- c( v2 , paste0( colnames(dat),"_a") )
dimnames(Xdes)
# define theta design and item discriminations
for (ii in 1:I){
for (hh in 1:(maxK-1) ){
    Xdes[ii , hh + 1 , , I*(maxK-1) + ii ] <- hh * theta.k
                    }
            }
# item intercepts
for (ii in 1:I){
for (hh in 1:(maxK-1) ){
    Xdes[ii,hh+1,, ( ii - 1)*(maxK-1) + hh] <- 1
                        }
                    }
#****
# skill space designmatrix
delta.designmatrix <- cbind( 1, theta.k,theta.k^2 )
# initial lambda parameters from partial credit model
Xlambda.init <- mod1$Xlambda
Xlambda.init <- c( mod1$Xlambda[ - length(Xlambda.init) ] , 
         rep( Xlambda.init[ length(Xlambda.init)  ] ,I) )

# estimate model
mod4 <- slca( dat , Xdes=Xdes , Xlambda.init=Xlambda.init , 
             delta.designmatrix=delta.designmatrix , decrease.increments=TRUE ,         
             maxiter = 300 )           

#############################################################################
# SIMULATED EXAMPLE 2: Latent class model with two classes
#############################################################################

set.seed(9876)
I <- 7	# number of items
# simulate response probabilities
a1 <- stats::runif(I, 0 , .4 )
a2 <- stats::runif(I, .6 , 1 )
N <- 1000	# sample size
# simulate data in two classes of proportions .3 and .7
N1 <- round(.3*N)
dat1 <- 1 * ( matrix(a1,N1,I,byrow=TRUE) > matrix( stats::runif( N1 * I) , N1 , I ) )
N2 <- round(.7*N)
dat2 <- 1 * ( matrix(a2,N2,I,byrow=TRUE) > matrix( stats::runif( N2 * I) , N2 , I ) )
dat <- rbind( dat1 , dat2 )
colnames(dat) <- paste0("I" , 1:I)

# define design matrices
TP <- 2     # two classes
# The idea is that latent classes refer to two different "dimensions".
# Items load on latent class indicators 1 and 2, see below.
Xdes <- array(0 , dim=c(I,2,2,2*I) )
items <- colnames(dat)
dimnames(Xdes)[[4]] <- c(paste0( colnames(dat) , "Class" , 1),
          paste0( colnames(dat) , "Class" , 2) )
    # items, categories , classes , parameters
# probabilities for correct solution
for (ii in 1:I){
    Xdes[ ii , 2 , 1 , ii ] <- 1    # probabilities class 1
    Xdes[ ii , 2 , 2 , ii+I ] <- 1  # probabilities class 2
                    }
# estimate model
mod1 <- slca( dat , Xdes=Xdes )            
summary(mod1)

#############################################################################
# SIMULATED EXAMPLE 3: Mixed Rasch model with two classes
#############################################################################

set.seed(987)
library(sirt)
# simulate two latent classes of Rasch populations
I <- 15  # 6 items
b1 <- seq( -1.5 , 1.5 , len=I)      # difficulties latent class 1
b2 <- b1    # difficulties latent class 2
b2[ c(4,7, 9 , 11 , 12, 13) ] <- c(1 , -.5 , -.5 , .33 , .33 , -.66 )
N <- 3000    # number of persons
wgt <- .25       # class probability for class 1
# class 1
dat1 <- sirt::sim.raschtype( stats::rnorm( wgt*N ) , b1 )
# class 2
dat2 <- sirt::sim.raschtype( stats::rnorm( (1-wgt)*N , mean= 1 , sd=1.7) , b2 )
dat <- rbind( dat1 , dat2 )
# theta grid
theta.k <- seq( -5 , 5 , len=9 )
TP <- length(theta.k)

#*** Model 1: Rasch model with normal distribution
maxK <- 2
NXlam <- I +1
Xdes <- array( 0 , dim=c(I , maxK , TP ,  NXlam ) )
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[4]] <- c( paste0( "b_" , colnames(dat)[1:I] ) , "a" )
# define item difficulties
for (ii in 1:I){  Xdes[ii , 2 , , ii ] <- -1        }
# theta design
for (tt in 1:TP){  Xdes[1:I , 2 , tt , I + 1] <- theta.k[tt]       }

# skill space definition
delta.designmatrix <- cbind( 1 , theta.k^2 )
delta.fixed <- NULL
Xlambda.init <- c( stats::runif( I  , -.8 , .8 ) , 1 )
Xlambda.fixed <- cbind( I+1 , 1 )
# estimate model
mod1 <- slca( dat , Xdes=Xdes , delta.designmatrix=delta.designmatrix , 
            delta.fixed=delta.fixed , Xlambda.fixed=Xlambda.fixed , 
            Xlambda.init=Xlambda.init, decrease.increments=TRUE , maxiter = 200 )
summary(mod1)

#*** Model 1b: Constraint the sum of item difficulties to zero

# change skill space definition
delta.designmatrix <- cbind( 1 , theta.k , theta.k^2 )
delta.fixed <- NULL
# constrain sum of difficulties Xlambda parameters to zero
Xlambda.constr.V <- matrix( 1 , nrow=I+1 , ncol=1 )
Xlambda.constr.V[I+1,1] <- 0
Xlambda.constr.c <- c(0)
# estimate model
mod1b <- slca( dat , Xdes=Xdes , delta.designmatrix=delta.designmatrix ,         
            Xlambda.fixed=Xlambda.fixed , Xlambda.constr.V=Xlambda.constr.V ,
            Xlambda.constr.c=Xlambda.constr.c  )
summary(mod1b)

#*** Model 2: Mixed Rasch model with two latent classes
NXlam <- 2*I +2
Xdes <- array( 0 , dim=c(I , maxK , 2*TP ,  NXlam ) )
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[4]] <- c( paste0( "bClass1_" , colnames(dat)[1:I] ) , 
        paste0( "bClass2_" , colnames(dat)[1:I] ) , "aClass1" , "aClass2" )
# define item difficulties
for (ii in 1:I){
    Xdes[ii , 2 , 1:TP , ii ] <- -1  # first class
    Xdes[ii , 2 , TP + 1:TP , I+ii ] <- -1  # second class    
                    }                                                                               
# theta design
for (tt in 1:TP){
    Xdes[1:I , 2 , tt , 2*I+1 ] <- theta.k[tt]
    Xdes[1:I , 2 , TP+tt , 2*I+2 ] <- theta.k[tt]
                    }
# skill space definition
delta.designmatrix <- matrix( 0 , nrow=2*TP , ncol=4 )
delta.designmatrix[1:TP,1] <- 1
delta.designmatrix[1:TP,2] <- theta.k^2
delta.designmatrix[TP + 1:TP,3] <- 1
delta.designmatrix[TP+ 1:TP,4] <- theta.k^2
b1 <- stats::qnorm( colMeans( dat ) )
Xlambda.init <- c( b1 , b1 + stats::runif(I, -2  , 2 ) , 1 , 1 )
Xlambda.init <- c( stats::runif( 2*I  , -1.8 , 1.8 ) , 1 ,1 )
Xlambda.fixed <- cbind( c(2*I+1 , 2*I+2) , 1 )
# estimate model
mod2 <- slca( dat, Xdes=Xdes,  delta.designmatrix=delta.designmatrix ,          
            Xlambda.fixed=Xlambda.fixed, decrease.increments=TRUE ,
            Xlambda.init=Xlambda.init, maxiter = 1000 )
summary(mod2)
summary(mod1)
# latent class proportions
stats::aggregate( mod2$pi.k , list( rep(1:2, each=TP)) , sum )

#*** Model 2b: Different parametrization with sum constraint on item difficulties

# skill space definition
delta.designmatrix <- matrix( 0 , nrow=2*TP , ncol=6 )
delta.designmatrix[1:TP,1] <- 1
delta.designmatrix[1:TP,2] <- theta.k
delta.designmatrix[1:TP,3] <- theta.k^2
delta.designmatrix[TP+ 1:TP,4] <- 1
delta.designmatrix[TP+ 1:TP,5] <- theta.k
delta.designmatrix[TP+ 1:TP,6] <- theta.k^2
Xlambda.fixed <- cbind( c(2*I+1,2*I+2) , c(1,1) )
b1 <- stats::qnorm( colMeans( dat ) )
Xlambda.init <- c( b1 , b1 + stats::runif(I, -1  , 1 ) , 1 , 1 )
# constraints on item difficulties
Xlambda.constr.V <- matrix( 0 , nrow=NXlam , ncol=2)
Xlambda.constr.V[1:I , 1 ] <- 1
Xlambda.constr.V[I + 1:I , 2 ] <- 1
Xlambda.constr.c <- c(0,0)
# estimate model
mod2b <- slca( dat, Xdes=Xdes, delta.designmatrix=delta.designmatrix ,          
            Xlambda.fixed=Xlambda.fixed ,  Xlambda.init=Xlambda.init ,  
            Xlambda.constr.V=Xlambda.constr.V , Xlambda.constr.c=Xlambda.constr.c , 
            decrease.increments=TRUE , maxiter = 1000 )
summary(mod2b)
stats::aggregate( mod2b$pi.k , list( rep(1:2, each=TP)) , sum )

#*** Model 2c: Estimation with mRm package
library(mRm)
mod2c <- mRm::mrm(data.matrix=dat, cl=2)
plot(mod2c)
print(mod2c)

#*** Model 2d: Estimation with psychomix package
library(psychomix)
mod2d <- psychomix::raschmix(data=dat, k=2 , verbose=TRUE )
summary(mod2d)
plot(mod2d)

#############################################################################
# SIMULATED EXAMPLE 4: Located latent class model, Rasch model
#############################################################################

set.seed(487)
library(sirt)
I <- 15  # I items
b1 <- seq( -2 , 2 , len=I)      # item difficulties
N <- 4000    # number of persons
# simulate 4 theta classes
theta0 <- c( -2.5 , -1 , 0.3 , 1.3 )  # skill classes
probs0 <- c( .1 , .4 , .2 , .3 )
TP <- length(theta0)
theta <- theta0[ rep(1:TP, round(probs0*N)  ) ]
dat <- sirt::sim.raschtype( theta , b1 )

#*** Model 1: Located latent class model with 4 classes
maxK <- 2
NXlam <- I + TP
Xdes <- array( 0 , dim=c(I , maxK , TP ,  NXlam ) )
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[3]] <- paste0("Class" , 1:TP )
dimnames(Xdes)[[4]] <- c( paste0( "b_" , colnames(dat)[1:I] ) , paste0("theta" , 1:TP) )
# define item difficulties
for (ii in 1:I){  Xdes[ii , 2 , , ii ] <- -1        }
# theta design
for (tt in 1:TP){  Xdes[1:I , 2 , tt , I + tt] <- 1       }

# skill space definition
delta.designmatrix <- diag(TP)
Xlambda.init <- c( - stats::qnorm( colMeans(dat) ) , seq(-2,1,len=TP)  )
# constraint on item difficulties
Xlambda.constr.V <- matrix( 0 , nrow=NXlam , ncol=1)
Xlambda.constr.V[1:I,1] <- 1
Xlambda.constr.c <- c(0)
delta.init <- matrix( c(1,1,1,1) , TP , 1 )
# estimate model
mod1 <- slca( dat , Xdes=Xdes, delta.designmatrix=delta.designmatrix, 
            delta.init=delta.init, Xlambda.init=Xlambda.init , 
            Xlambda.constr.V=Xlambda.constr.V , Xlambda.constr.c=Xlambda.constr.c , 
            decrease.increments=TRUE,  maxiter = 400 )
summary(mod1)
# compare estimated and simulated theta class locations
cbind( mod1$Xlambda[ - c(1:I) ] , theta0 )
# compare estimated and simulated latent class proportions
cbind( mod1$pi.k , probs0 )

#############################################################################
# SIMULATED EXAMPLE 5: DINA model with two skills
#############################################################################

set.seed(487)
N <- 3000   # number of persons
# define Q-matrix
I <- 9  # 9 items
NS <- 2 # 2 skills 
TP <- 4 # number of skill classes
Q <- scan( nlines=3)
  1 0   1 0   1 0
  0 1   0 1   0 1
  1 1   1 1   1 1
Q <- matrix(Q , I, ncol=NS,byrow=TRUE)
# define skill distribution
alpha0 <- matrix( c(0,0,1,0,0,1,1,1) , nrow=4,ncol=2,byrow=TRUE)
prob0 <- c( .2 , .4 , .1 , .3 )
alpha <- alpha0[ rep( 1:TP , prob0*N) ,]
# define guessing and slipping parameters
guess <- round( stats::runif(I, 0 , .4 ) , 2 )
slip <- round( stats::runif(I , 0 , .3 ) , 2 )
# simulate data according to the DINA model
dat <- CDM::sim.din( q.matrix=Q , alpha=alpha , slip=slip , guess=guess )$dat

# define Xlambda design matrix
maxK <- 2
NXlam <- 2*I
Xdes <- array( 0 , dim=c(I , maxK , TP ,  NXlam ) )
dimnames(Xdes)[[1]] <- colnames(dat)
dimnames(Xdes)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(Xdes)[[3]] <- c("S00","S10","S01","S11")
dimnames(Xdes)[[4]] <- c( paste0("guess",1:I ) , paste0( "antislip" , 1:I ) )
dimnames(Xdes)
# define item difficulties
for (ii in 1:I){  
        # define latent responses
        latresp <- 1*( alpha0 %*% Q[ii,]  == sum(Q[ii,]) )[,1]
        # model slipping parameters
        Xdes[ii , 2 , latresp == 1, I+ii ] <- 1
        # guessing parameters
        Xdes[ii , 2 , latresp == 0, ii ] <- 1        
                 }
Xdes[1,2,,] ; Xdes[7,2,,]
# skill space definition
delta.designmatrix <- diag(TP)
Xlambda.init <- c( rep( stats::qlogis( .2 ) , I ) , rep( stats::qlogis( .8 ) , I ) )

# estimate DINA model with slca function
mod1 <- slca( dat , Xdes=Xdes, delta.designmatrix=delta.designmatrix , 
            Xlambda.init=Xlambda.init, decrease.increments=TRUE, maxiter=400 )
summary(mod1)

# compare estimated and simulated latent class proportions
cbind( mod1$pi.k , probs0 )
# compare estimated and simulated guessing parameters
cbind( mod1$pjk[1,,2] , guess )
# compare estimated and simulated slipping parameters
cbind( 1 - mod1$pjk[4,,2] , slip )

Run the code above in your browser using DataLab