# Estimate Rasch Model. J.-P. Fox, University of Twente, Augustus 2010. ## K:: number of items ## N:: number persons ## DrawZ. Data Augmentation Z: Normal Ogive Model ## DrawTheta. Draw values latent variable theta ## DrawBeta. Draw values difficulty parameter ## IRT. Main Program. DrawZ <- function(beta0,theta0,Y){ N <- nrow(Y) K <- ncol(Y) eta <- matrix(theta0,ncol=K,nrow=N) - t(matrix(beta0, ncol = N, nrow = K)) BB <- matrix(pnorm(-eta),ncol = K, nrow = N) u <- matrix(runif(N*K), ncol = K, nrow = N) tt <- matrix( ( BB*(1-Y) + (1-BB)*Y )*u + BB*Y, ncol = K, nrow = N) Z <- matrix(qnorm(tt),ncol = K, nrow = N) + matrix(eta,ncol = K, nrow = N) return(Z) } DrawTheta <- function(beta0,Z) { #prior theta N(0,1) N <- nrow(Z) K <- ncol(Z) pvar <- K + 1 thetahat <- ((Z + t(matrix(beta0,ncol= N, nrow = K))) %*% matrix(1,ncol = 1,nrow = K)) mu <- thetahat/pvar theta <- matrix(rnorm(N,mean=mu,sd=sqrt(1/pvar)),ncol=1,nrow=N) return(theta) } DrawBeta <- function(theta0,Z) { # diffuse prior N <- nrow(Z) K <- ncol(Z) Zp <- matrix(theta0,ncol=K,nrow=N) - matrix(Z,ncol=K,nrow=N) mu <- apply(Zp,2,sum)/N beta <- matrix(rnorm(K,mean=mu,sd=sqrt(1/N)),ncol=1,nrow=K) return(beta) } IRT <- function(Y,XG) { ## Initialise N <- nrow(Y) K <- ncol(Y) beta0 <- matrix(0,ncol=1,nrow=1) theta0 <- matrix(rnorm(N),ncol=1,nrow=N) EAPtheta <- matrix(0,ncol=1,nrow=N) Mbeta <- matrix(0,ncol=K,nrow=XG) for(ii in 1:XG){ Z <-DrawZ(beta0,theta0,Y) theta0<-DrawTheta(beta0,Z) beta0 <-DrawBeta(theta0,Z) EAPtheta <- ((ii - 1)*EAPtheta + theta0)/ii Mbeta[ii,1:K] <- beta0[1:K] } return(list(Mbeta,EAPtheta)) } ## Generate data N= 1000 K= 10 sim <- simRasch(N,K) Y <- matrix(sim[[1]],ncol=K,nrow=N) theta <- matrix(sim[[2]],ncol=1,nrow=N) beta <- matrix(sim[[3]],ncol=1,nrow=K) out <- IRT(Y,XG=1000)## run program apply(out[[1]],2,mean) ## estimate item difficulties plot(sim[[2]],out[[2]]) ## plot simimulated theta and estimated theta ### Use estMLIRT program S <- c(0,0,0,0) out1 <- estMLIRT(Y, S, nll=N, XG=1000, rasch=1) MLIRTout(500,out1)