# NOT RUN {
#############################################################################
## EXAMPLE 1: Unidimensional item response functions
#############################################################################
data(data.read)
dat <- data.read
#------ Definition of item response functions
#*** IRF 2PL
P_2PL <- function( par, Theta, ncat){
a <- par[1]
b <- par[2]
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,1] <- 1
for (cc in 2:ncat){
P[,cc] <- exp( (cc-1) * a * Theta[,1] - b )
}
P <- P / rowSums(P)
return(P)
}
#*** IRF 1PL
P_1PL <- function( par, Theta, ncat){
b <- par[1]
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,1] <- 1
for (cc in 2:ncat){
P[,cc] <- exp( (cc-1) * Theta[,1] - b )
}
P <- P / rowSums(P)
return(P)
}
#** created item classes of 1PL and 2PL models
par <- c( "a"=1, "b"=0 )
# define some slightly informative prior of 2PL
item_2PL <- sirt::xxirt_createDiscItem( name="2PL", par=par, est=c(TRUE,TRUE),
P=P_2PL, prior=c(a="dlnorm"), prior_par1=c( a=0 ),
prior_par2=c(a=5) )
item_1PL <- sirt::xxirt_createDiscItem( name="1PL", par=par[2], est=c(TRUE),
P=P_1PL )
customItems <- list( item_1PL, item_2PL )
#---- definition theta distribution
#** theta grid
Theta <- matrix( seq(-6,6,length=21), ncol=1 )
#** theta distribution
P_Theta1 <- function( par, Theta, G){
mu <- par[1]
sigma <- max( par[2], .01 )
TP <- nrow(Theta)
pi_Theta <- matrix( 0, nrow=TP, ncol=G)
pi1 <- dnorm( Theta[,1], mean=mu, sd=sigma )
pi1 <- pi1 / sum(pi1)
pi_Theta[,1] <- pi1
return(pi_Theta)
}
#** create distribution class
par_Theta <- c( "mu"=0, "sigma"=1 )
customTheta <- sirt::xxirt_createThetaDistribution( par=par_Theta, est=c(FALSE,TRUE),
P=P_Theta1 )
#****************************************************************************
#******* Model 1: Rasch model
#-- create parameter table
itemtype <- rep( "1PL", 12 )
partable <- sirt::xxirt_createParTable( dat, itemtype=itemtype,
customItems=customItems )
# estimate model
mod1 <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable,
customItems=customItems, customTheta=customTheta)
summary(mod1)
# estimate Rasch model by providing starting values
partable1 <- sirt::xxirt_modifyParTable( partable, parname="b",
value=- stats::qlogis( colMeans(dat) ) )
# estimate model again
mod1b <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable1,
customItems=customItems, customTheta=customTheta )
summary(mod1b)
# extract coefficients, covariance matrix and standard errors
coef(mod1b)
vcov(mod1b)
IRT.se(mod1b)
#****************************************************************************
#******* Model 2: 2PL Model with three groups of item discriminations
#-- create parameter table
itemtype <- rep( "2PL", 12 )
partable <- sirt::xxirt_createParTable( dat, itemtype=itemtype, customItems=customItems)
# modify parameter table: set constraints for item groups A, B and C
partable1 <- sirt::xxirt_modifyParTable(partable, item=paste0("A",1:4),
parname="a", parindex=111)
partable1 <- sirt::xxirt_modifyParTable(partable1, item=paste0("B",1:4),
parname="a", parindex=112)
partable1 <- sirt::xxirt_modifyParTable(partable1, item=paste0("C",1:4),
parname="a", parindex=113)
# delete prior distributions
partable1 <- sirt::xxirt_modifyParTable(partable1, parname="a", prior=NA)
#-- fix sigma to 1
customTheta1 <- customTheta
customTheta1$est <- c("mu"=FALSE,"sigma"=FALSE )
# estimate model
mod2 <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable1,
customItems=customItems, customTheta=customTheta1 )
summary(mod2)
#****************************************************************************
#******* Model 3: Cloglog link function
#*** IRF cloglog
P_1N <- function( par, Theta, ncat){
b <- par
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,2] <- 1 - exp( - exp( Theta - b ) )
P[,1] <- 1 - P[,2]
return(P)
}
par <- c("b"=0)
item_1N <- sirt::xxirt_createDiscItem( name="1N", par=par, est=c(TRUE),
P=P_1N )
customItems <- list( item_1N )
itemtype <- rep( "1N", I )
partable <- sirt::xxirt_createParTable( dat[,items], itemtype=itemtype,
customItems=customItems )
partable <- sirt::xxirt_modifyParTable( partable=partable, parname="b",
value=- stats::qnorm( colMeans(dat[,items] )) )
#*** estimate model
mod3 <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable, customItems=customItems,
customTheta=customTheta )
summary(mod3)
IRT.compareModels(mod1,mod3)
#****************************************************************************
#******* Model 4: Latent class model
K <- 3 # number of classes
Theta <- diag(K)
#*** Theta distribution
P_Theta1 <- function( par, Theta, G ){
logitprobs <- par[1:(K-1)]
l1 <- exp( c( logitprobs, 0 ) )
probs <- matrix( l1/sum(l1), ncol=1)
return(probs)
}
par_Theta <- stats::qlogis( rep( 1/K, K-1 ) )
names(par_Theta) <- paste0("pi",1:(K-1) )
customTheta <- sirt::xxirt_createThetaDistribution( par=par_Theta,
est=rep(TRUE,K-1), P=P_Theta1)
#*** IRF latent class
P_lc <- function( par, Theta, ncat){
b <- par
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,1] <- 1
for (cc in 2:ncat){
P[,cc] <- exp( Theta %*% b )
}
P <- P / rowSums(P)
return(P)
}
par <- seq( -1.5, 1.5, length=K )
names(par) <- paste0("b",1:K)
item_lc <- sirt::xxirt_createDiscItem( name="LC", par=par,
est=rep(TRUE,K), P=P_lc )
customItems <- list( item_lc )
# create parameter table
itemtype <- rep( "LC", 12 )
partable <- sirt::xxirt_createParTable( dat, itemtype=itemtype, customItems=customItems)
partable
#*** estimate model
mod4 <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable, customItems=customItems,
customTheta=customTheta)
summary(mod4)
# class probabilities
mod4$probs_Theta
# item response functions
imod4 <- IRT.irfprob( mod5 )
round( imod4[,2,], 3 )
#****************************************************************************
#******* Model 5: Ordered latent class model
K <- 3 # number of classes
Theta <- diag(K)
Theta <- apply( Theta, 1, cumsum )
#*** Theta distribution
P_Theta1 <- function( par, Theta, G ){
logitprobs <- par[1:(K-1)]
l1 <- exp( c( logitprobs, 0 ) )
probs <- matrix( l1/sum(l1), ncol=1)
return(probs)
}
par_Theta <- stats::qlogis( rep( 1/K, K-1 ) )
names(par_Theta) <- paste0("pi",1:(K-1) )
customTheta <- sirt::xxirt_createThetaDistribution( par=par_Theta,
est=rep(TRUE,K-1), P=P_Theta1 )
#*** IRF ordered latent class
P_olc <- function( par, Theta, ncat){
b <- par
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,1] <- 1
for (cc in 2:ncat){
P[,cc] <- exp( Theta %*% b )
}
P <- P / rowSums(P)
return(P)
}
par <- c( -1, rep( .5,, length=K-1 ) )
names(par) <- paste0("b",1:K)
item_olc <- sirt::xxirt_createDiscItem( name="OLC", par=par, est=rep(TRUE,K),
P=P_olc, lower=c( -Inf, 0, 0 ) )
customItems <- list( item_olc )
itemtype <- rep( "OLC", 12 )
partable <- sirt::xxirt_createParTable( dat, itemtype=itemtype, customItems=customItems)
partable
#*** estimate model
mod5 <- sirt::xxirt( dat=dat, Theta=Theta, partable=partable, customItems=customItems,
customTheta=customTheta )
summary(mod5)
# estimated item response functions
imod5 <- IRT.irfprob( mod5 )
round( imod5[,2,], 3 )
#############################################################################
## EXAMPLE 2: Multiple group models with xxirt
#############################################################################
data(data.math)
dat <- data.math$data
items <- grep( "M[A-Z]", colnames(dat), value=TRUE )
I <- length(items)
Theta <- matrix( seq(-8,8,len=31), ncol=1 )
#****************************************************************************
#******* Model 1: Rasch model, single group
#*** Theta distribution
P_Theta1 <- function( par, Theta, G ){
mu <- par[1]
sigma <- max( par[2], .01 )
p1 <- stats::dnorm( Theta[,1], mean=mu, sd=sigma)
p1 <- p1 / sum(p1)
probs <- matrix( p1, ncol=1)
return(probs)
}
par_Theta <- c(0,1)
names(par_Theta) <- c("mu","sigma")
customTheta <- sirt::xxirt_createThetaDistribution( par=par_Theta,
est=c(FALSE,TRUE), P=P_Theta1 )
customTheta
#*** IRF 1PL logit
P_1PL <- function( par, Theta, ncat){
b <- par
TP <- nrow(Theta)
P <- matrix( NA, nrow=TP, ncol=ncat)
P[,2] <- plogis( Theta - b )
P[,1] <- 1 - P[,2]
return(P)
}
par <- c("b"=0)
item_1PL <- sirt::xxirt_createDiscItem( name="1PL", par=par, est=c(TRUE), P=P_1PL)
customItems <- list( item_1PL )
itemtype <- rep( "1PL", I )
partable <- sirt::xxirt_createParTable( dat[,items], itemtype=itemtype,
customItems=customItems )
partable <- sirt::xxirt_modifyParTable( partable=partable, parname="b",
value=- stats::qlogis( colMeans(dat[,items] )) )
#*** estimate model
mod1 <- sirt::xxirt( dat=dat[,items], Theta=Theta, partable=partable,
customItems=customItems, customTheta=customTheta )
summary(mod1)
#****************************************************************************
#******* Model 2: Rasch model, multiple groups
#*** Theta distribution
P_Theta2 <- function( par, Theta, G ){
mu1 <- par[1]
mu2 <- par[2]
sigma1 <- max( par[3], .01 )
sigma2 <- max( par[4], .01 )
TP <- nrow(Theta)
probs <- matrix( NA, nrow=TP, ncol=G)
p1 <- stats::dnorm( Theta[,1], mean=mu1, sd=sigma1)
probs[,1] <- p1 / sum(p1)
p1 <- stats::dnorm( Theta[,1], mean=mu2, sd=sigma2)
probs[,2] <- p1 / sum(p1)
return(probs)
}
par_Theta <- c(0,0,1,1)
names(par_Theta) <- c("mu1","mu2","sigma1","sigma2")
customTheta2 <- sirt::xxirt_createThetaDistribution( par=par_Theta,
est=c(FALSE,TRUE,TRUE,TRUE), P=P_Theta2 )
customTheta2
#*** estimate model
mod2 <- sirt::xxirt( dat=dat[,items], group=dat$female, Theta=Theta, partable=partable,
customItems=customItems, customTheta=customTheta2, maxit=40)
summary(mod2)
IRT.compareModels(mod1, mod2)
#*** compare results with TAM package
library(TAM)
mod2b <- TAM::tam.mml( resp=dat[,items], group=dat$female )
summary(mod2b)
IRT.compareModels(mod1, mod2, mod2b)
# }
Run the code above in your browser using DataLab