#############################################################################
# EXAMPLE 1: Reading data set
#############################################################################
data(data.read)
dat <- data.read
#******
# Model 1: Rasch model with PML estimation
mod1 <- rasch.pml3( dat )
summary(mod1)
#******
# Model 2: Excluding item pairs with local dependence
# from bivariate composite likelihood
itemcluster <- rep( 1:3 , each=4)
mod2 <- rasch.pml3( dat , itemcluster = itemcluster )
summary(mod2)
#*****
# Model 3: Modelling error correlations:
# joint residual correlations for each itemcluster
error.corr <- diag(1,ncol(dat))
for ( ii in 1:3){
ind.ii <- which( itemcluster == ii )
error.corr[ ind.ii , ind.ii ] <- ii
}
# estimate the model with error correlations
mod3 <- rasch.pml3( dat , error.corr = error.corr )
summary(mod3)
#****
# Model 4: model separate residual correlations
I <- ncol(error.corr)
error.corr1 <- matrix( 1:(I*I) , ncol= I )
error.corr <- error.corr1 * ( error.corr > 0 )
# estimate the model with error correlations
mod4 <- rasch.pml3( dat , error.corr = error.corr )
summary(mod4)
#****
# Model 5: assume equal item difficulties:
# b_1 = b_7 and b_2 = b_12
# fix item difficulty of the 6th item to .1
est.b <- 1:I
est.b[7] <- 1; est.b[12] <- 2 ; est.b[6] <- 0
b.init <- rep( 0, I ) ; b.init[6] <- .1
mod5 <- rasch.pml3( dat , est.b =est.b , b.init=b.init)
summary(mod5)
#****
# Model 6: estimate three item slope groups
est.a <- rep(1:3 , each=4 )
mod6 <- rasch.pml3( dat , est.a =est.a , est.sigma=0)
summary(mod6)
#############################################################################
# EXAMPLE 2: PISA reading
#############################################################################
data(data.pisaRead)
dat <- data.pisaRead$data
# select items
dat <- dat[ , substring(colnames(dat),1,1)=="R" ]
#******
# Model 1: Rasch model with PML estimation
mod1 <- rasch.pml3( as.matrix(dat) )
## Trait SD (Logit Link) : 1.419
#******
# Model 2: Model correlations within testlets
error.corr <- diag(1,ncol(dat))
testlets <- paste( data.pisaRead$item$testlet )
itemcluster <- match( testlets , unique(testlets ) )
for ( ii in 1:(length(unique(testlets))) ){
ind.ii <- which( itemcluster == ii )
error.corr[ ind.ii , ind.ii ] <- ii
}
# estimate the model with error correlations
mod2 <- rasch.pml3( dat , error.corr = error.corr )
## Trait SD (Logit Link) : 1.384
#****
# Model 3: model separate residual correlations
I <- ncol(error.corr)
error.corr1 <- matrix( 1:(I*I) , ncol= I )
error.corr <- error.corr1 * ( error.corr > 0 )
# estimate the model with error correlations
mod3 <- rasch.pml3( dat , error.corr = error.corr )
## Trait SD (Logit Link) : 1.384
#############################################################################
# SIMULATED EXAMPLE 3: 10 locally independent items
#############################################################################
#**********
# simulate some data
set.seed(554)
N <- 500 # persons
I <- 10 # items
theta <- rnorm(N,sd=1.3 ) # trait SD of 1.3
b <- seq(-2 , 2 , length=I) # item difficulties
# simulate data from the Rasch model
dat <- sim.raschtype( theta = theta , b = b )
# estimation with rasch.pml and probit link
mod1 <- rasch.pml3( dat )
summary(mod1)
# estimation with rasch.mml2 function
mod2 <- rasch.mml3( dat )
# estimate item parameters for groups with five item parameters each
est.b <- rep( 1:(I/2) , each=2 )
mod3 <- rasch.pml3( dat , est.b=est.b )
summary(mod3)
# compare parameter estimates
summary(mod1)
summary(mod2)
summary(mod3)
#############################################################################
# SIMULATED EXAMPLE 4: 11 items and 2 item clusters with 2 and 3 items
#############################################################################
set.seed(5698)
I <- 11 # number of items
n <- 5000 # number of persons
b <- seq(-2,2, len=I) # item difficulties
theta <- rnorm( n , sd = 1 ) # person abilities
# itemcluster
itemcluster <- rep(0,I)
itemcluster[c(3,5)] <- 1
itemcluster[c(2,4,9)] <- 2
# residual correlations
rho <- c( .7 , .5 )
# simulate data (under the logit link)
dat <- sim.rasch.dep( theta , b , itemcluster , rho )
colnames(dat) <- paste("I" , seq(1,ncol(dat)) , sep="")
#***
# Model 1: estimation using the Rasch model (with probit link)
mod1 <- rasch.pml3( dat )
#***
# Model 2: estimation when pairs of locally dependent items are eliminated
mod2 <- rasch.pml3( dat , itemcluster=itemcluster)
#***
# Model 3: Positive correlations within testlets
est.corrs <- diag( 1 , I )
est.corrs[ c(3,5) , c(3,5) ] <- 2
est.corrs[ c(2,4,9) , c(2,4,9) ] <- 3
mod3 <- rasch.pml3( dat , error.corr=est.corrs )
#***
# Model 4: Negative correlations between testlets
est.corrs <- diag( 1 , I )
est.corrs[ c(3,5) , c(2,4,9) ] <- 2
est.corrs[ c(2,4,9) , c(3,5) ] <- 2
mod4 <- rasch.pml3( dat , error.corr=est.corrs )
#***
# Model 5: sum constraint of zero within and between testlets
est.corrs <- matrix( 1:(I*I) , I , I )
cluster2 <- c(2,4,9)
est.corrs[ setdiff( 1:I , c(cluster2)) , ] <- 0
est.corrs[ , setdiff( 1:I , c(cluster2)) ] <- 0
# define an error constraint matrix
itempairs0 <- mod4$itempairs
IP <- nrow(itempairs0)
err.constraint <- matrix( 0 , IP , 1 )
err.constraint[ ( itempairs0$item1 %in% cluster2 )
& ( itempairs0$item2 %in% cluster2 ) , 1 ] <- 1
# set sum of error covariances to 1.2
err.constraintV <- matrix(3*.4,1,1)
mod5 <- rasch.pml3( dat , error.corr=est.corrs ,
err.constraintM=err.constraint, err.constraintV=err.constraintV)
#****
# Model 6: Constraint on sum of all correlations
est.corrs <- matrix( 1:(I*I) , I , I )
# define an error constraint matrix
itempairs0 <- mod4$itempairs
IP <- nrow(itempairs0)
# define two side conditions
err.constraint <- matrix( 0 , IP , 2 )
err.constraintV <- matrix( 0 , 2 , 1)
# sum of all correlations is zero
err.constraint[ , 1 ] <- 1
err.constraintV[1,1] <- 0
# sum of items cluster c(1,2,3) is 0
cluster2 <- c(1,2,3)
err.constraint[ ( itempairs0$item1 %in% cluster2 )
& ( itempairs0$item2 %in% cluster2 ) , 2 ] <- 1
err.constraintV[2,1] <- 0
mod6 <- rasch.pml3( dat , error.corr=est.corrs ,
err.constraintM=err.constraint, err.constraintV=err.constraintV)
summary(mod6)
#############################################################################
# SIMULATED EXAMPLE 5: 10 Items: Cluster 1 -> Items 1,2
# Cluster 2 -> Items 3,4,5; Cluster 3 -> Items 7,8,9
#############################################################################
set.seed(7650)
I <- 10 # number of items
n <- 5000 # number of persons
b <- seq(-2,2, len=I) # item difficulties
bsamp <- b <- sample(b) # sample item difficulties
theta <- rnorm( n , sd = 1 ) # person abilities
# define itemcluster
itemcluster <- rep(0,I)
itemcluster[ 1:2 ] <- 1
itemcluster[ 3:5 ] <- 2
itemcluster[ 7:9 ] <- 3
# define residual correlations
rho <- c( .55 , .35 , .45)
# simulate data
dat <- sim.rasch.dep( theta , b , itemcluster , rho )
colnames(dat) <- paste("I" , seq(1,ncol(dat)) , sep="")
#***
# Model 1: residual correlation (equal within item clusters)
# define a matrix of integers for estimating error correlations
error.corr <- diag(1,ncol(dat))
for ( ii in 1:3){
ind.ii <- which( itemcluster == ii )
error.corr[ ind.ii , ind.ii ] <- ii
}
# estimate the model
mod1 <- rasch.pml3( dat , error.corr = error.corr )
#***
# Model 2: residual correlation (different within item clusters)
# define again a matrix of integers for estimating error correlations
error.corr <- diag(1,ncol(dat))
for ( ii in 1:3){
ind.ii <- which( itemcluster == ii )
error.corr[ ind.ii , ind.ii ] <- ii
}
I <- ncol(error.corr)
error.corr1 <- matrix( 1:(I*I) , ncol= I )
error.corr <- error.corr1 * ( error.corr > 0 )
# estimate the model
mod2 <- rasch.pml3( dat , error.corr = error.corr )
#***
# Model 3: eliminate item pairs within itemclusters for PML estimation
mod3 <- rasch.pml3( dat , itemcluster = itemcluster )
#***
# Model 4: Rasch model ignoring dependency
mod4 <- rasch.pml3( dat )
# compare different models
summary(mod1)
summary(mod2)
summary(mod3)
summary(mod4)
Run the code above in your browser using DataCamp Workspace