R/validation.R

In AQuality: Water and Measurements Quality

validation <-
function(datavalidation,numest=NULL,measurementunit=NULL,techniquename=NULL,date=Sys.Date(),graph=FALSE){

if(numest=='NULL'){
stop('argument numest was avoided')
}
#linear regression with data from day=1
curveday1<-lm(abs~concentration,data=datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),])

rcurveday1<-cor.test(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$abs,datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration)
Rsquare<- rcurveday1$estimate

day1<-datavalidation[datavalidation$day==1,]
meanabs<-mean(day1$abs)
meanconcentration<-mean(day1$concentration)
sdabs<-sd(day1$abs)
sdconcentration<-sd(day1$concentration)

sumabs<-c()
for(i in 1:nrow(day1)){
sumabs<-c(sumabs,(day1[i,3]-meanabs)^2)
}
sumabs<-sum(sumabs)
sumconcentration<-c()
for(i in 1:nrow(day1)){
sumconcentration<-c(sumconcentration,(day1[i,2]-meanconcentration)^2)
}
sumconcentration<-sum(sumconcentration)
SE<-sqrt(sumabs/(nrow(day1)-2))/sqrt(sumconcentration)

if(graph=='TRUE'){
fcurveday1<-function(x){curveday1$coefficients[1]+curveday1$coefficients[2]*x}
plot(0,0,type="n",xlim=c(min(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration),max(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration)),ylim=c(min(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$abs),max(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$abs)))
plot(fcurveday1,min(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration),max(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration),lwd=2,col="blue")
points(datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$concentration,datavalidation[datavalidation$day==1 & (datavalidation$tube=="s" | datavalidation$tube=="b"),]$abs,pch=20,col="red",cex=1.5)
}
#sensitivity calculation
sensitivity=SE*100/curveday1$coefficients[2]
# detection limit calculation
LOD=mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="b",]$abs)+3.3*sd(datavalidation[datavalidation$day==1 & datavalidation$tube=="b",]$abs)
LOD= (LOD -curveday1$coefficients[1])/curveday1$coefficients[2] 
# cuantification limit calculation
LOQ=mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="b",]$abs)+10*sd(datavalidation[datavalidation$day==1 & datavalidation$tube=="b",]$abs)
LOQ= (LOQ -curveday1$coefficients[1])/curveday1$coefficients[2]

# veracidad calculation
V= abs(((mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$abs)-curveday1$coefficients[1])/curveday1$coefficients[2])-mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$concentration))*100/mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$concentration)
 

# repetitivity calculation (CVr%)
CVr= sd(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$abs)*100/mean(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$abs)


#intermediate precision calculation (CVR%)
CVR= sd(c(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$abs,datavalidation[datavalidation$day==2 & datavalidation$tube=="qc",]$abs,datavalidation[datavalidation$day==3 & datavalidation$tube=="qc",]$abs))*100/mean(c(datavalidation[datavalidation$day==1 & datavalidation$tube=="qc",]$abs,datavalidation[datavalidation$day==2 & datavalidation$tube=="qc",]$abs,datavalidation[datavalidation$day==3 & datavalidation$tube=="qc",]$abs))


# summary
writeLines(paste("Validation parameters technique: ", techniquename))
writeLines(paste("slope: ", round(curveday1$coefficients[2],digits=2)))
writeLines(paste("Linearity, R square = ", (trunc(Rsquare*1000))/1000))
writeLines(paste("Sensitivity ", (trunc(sensitivity*100))/100, " %")) 
writeLines(paste("LOD ", (trunc(LOD*100))/100, measurementunit))
writeLines(paste("LOQ ", (trunc(LOQ*100))/100, measurementunit))
writeLines(paste("Accuracy ", 100-(trunc(V*100))/100, "%"))
writeLines(paste("CVr% :intraassay repetitivity ", (trunc(CVr*100))/100, "%"))
writeLines(paste("CVR% :interassay repetitivity ", (trunc(CVR*100))/100, "%"))
writeLines(paste("linear range: ", (trunc(LOQ*100))/100, "-", max(day1$concentration), measurementunit))
writeLines(paste("Uncertainty: Urel(%): ", trunc(2*sqrt((CVR/sqrt(3))^2+(V/sqrt(3))^2)*100)/100, "%"))
}