irf<-function(theta=0,a=1,b=0,cc=0,type="3pl"){
n<-length(a)
prob<-0
if (type=="3pl"){
  for (i in 1:n){
    prob<-prob+cc[i]+(1-cc[i])/(1+exp(-1.7*a[i]*(theta-b[i])))
  }
}
else{
  for (i in 1:n){
  prob<-cc[i]+(1-cc)[i]*pnorm(a[i]*(theta-b[i]))
  }
}
return(prob)
}

irfplot<-function(a=1,b=0,cc=0,type="3pl",label="",dashed=F){
prob<-rep(0,71)
theta<-rep(0,71)
lim<-length(a)
if (lim==1){lab="Probability"}
else {lab="Expected Score"}
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="Ability",ylab=lab,type="l",
xlim=c(-3.5,3.5),ylim=c(0,lim),main=label)
}
if (dashed==T){
plot(theta,prob,xlab="Ability",ylab=lab,type="l",
xlim=c(-3.5,3.5),ylim=c(0,lim),main=label,lty=2)
}
}

items<-function(testdat,correct=T,dpct=0.27,digits=4){
n<-ncol(testdat)
nexmn<-nrow(testdat)
x<-apply(testdat,1,sum)
medx<-median(x)
lowx<-sort(x)[floor((nexmn+1)*dpct)]
highx<-sort(x)[ceiling((nexmn+1)*(1-dpct))]
pvalue<-as.numeric(apply(testdat,2,mean))
var<-as.numeric(apply(testdat,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-testdat[,i]
lowx<-sort(xx)[floor((nexmn+1)*dpct)]
highx<-sort(xx)[ceiling((nexmn+1)*(1-dpct))]
}
else {xx<-x}
pu<-mean(testdat[xx >= highx,i])
pl<-mean(testdat[xx <= lowx,i])
D[i]<-pu-pl
PBis[i]<-cor(xx,testdat[,i])
Bis[i]<-sqrt(pvalue[i]*(1-pvalue[i]))*PBis[i]/dnorm(qnorm(pvalue[i]))
}
round((10^digits)*cbind(pvalue,var,D,PBis,Bis))/(10^digits)}

resids<-function(dat,scores,pars,ngroup=10,qtype="axis"){
  thetas<-as.matrix(scores)[,1]
  pars<-as.matrix(pars)
  dat<-as.matrix(dat)
  avec<-pars[,1]
  bvec<-pars[,3]
  cvec<-pars[,5]
  n<-ncol(dat)
  N<-nrow(dat)
  pmat<-matrix(0,ncol=n,nrow=N)
  for (i in 1:n){
    pmat[,i]<-irf(thetas,avec[i],bvec[i],cvec[i])
  }
  rmat<-dat-pmat
  vmat<-pmat*(1-pmat)
  if (ngroup>0){
    o<-matrix(0,nrow=ngroup,ncol=n)
    r<-matrix(0,nrow=ngroup,ncol=n)
    s<-matrix(0,nrow=ngroup,ncol=n)
    t<-rep(0,ngroup)
    c<-rep(0,ngroup)
    if (qtype=="dat"){
      q<-quantile(thetas,(0:ngroup)/ngroup)}
    else q<-seq(-3,3,6/ngroup)
    for (i in 1:ngroup){
      c[i]<-sum((thetas>=q[i])&(thetas<=q[i+1]))
      if (c[i]>1){
        o[i,]<-apply(dat[(thetas>=q[i])&(thetas<=q[i+1]),],2,mean)
        r[i,]<-apply(rmat[(thetas>=q[i])&(thetas<=q[i+1]),],2,mean)
        c[i]<-sum((thetas>=q[i])&(thetas<=q[i+1]))
        s[i,]<-sqrt(apply(vmat[(thetas>=q[i])&(thetas<=q[i+1]),],2,mean)/c[i])
        t[i]<-mean(thetas[(thetas>=q[i])&(thetas<=q[i+1])])
      }
    }
    o<-o[c>1,]
    r<-r[c>1,]
    s<-s[c>1,]
    t<-t[c>1]
    c<-c[c>1]/sum(c[c>1])
  }
  if (ngroup==0){
    o<-dat
    r<-rmat
    s<-sqrt(pmat*(1-pmat))
    t<-thetas
  }
  list(obspct=o,resid=r,sresid=r/s,theta=t,weights=c)}


irfresidplot<-function(item,dat,scores,pars,ngroup=10,qtype="axis",label=""){
  temp<-resids(dat,scores,pars,ngroup,qtype)
  irfplot(a=pars[item,1],b=pars[item,3],c=pars[item,5])
  par(new=T)
  plot(temp$theta,temp$obspct[,item],xlim=c(-3.5,3.5),ylim=c(0,1),main=label,
    xlab="",ylab="")
}

itemresidplot<-function(item,dat,scores,pars,ngroup=10,qtype="axis",label="",xlim=c(-3.5,3.5),ylim=NULL){
  temp<-resids(dat,scores,pars,ngroup,qtype)
  plot(temp$theta,temp$sresid[,item],xlim=xlim,ylim=ylim,main=label,
    xlab="Ability",ylab="Standardized Residual")
  lines(xlim,c(0,0))}

bisresidplot<-function(dat,scores,pars,ngroup=10,qtype="axis",weight=T,label=""){
  tempx<-items(dat)[,5]
  if (weight==F){
    tempy<-apply(abs(resids(dat,scores,pars,ngroup,qtype)$sresid),2,mean)
  }
  if (weight==T){
    temp<-resids(dat,scores,pars,ngroup,qtype)
    wgt<-temp$weights
    tempy<-apply(abs(temp$sresid),2,weighted.mean,w=wgt)
  }
  plot(tempx,tempy,xlab="Item Biserial Correlations",ylab="Average Absolute SR",main=label)
}

pvalresidplot<-function(dat,scores,pars,ngroup=10,qtype="axis",weight=T,label="",xlim=NULL,ylim=NULL){
  tempx<-items(dat)[,1]
  if (weight==F){
    tempy<-apply(abs(resids(dat,scores,pars,ngroup,qtype)$sresid),2,mean)
  }
  if (weight==T){
    temp<-resids(dat,scores,pars,ngroup,qtype)
    wgt<-temp$weights
    tempy<-apply(abs(temp$sresid),2,weighted.mean,w=wgt)
  }
  plot(tempx,tempy,xlab="Item Percent Correct",ylab="Average Absolute SR",main=label,xlim=xlim,ylim=ylim)
}

residsummary<-function(dat,scores,pars,ngroup=10,qtype="axis",label="",nclass=NULL,...){
  temp<-c(resids(dat,scores,pars,ngroup,qtype)$sresid)
  hist(temp,main=label,nclass=nclass,xlab="Standardized Residuals",...)
  out<-matrix(0,ncol=4,nrow=1)
  temp<-abs(temp)
  n<-length(temp)
  out[1]<-sum(temp<=1)/n
  out[2]<-sum((temp>1)&(temp<=2))/n
  out[3]<-sum((temp>2)&(temp<=3))/n
  out[4]<-sum(temp>3)/n
  out<-round(100*out,1)
  colnames(out)<-list("[0,1]","(1,2]","(2,3]","(3,inf)")
  out}