#############################################################################
# SIMULATED EXAMPLE 1: Two-dimensional IRT model with 10 items
#############################################################################
#**** data simulation
set.seed(9776)
N <- 3400 # sample size
# define difficulties
f0 <- c( .5 , .25 , -.25 , -.5 , 0 , -.5 , -.25 , .25 , .5 , 0 )
I <- length(f0)
# define loadings
f1 <- matrix( 0 , I , 2 )
f1[ 1:5,1] <- c(.8,.7,.6,.5 , .5)
f1[ 6:10 ,2] <- c(.8,.7,.6,.5 , .5 )
# covariance matrix
Pval <- matrix( c(1,.5,.5,1) , 2 , 2 )
# simulate theta
library(MASS)
theta <- MASS::mvrnorm(N , mu=c(0,0) , Sigma = Pval )
# simulate item responses
dat <- matrix( NA , N , I )
for (ii in 1:I){ # ii <- 1
dat[,ii] <- 1*(pnorm(f0[ii]+theta[,1]*f1[ii,1]+theta[,2]*f1[ii,2])>runif(N))
}
colnames(dat) <- paste0("I" , 1:I)
#**** Model 1: Two-dimensional CFA with estimated item loadings
# define pattern matrices
Pval <- .3+0*Pval
Ppatt <- 1*(Pval>0)
diag(Ppatt) <- 0
diag(Pval) <- 1
Fval <- .7 * ( f1>0)
Fpatt <- 1 * ( Fval > 0 )
# estimate model
mod1 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod1)
# EAP ability estimates
pmod1 <- R2noharm.EAP(mod1 , theta.k = seq(-4,4,len=10) )
# model fit
summary( modelfit.sirt( mod1 ))
# estimate model based on tetrachoric correlations
mod1b <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
modesttype=2)
summary(mod1b)
#**** Model 2: Two-dimensional CFA with correlated residuals
# define pattern matrix for residual correlation
Psipatt <- 0*diag(I)
Psipatt[1,2] <- 1
Psival <- 0*Psipatt
# estimate model
mod2 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
Psival=Psival , Psipatt=Psipatt )
summary(mod2)
#**** Model 3: Two-dimensional Rasch model
# pattern matrices
Fval <- matrix(0,10,2)
Fval[1:5,1] <- Fval[6:10,2] <- 1
Fpatt <- 0*Fval
Ppatt <- Pval <- matrix(1,2,2)
Pval[1,2] <- Pval[2,1] <- 0
# estimate model
mod3 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod3)
# model fit
summary( modelfit.sirt( mod3 ))
#** compare fit with NOHARM
noharm.path <- "c:/NOHARM"
P.pattern <- Ppatt ; P.init <- Pval
F.pattern <- Fpatt ; F.init <- Fval
mod3b <- R2noharm( dat = dat , model.type="CFA" ,
F.pattern = F.pattern , F.init = F.init , P.pattern = P.pattern ,
P.init = P.init , writename = "example_sim_2dim_rasch" ,
noharm.path = noharm.path , dec = "," )
summary(mod3b)
#############################################################################
# EXAMPLE 2: data.read
#############################################################################
data(data.read)
dat <- data.read
I <- ncol(dat)
#**** Model 1: Unidimensional Rasch model
Fpatt <- matrix( 0 , I , 1 )
Fval <- 1 + 0*Fpatt
Ppatt <- Pval <- matrix(1,1,1)
# estimate model
mod1 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod1)
plot(mod1) # semPaths plot
#**** Model 2: Rasch model in which item pairs within a testlet are excluded
wgtm <- matrix( 1 , I , I )
wgtm[1:4,1:4] <- wgtm[5:8,5:8] <- wgtm[ 9:12, 9:12] <- 0
# estimation
mod2 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval , wgtm=wgtm)
summary(mod2)
#**** Model 3: Rasch model with correlated residuals
Psipatt <- Psival <- 0*diag(I)
Psipatt[1:4,1:4] <- Psipatt[5:8,5:8] <- Psipatt[ 9:12, 9:12] <- 1
diag(Psipatt) <- 0
Psival <- .6*(Psipatt>0)
# estimation
mod3 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
Psival=Psival , Psipatt=Psipatt )
summary(mod3)
# allow only positive residual correlations
mod3b <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
Psival=Psival , Psipatt=Psipatt , pos.residcorr=TRUE)
summary(mod3b)
plot(mod3 , fade=FALSE)
#**** Model 4: Rasch testlet model
Fval <- Fpatt <- matrix( 0 , I , 4 )
Fval[,1] <- Fval[1:4,2] <- Fval[5:8,3] <- Fval[9:12,4 ] <- 1
Ppatt <- Pval <- diag(4)
colnames(Ppatt) <- c("g" , "A", "B","C")
Pval <- .5*Pval
# estimation
mod4 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod4)
# allow only positive variance entries
mod4b <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
pos.variance=TRUE )
summary(mod4b)
#**** Model 5: Bifactor model
Fval <- matrix( 0 , I , 4 )
Fval[,1] <- Fval[1:4,2] <- Fval[5:8,3] <- Fval[9:12,4 ] <- .6
Fpatt <- 1 * ( Fval > 0 )
Pval <- diag(4)
Ppatt <- 0*Pval
colnames(Ppatt) <- c("g" , "A", "B","C")
# estimation
mod5 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod5)
# allow only positive loadings
mod5b <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval ,
pos.loading=TRUE )
summary(mod5b)
summary( modelfit.sirt(mod5b))
plot(x= mod5b ,what="stand" , ThreshAtSide=TRUE , fade=FALSE , residuals= TRUE )
#**** Model 6: 3-dimensional Rasch model
Fval <- matrix( 0 , I , 3 )
Fval[1:4,1] <- Fval[5:8,2] <- Fval[9:12,3 ] <- 1
Fpatt <- 0*Fval
Pval <- .6*diag(3)
diag(Pval) <- 1
Ppatt <- 1+0*Pval
colnames(Ppatt) <- c("A", "B","C")
# estimation
mod6 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod6)
summary( modelfit.sirt(mod6) ) # model fit
plot(x= mod6 ,what="stand" , edge.label.cex=.75 , ThreshAtSide=TRUE ,
fixedStyle = c("black",1) , freeStyle=c("black",1) , residuals= TRUE )
#**** Model 7: 3-dimensional 2PL model
Fval <- matrix( 0 , I , 3 )
Fval[1:4,1] <- Fval[5:8,2] <- Fval[9:12,3 ] <- 1
Fpatt <- Fval
Pval <- .6*diag(3)
diag(Pval) <- 1
Ppatt <- 1+0*Pval
diag(Ppatt) <- 0
colnames(Ppatt) <- c("A", "B","C")
# estimation
mod7 <- noharm.sirt( dat=dat , Ppatt=Ppatt,Fpatt=Fpatt , Fval=Fval , Pval=Pval )
summary(mod7)
summary( modelfit.sirt(mod7) )
plot(x= mod7 ,what="est" , edge.label.cex=.75 , ThreshAtSide=TRUE , fade=FALSE ,
residuals= TRUE )
#**** Model 8: Exploratory factor analysis with 3 dimensions
# estimation
mod8 <- noharm.sirt( dat=dat , dimensions=3 )
summary(mod8)
plot(mod8)
# plot only standardized loadings larger than .40
plot(mod8 , efa.load.min=.4 )Run the code above in your browser using DataLab