################################################################## # IRT Functions (see main # program below - click the "Run Program" # icon to submit program) irf<-function(theta=0,a=1,b=0,cc=0,type="3pl"){ if (type=="3pl"){ prob<-cc+(1-cc)/(1+exp(-1.7*a*(theta-b))) } if (type=="norm"){ prob<-cc+(1-cc)*pnorm(a*(theta-b)) } return(prob) } ################################################################### irfplot<-function(a=1,b=0,cc=0,type="3pl",label="",dashed="F"){ prob<-rep(0,71) theta<-rep(0,71) for (i in 1:71){ theta[i]<-(i-36)/10 prob[i]<-irf(theta[i],a,b,cc,type) } if (dashed=="F"){ plot(theta,prob,xlab="theta",ylab="prob",type="l", xlim=c(-3.5,3.5),ylim=c(0,1),main=label) } if (dashed=="T"){ plot(theta,prob,xlab="theta",ylab="prob",type="l", xlim=c(-3.5,3.5),ylim=c(0,1),main=label,lty=2) } } ############################################################## irtgen<-function(nexmn=10,avec=c(1),bvec=c(0),cvec=c(0),type="3pl", mthet=0,sdthet=1){ n<-length(avec) data<-matrix(0,nrow=nexmn,ncol=n) ability<-rnorm(nexmn,mean=mthet,sd=sdthet) p<-t(apply(as.array(ability),1,irf,avec,bvec,cvec,type)) try<-runif(nexmn*n,0,1) data[try < p]<-1 return(data) } ###################################################################### items<-function(testdata,correct=T,dpct=0.27,digits=4){ n<-ncol(testdata) nexmn<-nrow(testdata) x<-apply(testdata,1,sum) medx<-median(x) lowx<-sort(x)[floor((nexmn+1)*dpct)] highx<-sort(x)[ceiling((nexmn+1)*(1-dpct))] pvalue<-as.numeric(apply(testdata,2,mean)) var<-as.numeric(apply(testdata,2,var)*(nexmn-1)/nexmn) D<-rep(0,n) PBis<-rep(0,n) Bis<-rep(0,n) for (i in 1:n){ if (correct==T){ xx<-x-testdata[,i] lowx<-sort(xx)[floor((nexmn+1)*dpct)] highx<-sort(xx)[ceiling((nexmn+1)*(1-dpct))] } else {xx<-x} pu<-mean(testdata[xx >= highx,i]) pl<-mean(testdata[xx <= lowx,i]) D[i]<-pu-pl PBis[i]<-cor(xx,testdata[,i]) Bis[i]<-sqrt(pvalue[i]*(1-pvalue[i]))*PBis[i]/dnorm(qnorm(pvalue[i])) } return(round((10^digits)*cbind(pvalue,var,D,PBis,Bis))/(10^digits))} ######################################################################## irfinfo<-function(theta=0,a=1,b=0,cc=0,type="3pl"){ n<-length(a) info<-0 if (type=="3pl"){ for (i in 1:n){ info<-info+2.89*a[i]^2*(1-cc[i])/( (cc[i]+exp(1.7*a[i]*(theta-b[i])))* ((1+exp(-1.7*a[i]*(theta-b[i])))^2) ) } } else{ for (i in 1:n){ info<-info+ a[i]^2*(1-cc[i])^2*dnorm(a[i]*(theta-b[i]))^2/ (irf(theta=theta,a=a[i],b=b[i],cc=cc[i],type="norm")* (1-irf(theta=theta,a=a[i],b=b[i],cc=cc[i],type="norm"))) } } return(info) } ############################################################################ irfinfoplot<-function(a=1,b=0,cc=0,type="3pl",label="",dashed="F",maxy=0){ info<-rep(0,71) theta<-rep(0,71) for (i in 1:71){ theta[i]<-(i-36)/10 info[i]<-irfinfo(theta=theta[i],a=a,b=b,cc=cc,type=type) } if (maxy==0){maxy<-1.1*max(info)} if (dashed=="F"){ plot(theta,info,xlab="theta",ylab="information",type="l", xlim=c(-3.5,3.5),ylim=c(0,maxy),main=label) } if (dashed=="T"){ plot(theta,info,xlab="theta",ylab="information",type="l", xlim=c(-3.5,3.5),ylim=c(0,maxy),main=label,lty=2) } } ######################################################################## stratify<-function(respmat,it1,it2,include=F){ nex<-nrow(respmat) nit<-ncol(respmat) outarray<-array(0,dim=c(2,2,nit+1)) numatscore<-rep(0,nit+1) if (include==T){ score<-apply(respmat,1,sum)+1 } else{ score<-apply(respmat,1,sum)-respmat[,it1]-respmat[,it2]+3 } for (j in 1:nex){ outarray[respmat[j,it1]+1,respmat[j,it2]+1,score[j]]<- outarray[respmat[j,it1]+1,respmat[j,it2]+1,score[j]]+1 numatscore[score[j]]<-numatscore[score[j]]+1 } return(outarray[,,(1:(nit+1))[(numatscore>1)]]) } ######################################################################### irtmh<-function(respmat,it1,it2,include=F){ temp<-mantelhaen.test(stratify(respmat,it1,it2,include)) zmh<-as.numeric(sqrt(temp$statistic)*sign(temp$estimate-1)) return(zmh=zmh,p.equal=2*min(pnorm(zmh),1-pnorm(zmh)) ,p.less=pnorm(zmh)) } ############################################################################ accov<-function(respmat,include=F){ nit<-ncol(respmat) nex<-nrow(respmat) ccov<-matrix(0,nrow=nit,ncol=nit) tscore<-apply(respmat,1,sum) for (i in 1:(nit-1)) { for (l in (i+1):nit) { sil<-tscore-respmat[,i]-respmat[,l]+1 irfi<-tapply(respmat[,i],sil,sum) irfl<-tapply(respmat[,l],sil,sum) irfil<-tapply(respmat[,i]*respmat[,l],sil,sum) numatscore<-table(sil) numatscore2<-numatscore[numatscore>1] irfi<-irfi[numatscore>1]/numatscore2 irfl<-irfl[numatscore>1]/numatscore2 irfil<-irfil[numatscore>1]/numatscore2 ccov[i,l]<-sum((irfil-irfi*irfl)*numatscore2)/sum(numatscore2) ccov[l,i]<-ccov[i,l] } } return(ccov) } ########################################################################## hcaccprox<-function(respmat){ library(mva) dmat<-as.dist(-1*accov(respmat)+1) plot(hclust(dmat,method="average")) } ############################################################################# difstrata<-function(refmat,focmat,it){ nref<-nrow(refmat) nfoc<-nrow(focmat) refscore<-apply(refmat,1,sum) focscore<-apply(focmat,1,sum) nit<-max(refscore,focscore) outarray<-array(0,dim=c(2,2,nit+1)) numatscore<-rep(0,nit+1) refscore<-apply(refmat,1,sum)+1 focscore<-apply(focmat,1,sum)+1 for (j in 1:nref){ outarray[refmat[j,it]+1,1,refscore[j]]<- outarray[refmat[j,it]+1,1,refscore[j]]+1 numatscore[refscore[j]]<-numatscore[refscore[j]]+1 } for (j in 1:nfoc){ outarray[focmat[j,it]+1,2,focscore[j]]<- outarray[focmat[j,it]+1,2,focscore[j]]+1 numatscore[focscore[j]]<-numatscore[focscore[j]]+1 } return(outarray[,,(1:(nit+1))[(numatscore>1)]]) } ######################################################################### difmh<-function(refmat,focmat,it){ temp<-mantelhaen.test(difstrata(refmat,focmat,it)) zmh<-as.numeric(sqrt(temp$statistic)*sign(temp$estimate-1)) return(zmh=zmh,p.equal=2*min(pnorm(zmh),1-pnorm(zmh)), D=as.numeric(-2.35*log(temp$estimate))) } ############################################################################# alpha<-function(testdata){ n<-ncol(testdata) nexmn<-nrow(testdata) x<-apply(testdata,1,sum) s2y<-diag(var(testdata))*(nexmn-1)/nexmn s2x<-var(x)*(nexmn-1)/nexmn alpha<-(n/(n-1))*(1-sum(s2y)/s2x) s2yy<-(((var(testdata)-diag(diag(var(testdata))))*(nexmn-1)/nexmn))^2 lambda2<-1-sum(s2y)/s2x+sqrt((n/(n-1))*sum(s2yy))/s2x return(alpha,lambda2)} ############################################################################## sem<-function(testdata){ nexmn<-nrow(testdata) x<-apply(testdata,1,sum) rhat<-alpha(testdata) sx<-sqrt(var(x)*(nexmn-1)/nexmn) sem1<-sx*sqrt(1-rhat$alpha) sem2<-sx*sqrt(1-rhat$lambda2) return(sem1,sem2)} ############################################################################## xandt<-function(testdata,rtype="lambda2",conf.level=0.95){ nexmn<-nrow(testdata) x<-apply(testdata,1,sum) avgx<-mean(x) varx<-var(x)*(nexmn-1)/nexmn sdx<-sqrt(varx) maxx<-max(x) xrange<-c(0:maxx) count<-rep(0,maxx) for (i in 0:maxx){ count[i+1]<-sum(rep(1,nexmn)[x==i])} x<-c(0:maxx) z<--1*qnorm((1-conf.level)/2) if (rtype=="alpha") { rhat<-alpha(testdata)$alpha sem.est<-sem(testdata)$sem1 } else { rhat<-alpha(testdata)$lambda2 sem.est<-sem(testdata)$sem2 } losem<-x-z*sem.est hisem<-x+z*sem.est treg<-rhat*x+(1-rhat)*avgx see.est<-sem.est*sqrt(rhat) loreg<-treg-z*see.est hireg<-treg+z*see.est return(avgx,varx,sdx,cbind(x,treg,count,losem,hisem,loreg,hireg))} ############################################################################# ### Main Program a<-rep(1,10) b<-rep(0,10) cc<-rep(0,10) sampdat<-irtgen(nexmn=100,avec=a,bvec=b,cvec=cc,type="norm",mthet=0,sdthet=1) sampdat testdata<-sampdat items(sampdat)