## RaschMeasurement provides functions to analyze Inspection experiments with ltm models
## This is the official library for general use.
## New functions must be compatible with InspectionRealm

library('eRm')
#library('ltm')

PACKAGE<-"eRm"

##################################################################################################################
########################## L T M   I M P L E M E N T A T I O N ###################################################
##################################################################################################################

ltmEstimateItemPar<-function(responseMatrix){
## Estimates the item parameters of a response matrix
## Returns a rasch object
	ltmfit<-rasch(data=responseMatrix)
	ltmfit
}

ltmPlotItemCurves<-function(ltmfit,icc=TRUE,iic=TRUE,tif=TRUE,se=TRUE){
## Plots ICC, IIC, TIF for estimated item parameters
## Argument: rasch object as returned by EstimateItemPar
## Remember, that three plots are drawn, e.g. par(mfrow=c(3,1))
## No return value
	if (iic){
		plot(ltmfit, legend = TRUE, cx = "bottomright", lwd = 3, cex.main = 1.5, cex.lab = 1.3, cex = 1.1)
	}
	if (icc){
		plot(ltmfit, type = "IIC", annot = FALSE, lwd = 3, cex.main = 1.5, cex.lab = 1.3)
	}
	if (tif){
		plot(ltmfit, type = "IIC", items = 0, lwd = 3, cex.main = 1.5,     cex.lab = 1.3)
	}
	if (se){
		vals <- plot(ltmfit, type = "IIC", items = 0, plot = FALSE)
		plot(vals[, "z"], 1 / sqrt(vals[, "info"]), type = "l", lwd = 2,xlab = "Ability", ylab = "Standard Error",main = "Standard Error of Measurement")

	}
}

ltmEstimatePersonPar<-function(ltmfit, responseMatrix, plot=FALSE){
## Estimates the person parameters from a given ltm and the original response matrix
## Arguments: 
##	an rasch object as returned by EstimateItemPar
##	the original response matrix
## Returns a vector with person parameters
## plots a boxplot for the distribution of person parameters
	scores<-factor.scores(ltmfit,responseMatrix)
	#regressAbil<-data.frame(estimated=scores[["score.dat"]][["z1"]],true=Person$Ability,row.names=1:noabil)
	#abil_lm<-lm(estimated~true,regressAbil)
	#abil_fitted=fitted(abil_lm)
	if (plot){
		PlotAbilityDist(scores)
	}
	scores
}

ltmPlotAbilityDist<-function(abilities){
## Plots regarding the distribution of person parameters in the sample
## Argument: fscores object as returned by EstimatePersonPar
## No return value
	boxplot(abilities[["score.dat"]][["z1"]])
}

ltmModelTest<-function(ltmfit,LR=FALSE,unidim=TRUE,plot=TRUE){
## Conduct a model test
if (LR){
	warning("Not implemented with ltm")
}
if (unidim){
	unidimresult<-unidimTest(ltmfit)
}
if (plot){
	if (unidim){
		plot(unidimresult, type = "b", pch = 1:2)
		legend("topright", c("Real Data", "Average Simulated Data"), lty = 1, pch = 1:2, col = 1:2, bty = "n")
	}
}
unidimresult
}

CorrelationTest<-function(simvector, estvector, plot=TRUE){
## Correlates two vectors: the input vector of the simulation and the estimated parameter vector
## Arguments: two vectors
## Plots a scatter plot with a linear curve by default
## Returns and prints a correlation object
}

ltmDIFTest<-function(ltmfitA,ltmfitB,plotbar=TRUE,plotpoints=TRUE){
## Performs a corrrelation test between two item parameter sets
## Two plots for graphical checking off differences
## Returns the result of the correlation test
	coefA<-coef(ltmfitA)	
	coefB<-coef(ltmfitB)
	rasch_lm<-lm(coefA[,1]~coefB[,1])
	cortest<-cor.test(coefA[,1],coefB[,1], method = c("pearson"))
	if (plotpoints){
		plot(coefA[,1],coefB[,1],xlab="Method A",ylab="Method B", main="Correlation of item parameters")
		abline(rasch_lm)
	}
	if (plotbar){
		difficulties<-t(data.frame(DifficultyA=coefA[,1],DifficultyB=coefB[,1]))
		barplot(as.matrix(difficulties),beside=TRUE)
	}
	cortest
}

##################################################################################################################
########################## e R m   I M P L E M E N T A T I O N ###################################################
##################################################################################################################

eRmEstimateItemPar<-function(responseMatrix,sum0=FALSE){
## Estimates the item parameters of a response matrix
## Returns a rasch object
	ltmfit<-RM(responseMatrix,sum0=TRUE)
	ltmfit
}

eRmModelTest<-function(ltmfit,unidim=FALSE,LR=TRUE,plot=TRUE){
## Conduct a model test
if (unidim){
	warning("Not implemented with eRm")
}
if (LR){
	lrres<-LRtest(ltmfit)
}
if (plot){
	if(LR){plotGOF(lrres)}
}
lrres
}

eRmPlotItemCurves<-function(ltmfit,icc=TRUE,iic=FALSE,tif=FALSE,se=FALSE){
## Plots ICC, IIC, TIF for estimated item parameters
## Argument: rasch object as returned by EstimateItemPar
## Remember, that three plots are drawn, e.g. par(mfrow=c(3,1))
## No return value
	if (iic){
		warning("IIC plots not supported by eRm")
	}
	if (icc){
		plotjointICC(ltmfit)
	}
	if (tif){
		warning("TIF plots not supported by eRm")
	}
	if (se){
		warning("Standard Error plots not supported by eRm")
	}
}

eRmEstimatePersonPar<-function(ltmfit, responseMatrix, plot=FALSE){
## Estimates the person parameters from a given ltm and the original response matrix
## Arguments: 
##	an rasch object as returned by EstimateItemPar
##	the original response matrix
## Returns a vector with person parameters
## plots a boxplot for the distribution of person parameters
	personpar<-person.parameter(ltmfit)
	scores<-c()
	for (group in personpar[["thetapar"]][1:length(personpar[["thetapar"]])]){scores<-c(x,group)}
	if (plot){
		PlotAbilityDist(scores)
	}
	scores
}

eRmPlotAbilityDist<-function(abilities){
## Plots regarding the distribution of person parameters in the sample
## Argument: fscores object as returned by EstimatePersonPar
## No return value
	boxplot(abilities,main=paste("Distribution of theta"))
}

if(PACKAGE=="ltm"){
	EstimateItemPar<-ltmEstimateItemPar
	PlotItemCurves<-ltmPlotItemCurves
	EstimatePersonPar<-ltmEstimatePersonPar
	PlotAbilityDist<-ltmPlotAbilityDist
	ModelTest<-ltmModelTest
	DIFTest<-ltmDIFTest
}
if(PACKAGE=="eRm"){
	EstimateItemPar<-eRmEstimateItemPar
	PlotItemCurves<-eRmPlotItemCurves
	EstimatePersonPar<-eRmEstimatePersonPar
	PlotAbilityDist<-eRmPlotAbilityDist
	ModelTest<-eRmModelTest
	#DIFTest<-eRmDIFTest
}

##################################################################################################################
########################## O W N   I M P L E M E N T A T I O N S #################################################
##################################################################################################################

TwoPLFormula<-function(theta,epsilon,alpha=1){
## implements the Birnbaum formula
	term<-exp(alpha*(theta-epsilon))
	out=term/(1+term)
}

plotICC<-function(e,a,col='Black',annot=TRUE,text=e,lty=1) {
## plot a nice ICC curve (like in the IRT inspection paper)
	curve(TwoPLFormula(x,e,a),add=TRUE,col=col,lty=lty)
	if(annot){
		lines(x=c(e,e),y=c(-0.03,0.5),col=col,lty=2)
		text(e,0.5,text,pos=4,col=col)
	}
}