#############################################################################
# EXAMPLE 1: Dichotomous data | sim.rasch
#############################################################################
data(sim.rasch)
dat <- sim.rasch
# some control arguments
ctl.list <- list(maxiter= 20 ) # increase the number of iterations in applications!
#*** Model 1: Rasch model, normal trait distribution
mod1 <- tam.mml.3pl(resp=dat , skillspace= "normal" , est.some.slopes=FALSE ,
control=ctl.list)
summary(mod1)
#*** Model 2: Rasch model, discrete trait distribution
# choose theta grid
theta.k <- seq( -3 , 3 , len=7 ) # discrete theta grid distribution
# define symmetric trait distribution
delta.designmatrix <- matrix( 0 , nrow=7 , ncol= 4 )
delta.designmatrix[4,1] <- 1
delta.designmatrix[c(3,5),2] <- 1
delta.designmatrix[c(2,6),3] <- 1
delta.designmatrix[c(1,7),4] <- 1
mod2 <- tam.mml.3pl(resp=dat , skillspace= "discrete" , est.some.slopes=FALSE ,
theta.k=theta.k, delta.designmatrix=delta.designmatrix , control=ctl.list)
summary(mod2)
#*** Model 3: 2PL model
mod3 <- tam.mml.3pl(resp=dat , skillspace= "normal" , gammaslope.des="2PL" ,
control=ctl.list, est.variance=FALSE )
summary(mod3)
#*** Model 4: 3PL model
# estimate guessing parameters for items 3,7,9 and 12
I <- ncol(dat)
est.guess <- rep(0, I )
# set parameters 9 and 12 equal -> same integers
est.guess[ c(3,7,9,12) ] <- c( 1 , 3 , 2 , 2 )
# starting values guessing parameter
guess <- .2*(est.guess > 0)
# estimate model
mod4 <- tam.mml.3pl(resp=dat , skillspace= "normal" , gammaslope.des="2PL" ,
control=ctl.list , est.guess = est.guess , guess=guess, est.variance=FALSE)
summary(mod4)
#--- specification in tamaan
tammodel <- "
LAVAAN MODEL:
F1 =~ I1__I40
F1 ~~ 1*F1
I3 + I7 ?= g1
I9 + I12 ?= g912 * g1
"
mod4a <- tamaan( tammodel , resp=dat , control=list(maxiter=20))
summary(mod4a)
#*** Model 5: 3PL model, add some prior Beta distribution
guess.prior <- matrix( 0 , nrow=I , ncol=2 )
guess.prior[ est.guess > 0 , 1] <- 5
guess.prior[ est.guess > 0 , 2] <- 17
mod5 <- tam.mml.3pl(resp=dat , skillspace= "normal" , gammaslope.des="2PL" ,
control=ctl.list , est.guess = est.guess , guess=guess , guess.prior=guess.prior)
summary(mod5)
#--- specification in tamaan
tammodel <- "
LAVAAN MODEL:
F1 =~ I1__I40
F1 ~~ 1*F1
I3 + I7 ?= g1
I9 + I12 ?= g912 * g1
MODEL PRIOR:
g912 ~ Beta(5,17)
I3_guess ~ Beta(5,17)
I7_guess ~ Beta(5,17)
"
mod5a <- tamaan( tammodel , resp=dat , control=list(maxiter=20))
#*** Model 6: 2PL model with design matrix for item slopes
I <- 40 # number of items
D <- 1 # dimensions
maxK <- 2 # maximum number of categories
Ngam <- 13 # number of different slope parameters
E <- array( 0 , dim=c(I,maxK,D,Ngam) )
# joint slope parameters for items 1 to 10, 11 to 20, 21 to 30
E[ 1:10 , 2 , 1 , 2 ] <- 1
E[ 11:20 , 2 , 1 , 1 ] <- 1
E[ 21:30 , 2 , 1 , 3 ] <- 1
for (ii in 31:40){ E[ii,2,1,ii - 27 ] <- 1 }
# estimate model
mod6 <- tam.mml.3pl(resp= dat , control= ctl.list , E=E , est.variance=FALSE )
summary(mod6)
#*** Model 6b: Truncated normal prior distribution for slope parameters
gammaslope.prior <- matrix( 0 , nrow = Ngam , ncol=4 )
gammaslope.prior[,1] <- 2 # mean
gammaslope.prior[,2] <- 10 # standard deviation
gammaslope.prior[,3] <- -Inf # lower bound
gammaslope.prior[ 4:13,3] <- 1.2
gammaslope.prior[,4] <- Inf # upper bound
# estimate model
mod6b <- tam.mml.3pl(resp= dat , E=E , est.variance=FALSE ,
gammaslope.prior=gammaslope.prior , control=ctl.list )
summary(mod6b)
#*** Model 7: 2PL model with design matrix of slopes and slope constraints
Ngam <- dim(E)[4] # number of gamma parameters
# define two constraint equations
gammaslope.constr.V <- matrix( 0 , nrow=Ngam , ncol=2 )
gammaslope.constr.c <- rep(0,2)
# set sum of first two xlambda entries to 1.8
gammaslope.constr.V[1:2,1] <- 1
gammaslope.constr.c[1] <- 1.8
# set sum of entries 4, 5 and 6 to 3
gammaslope.constr.V[4:6,2] <- 1
gammaslope.constr.c[2] <- 3
mod7 <- tam.mml.3pl(resp= dat, control= ctl.list , E=E , est.variance=FALSE ,
gammaslope.constr.V=gammaslope.constr.V , gammaslope.constr.c=gammaslope.constr.c)
summary(mod7)
#**** Model 8: Located latent class Rasch model with estimated three skill points
# three classes of theta's are estimated
TP <- 3
theta.k <- diag(TP)
# because item difficulties are unrestricted, we define the sum of the estimated
# theta points equal to zero
Ngam <- TP # estimate three gamma loading parameters which are discrete theta points
E <- array( 0 , dim=c(I,2,TP,Ngam) )
E[,2,1,1] <- E[,2,2,2] <- E[,2,3,3] <- 1
gammaslope.constr.V <- matrix( 1 , nrow=3 , ncol=1 )
gammaslope.constr.c <- c(0)
# initial gamma values
gammaslope <- c( -2 , 0 , 2 )
# estimate model
mod8 <- tam.mml.3pl(resp= dat, control= ctl.list , E=E , skillspace="discrete" ,
theta.k=theta.k , gammaslope=gammaslope , gammaslope.constr.V=gammaslope.constr.V ,
gammaslope.constr.c=gammaslope.constr.c )
summary(mod8)
#############################################################################
# EXAMPLE 2: Polytomous data
#############################################################################
data( data.mg , package="CDM")
dat <- data.mg[1:1000, paste0("I",1:11)]
#*******************************************************
#*** Model 1: 1-dimensional 1PL estimation, normal skill distribution
mod1 <- tam.mml.3pl(resp=dat , skillspace="normal" ,
control= list( maxiter= 30 ) , gammaslope.des = "2PL" ,
est.some.slopes=FALSE , est.variance=TRUE )
summary(mod1)
#*******************************************************
#*** Model 2: 1-dimensional 2PL estimation, discrete skill distribution
# define skill space
theta.k <- matrix( seq(-5,5,len=21) )
# allow skew skill distribution
delta.designmatrix <- cbind( 1 , theta.k , theta.k^2 , theta.k^3 )
# fix 13th xsi item parameter to zero
xsi.fixed <- cbind( 13 , 0 )
# fix 10th slope paremeter to one
gammaslope.fixed <- cbind( 10 , 1 )
# estimate model
mod2 <- tam.mml.3pl(resp=dat , skillspace="discrete" , theta.k=theta.k ,
delta.designmatrix=delta.designmatrix, control= list(maxiter=30),
gammaslope.des = "2PL" , xsi.fixed=xsi.fixed ,
gammaslope.fixed=gammaslope.fixed)
summary(mod2)
#*******************************************************
#*** Model 3: 2-dimensional 2PL estimation, normal skill distribution
# define loading matrix
Q <- matrix(0,11,2)
Q[1:6,1] <- 1 # items 1 to 6 load on dimension 1
Q[7:11,2] <- 1 # items 7 to 11 load on dimension 2
# estimate model
mod3 <- tam.mml.3pl(resp=dat , control= list(maxiter=30),
gammaslope.des = "2PL" , Q=Q )
summary(mod3)
#############################################################################
# SIMULATED EXAMPLE 3: Dichotomous data with guessing
#############################################################################
#*** simulate data
set.seed(9765)
N <- 4000 # number of persons
I <- 20 # number of items
b <- seq( -1.5 , 1.5 , len=I )
guess <- rep(0 , I )
guess.items <- c(6,11,16)
guess[ guess.items ] <- .33
library(sirt)
dat <- sirt::sim.raschtype( rnorm(N) , b=b , fixed.c=guess )
#*******************************************************
#*** Model 1: Difficulty + guessing model, i.e. fix slopes to 1
est.guess <- rep(0,I)
est.guess[ guess.items ] <- seq(1, length(guess.items) )
# define prior distribution
guess.prior <- matrix( cbind( 5 , 17 ) , I , 2 , byrow=TRUE )
guess.prior[ ! est.guess , ] <- 0
# estimate model
mod1 <- tam.mml.3pl(resp= dat , guess = guess , est.guess = est.guess ,
guess.prior=guess.prior , control= ctl.list ,est.variance=TRUE ,
est.some.slopes=FALSE )
summary(mod1)
#*******************************************************
#*** Model 2: estimate a joint guessing parameter
est.guess <- rep(0,I)
est.guess[ guess.items ] <- 1
# estimate model
mod2 <- tam.mml.3pl(resp= dat , guess = guess , est.guess = est.guess ,
guess.prior=guess.prior , control= ctl.list ,est.variance=TRUE ,
est.some.slopes=FALSE )
summary(mod2)
#############################################################################
# SIMULATED EXAMPLE 4: Latent class model with two classes
# See slca Simulated Example 2 in the CDM package
#############################################################################
#*******************************************************
#*** simulate data
set.seed(9876)
I <- 7 # number of items
# simulate response probabilities
a1 <- round(runif(I,0 , .4 ),4)
a2 <- round(runif(I , .6 , 1 ),4)
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( runif( N1 * I) , N1 , I ) )
N2 <- round(.7*N)
dat2 <- 1 * ( matrix(a2,N2,I,byrow=TRUE) > matrix( runif( N2 * I) , N2 , I ) )
dat <- rbind( dat1 , dat2 )
colnames(dat) <- paste0("I" , 1:I)
#*******************************************************
# estimation using tam.mml.3pl function
# define design matrices
TP <- 2 # two classes
theta.k <- diag(TP) # there are theta dimensions -> two classes
# The idea is that latent classes refer to two different "dimensions".
# Items load on latent class indicators 1 and 2, see below.
E <- array(0 , dim=c(I,2,2,2*I) )
items <- colnames(dat)
dimnames(E)[[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){
E[ ii , 2 , 1 , ii ] <- 1 # probabilities class 1
E[ ii , 2 , 2 , ii+I ] <- 1 # probabilities class 2
}
# estimation command
mod1 <- tam.mml.3pl(resp= dat , E=E , control= list(maxit=20) , skillspace="discrete" ,
theta.k=theta.k , notA=TRUE)
summary(mod1)
# compare simulated and estimated data
cbind( mod1$rprobs[,2,1] , a2 ) # Simulated class 2
cbind( mod1$rprobs[,2,2] , a1 ) # Simulated class 1
#*******************************************************
#** specification with tamaan
tammodel <- "
ANALYSIS:
TYPE=LCA;
NCLASSES(2); # 2 classes
NSTARTS(5,20); # 5 random starts with 20 iterations
LAVAAN MODEL:
F =~ I1__I7
"
mod1b <- tamaan( tammodel , resp=dat )
summary(mod1b)
# compare with mod1
logLik(mod1)
logLik(mod1b)
#############################################################################
# SIMULATED EXAMPLE 5: Located latent class model, Rasch model
# See slca Simulated Example 4 in the CDM package
#############################################################################
#*** simulate data
set.seed(487)
I <- 15 # I items
b1 <- seq( -2 , 2 , len=I) # item difficulties
N <- 2000 # number of persons
# simulate 4 theta classes
theta0 <- c( -2.5 , -1 , 0.3 , 1.3 ) # skill classes
probs0 <- c( .1 , .4 , .2 , .3 ) # skill class probabilities
TP <- length(theta0)
theta <- theta0[ rep(1:TP, round(probs0*N) ) ]
library(sirt)
dat <- sirt::sim.raschtype( theta , b1 )
colnames(dat) <- paste0("I",1:I)
#*******************************************************
#*** Model 1: Located latent class model with 4 classes
maxK <- 2
Ngam <- TP
E <- array( 0 , dim=c(I , maxK , TP , Ngam ) )
dimnames(E)[[1]] <- colnames(dat)
dimnames(E)[[2]] <- paste0("Cat", 1:(maxK) )
dimnames(E)[[3]] <- paste0("Class", 1:TP)
dimnames(E)[[4]] <- paste0("theta", 1:TP)
# theta design
for (tt in 1:TP){ E[1:I , 2 , tt , tt] <- 1 }
theta.k <- diag(TP)
# set eighth item difficulty to zero
xsi.fixed <- cbind( 8 , 0 )
# initial gamma parameter
gammaslope <- seq( -1.5 , 1.5 , len=TP)
# estimate model
mod1 <- tam.mml.3pl(resp= dat , E=E , xsi.fixed=xsi.fixed ,
control= list(maxiter=100) , skillspace="discrete" ,
theta.k=theta.k , gammaslope=gammaslope)
summary(mod1)
# compare estimated and simulated theta class locations
cbind( mod1$gammaslope , theta0 )
# compare estimated and simulated latent class proportions
cbind( mod1$pi.k , probs0 )
#----- specification using tamaan
tammodel <- "
ANALYSIS:
TYPE=LOCLCA;
NCLASSES(4)
LAVAAN MODEL:
F =~ I1__I15
I8 | 0*t1
"
mod1b <- tamaan( tammodel , resp=dat )
summary(mod1b)
#############################################################################
# SIMULATED EXAMPLE 6: DINA model with two skills
# See slca Simulated Example 5 in the CDM package
#############################################################################
#*** simulate data
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, text=
"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( runif(I, 0 , .4 ) , 2 )
slip <- round( 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
#*** Model 1: Estimate DINA model
# define E matrix which contains the anti-slipping parameters
maxK <- 2
Ngam <- I
E <- array( 0 , dim=c(I , maxK , TP , Ngam ) )
dimnames(E)[[1]] <- colnames(dat)
dimnames(E)[[2]] <- paste0("Cat" , 1:(maxK) )
dimnames(E)[[3]] <- c("S00","S10","S01","S11")
dimnames(E)[[4]] <- paste0( "antislip" , 1:I )
# define anti-slipping parameters in E
for (ii in 1:I){
# define latent responses
latresp <- 1*( alpha0 %*% Q[ii,] == sum(Q[ii,]) )[,1]
# model slipping parameters
E[ii , 2 , latresp == 1, ii ] <- 1
}
# skill space definition
theta.k <- diag(TP)
gammaslope <- rep( qlogis( .8 ) , I )
# estimate model
mod1 <- tam.mml.3pl(resp= dat , E=E , control= list(maxiter=100) , skillspace="discrete" ,
theta.k=theta.k , gammaslope=gammaslope)
summary(mod1)
# compare estimated and simulated latent class proportions
cbind( mod1$pi.k , probs0 )
# compare estimated and simulated guessing parameters
cbind( mod1$rprobs[,2,1] , guess )
# compare estimated and simulated slipping parameters
cbind( 1 - mod1$rprobs[,2,4] , slip )
#############################################################################
# SIMULATED EXAMPLE 7: Mixed Rasch model with two classes
# See slca Simulated Example 3 in the CDM package
#############################################################################
#*** simulate data
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 )
b2 <- b2 - mean(b2)
N <- 3000 # number of persons
wgt <- .25 # class probability for class 1
# class 1
dat1 <- sirt::sim.raschtype( rnorm( wgt*N ) , - b1 )
# class 2
dat2 <- sirt::sim.raschtype( rnorm( (1-wgt)*N , mean= 1 , sd=1.7) , - b2 )
dat <- rbind( dat1 , dat2 )
# The idea is that each grid point class x theta is defined as new
# dimension. If we approximate the trait distribution by 7 theta points
# and are interested in estimating 2 latent classes, then we need
# 7*2=14 dimensions.
#*** Model 1: Rasch model
# theta grid
theta.k1 <- seq( -5 , 5 , len=7 )
TT <- length(theta.k1)
#-- define theta design matrix
theta.k <- diag(TT)
#-- delta designmatrix
delta.designmatrix <- matrix( 0 , TT , ncol=3 )
delta.designmatrix[ , 1] <- 1
delta.designmatrix[, 2:3] <- cbind( theta.k1 , theta.k1^2 )
#-- define loading matrix E
E <- array( 0 , dim=c(I,2,TT,I + 1) ) # last parameter is constant 1
for (ii in 1:I){
E[ ii , 2 , 1:TT , ii ] <- -1 # '-b' in '1*theta - b'
E[ ii , 2 , 1:TT , I+1] <- theta.k1 # '1*theta' in '1*theta - b'
}
# initial gammaslope parameters
par1 <- qlogis( colMeans( dat ) )
gammaslope <- c( par1 , 1 )
# sum constraint of zero on item difficulties
gammaslope.constr.V <- matrix( 0 , I+1 , 1 )
gammaslope.constr.V[ 1:I , 1] <- 1 # Class 1
gammaslope.constr.c <- c(0)
# fixed gammaslope parameter
gammaslope.fixed <- cbind( I+1 , 1 )
# estimate model
mod1 <- tam.mml.3pl(resp= dat1 , E=E , control= list(maxiter=10) , skillspace="discrete" ,
theta.k=theta.k , delta.designmatrix=delta.designmatrix ,
gammaslope=gammaslope , gammaslope.constr.V=gammaslope.constr.V ,
gammaslope.constr.c=gammaslope.constr.c , gammaslope.fixed=gammaslope.fixed ,
notA=TRUE , est.variance=FALSE)
summary(mod1)
#*** Model 2: Mixed Rasch model with two latent classes
# theta grid
theta.k1 <- seq( -4 , 4 , len=7 )
TT <- length(theta.k1)
#-- define theta design matrix
theta.k <- diag(2*TT) # 2*7=14 classes
#-- delta designmatrix
delta.designmatrix <- matrix( 0 , 2*TT , ncol=6 )
# Class 1
delta.designmatrix[1:TT , 1] <- 1
delta.designmatrix[1:TT , 2:3] <- cbind( theta.k1 , theta.k1^2 )
# Class 2
delta.designmatrix[TT+1:TT , 4] <- 1
delta.designmatrix[TT+1:TT , 5:6] <- cbind( theta.k1 , theta.k1^2 )
#-- define loading matrix E
E <- array( 0 , dim=c(I,2,2*TT,2*I + 1) ) # last parameter is constant 1
dimnames(E)[[1]] <- colnames(dat)
dimnames(E)[[2]] <- c("Cat0","Cat1")
dimnames(E)[[3]] <- c( paste0("Class1_theta" , 1:TT) , paste0("Class2_theta" , 1:TT) )
dimnames(E)[[4]] <- c( paste0("b_Class1_" , colnames(dat)) ,
paste0("b_Class2_" , colnames(dat)) , "One")
for (ii in 1:I){
# Class 1 item parameters
E[ ii , 2 , 1:TT , ii ] <- -1 # '-b' in '1*theta - b'
E[ ii , 2 , 1:TT , 2*I+1] <- theta.k1 # '1*theta' in '1*theta - b'
# Class 2 item parameters
E[ ii , 2 , TT + 1:TT , I + ii ] <- -1
E[ ii , 2 , TT + 1:TT , 2*I+1] <- theta.k1
}
# initial gammaslope parameters
par1 <- qlogis( colMeans( dat ) )
gammaslope <- c( par1 , par1 + runif(I, -2 ,2 ) , 1 )
# sum constraint of zero on item difficulties within a class
gammaslope.center.index <- c( rep( 1, I ) , rep(2,I) , 0 )
gammaslope.center.value <- c(0,0)
# estimate model
mod1 <- tam.mml.3pl(resp= dat , E=E , control= list(maxiter=30) , skillspace="discrete" ,
theta.k=theta.k , delta.designmatrix=delta.designmatrix ,
gammaslope=gammaslope , gammaslope.center.index=gammaslope.center.index ,
gammaslope.center.value=gammaslope.center.value, gammaslope.fixed=gammaslope.fixed ,
notA=TRUE)
summary(mod1)
# latent class proportions
aggregate( mod1$pi.k , list( rep(1:2, each=TT)) , sum )
# compare simulated and estimated item parameters
cbind( b1 , b2 , - mod1$gammaslope[1:I] , - mod1$gammaslope[I + 1:I ] )
#--- specification in tamaan
tammodel <- "
ANALYSIS:
TYPE=MIXTURE;
NCLASSES(2)
NSTARTS(5,30)
LAVAAN MODEL:
F =~ I0001__I0015
ITEM TYPE:
ALL(Rasch);
"
mod1b <- tamaan( tammodel , resp=dat )
summary(mod1b)
#############################################################################
# SIMULATED EXAMPLE 8: 2PL mixture distribution model
#############################################################################
#*** simulate data
set.seed(9187)
library(sirt)
# simulate two latent classes of Rasch populations
I <- 20
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 , 16 , 18) ] <- c(1 , -.5 , -.5 , .33 , .33 , -.66 , -1 , .3)
# b2 <- scale( b2 , scale=FALSE)
b2 <- b2 - mean(b2)
N <- 4000 # number of persons
wgt <- .75 # class probability for class 1
# item slopes
a1 <- rep( 1 , I ) # first class
a2 <- rep( c(.5,1.5) , I/2 )
# class 1
dat1 <- sirt::sim.raschtype( rnorm( wgt*N ) , - b1 , fixed.a=a1)
# class 2
dat2 <- sirt::sim.raschtype( rnorm( (1-wgt)*N , mean= 1 , sd=1.4) , - b2 , fixed.a=a2)
dat <- rbind( dat1 , dat2 )
#*** Model 1: Mixed 2PL model with two latent classes
theta.k1 <- seq( -4 , 4 , len=7 )
TT <- length(theta.k1)
#-- define theta design matrix
theta.k <- diag(2*TT) # 2*7=14 classes
#-- delta designmatrix
delta.designmatrix <- matrix( 0 , 2*TT , ncol=6 )
# Class 1
delta.designmatrix[1:TT , 1] <- 1
delta.designmatrix[1:TT , 2:3] <- cbind( theta.k1 , theta.k1^2 )
# Class 2
delta.designmatrix[TT+1:TT , 4] <- 1
delta.designmatrix[TT+1:TT , 5:6] <- cbind( theta.k1 , theta.k1^2 )
#-- define loading matrix E
E <- array( 0 , dim=c(I,2,2*TT,4*I ) )
dimnames(E)[[1]] <- colnames(dat)
dimnames(E)[[2]] <- c("Cat0","Cat1")
dimnames(E)[[3]] <- c( paste0("Class1_theta" , 1:TT) , paste0("Class2_theta" , 1:TT) )
dimnames(E)[[4]] <- c( paste0("b_Class1_" , colnames(dat)) , paste0("a_Class1_" , colnames(dat)) ,
paste0("b_Class2_" , colnames(dat)) , paste0("a_Class2_" , colnames(dat)) )
for (ii in 1:I){
# Class 1 item parameters
E[ ii , 2 , 1:TT , ii ] <- -1 # '-b' in 'a*theta - b'
E[ ii , 2 , 1:TT , I + ii] <- theta.k1 # 'a*theta' in 'a*theta - b'
# Class 2 item parameters
E[ ii , 2 , TT + 1:TT , 2*I + ii ] <- -1
E[ ii , 2 , TT + 1:TT , 3*I + ii ] <- theta.k1
}
# initial gammaslope parameters
par1 <- scale( - qlogis( colMeans( dat ) ) , scale =FALSE )
gammaslope <- c( par1 , rep(1,I) , scale( par1 + runif(I, - 1.4 , 1.4 ),
scale=FALSE) , runif( I,.6,1.4) )
# constraint matrix
gammaslope.constr.V <- matrix( 0 , 4*I , 4 )
# sum of item intercepts equals zero
gammaslope.constr.V[ 1:I , 1] <- 1 # Class 1 (b)
gammaslope.constr.V[ 2*I + 1:I , 2] <- 1 # Class 2 (b)
# sum of item slopes equals number of items -> mean slope of 1
gammaslope.constr.V[ I + 1:I , 3] <- 1 # Class 1 (a)
gammaslope.constr.V[ 3*I + 1:I , 4] <- 1 # Class 2 (a)
gammaslope.constr.c <- c(0,0,I,I)
# estimate model
mod1 <- tam.mml.3pl(resp= dat , E=E , control= list(maxiter=80) , skillspace="discrete" ,
theta.k=theta.k , delta.designmatrix=delta.designmatrix ,
gammaslope=gammaslope , gammaslope.constr.V=gammaslope.constr.V ,
gammaslope.constr.c=gammaslope.constr.c , gammaslope.fixed=gammaslope.fixed ,
notA=TRUE)
# estimated item parameters
mod1$gammaslope
# summary
summary(mod1)
# latent class proportions
round( aggregate( mod1$pi.k , list( rep(1:2, each=TT)) , sum ) , 3 )
# compare simulated and estimated item intercepts
int <- cbind( b1*a1 , b2 * a2 , - mod1$gammaslope[1:I] , - mod1$gammaslope[2*I + 1:I ] )
round( int , 3 )
# simulated and estimated item slopes
slo <- cbind( a1 , a2 , mod1$gammaslope[I+1:I] , mod1$gammaslope[3*I + 1:I ] )
round(slo,3)
#--- specification in tamaan
tammodel <- "
ANALYSIS:
TYPE=MIXTURE;
NCLASSES(2)
NSTARTS(10,25)
LAVAAN MODEL:
F =~ I0001__I0020
"
mod1t <- tamaan( tammodel , resp=dat )
summary(mod1t)
#############################################################################
# EXAMPLE 9: Toy example: Exact representation of an item by a factor
#############################################################################
data(data.gpcm)
dat <- data.gpcm[,1,drop=FALSE ] # choose first item
# some descriptives
( t1 <- table(dat) )
# The idea is that we define an IRT model with one latent variable
# which extactly corresponds to the manifest item.
I <- 1 # 1 item
K <- 4 # 4 categories
TP <- 4 # 4 discrete theta points
# define skill space
theta.k <- diag(TP)
# define loading matrix E
E <- array( -99 , dim=c(I,K,TP,1 ) )
for (vv in 1:K){
E[ 1, vv , vv , 1 ] <- 9
}
# estimate model
mod1 <- tam.mml.3pl(resp= dat , E=E , skillspace="discrete" ,
theta.k=theta.k , notA=TRUE)
summary(mod1)
# -> the latent distribution corresponds to the manifest distribution, because ...
round( mod1$pi.k , 3 )
round( t1 / sum(t1) , 3 )Run the code above in your browser using DataLab