Los cálculos utilizados se han realizado en base a la siguiente
documentación:
- SANTE/11312/2021 https://www.eurl-pesticides.eu/docs/public/tmplt_article.asp?CntID=727
-Artículo científico https://www.mdpi.com/2077-0472/14/7/1014
- ISO
11843-1
- ISO 8466-1
- UNE-CEN/TS 17061
Aspectos a considerar:
Aunque por sencillez de calculos, en el ejemplo desarrollado se ha puesto los valores de la respuesta instrumental (area) de tres plaguicidas, lo normal es importar los datos directamente del software (sea cual sea el numero de plaguicidas)
#Function to calculate delta ISO 11843-1
doCalcDelta=function (gl){
if(gl>=2 && gl<=50){
delta=switch(gl, "2"=5.516, "3"=4.456, "4"=4.067,"5"=3.870,"6"=3.752,"7"=3.673,"8"=3.617,"9"=3.575,"10"=3.543,
"11"=3.517,"12"=3.496, "13"=3.479, "14"=3.464,"15"=3.451,"16"=3.440,"17"=3.431,"18"=3.422,"19"=3.415,
"20"=3.408,"21"=3.402,"22"=3.397,"23"=3.392, "24"=3.387,"25"=3.383,"26"=3.380,"27"=3.376,"28"=3.373,
"29"=3.370,"30"=3.367,"31"=3.365,"32"=3.362,"33"=3.360,"34"=3.358, "35"=3.356,"36"=3.354,"37"=3.352,
"38"=3.350,"39"=3.349,"40"=3.347,"41"=3.346,"42"=3.344,"43"=3.343,"44"=3.342,"45"=3.341, "46"=3.339,
"47"=3.338,"48"=3.337,"49"=3.336,"50"=3.335)
}else{
delta=2*(qt(.05, df = gl,lower.tail=FALSE))
}
return(delta)
}
#Function to evaluate the linearity ISO 8466-1
doTestLin=function(sylin,sypol,N){
DS2=((N-2)*sylin^2)-((N-3)*sypol^2)
fcalc=DS2/sypol^2
fcrit=qf(p=.01, df1=1, df2=(N-3), lower.tail=FALSE)
if(fcalc<fcrit){
calibration='Linear'
}else{
calibration='Other calibration'
}
return (calibration)
}
#Function to score
doTestLinT=function(sylin,sylinW,sypol,N){
DS2=((N-2)*sylin^2)-((N-3)*sypol^2)
fcalc=DS2/sypol^2
fcrit=qf(p=.01, df1=1, df2=(N-3), lower.tail=FALSE)
if(fcalc<fcrit){
#Comparing between linear models
Weightcalibration=doTestLin(sylin,sylinW,N)
if(Weightcalibration=='Linear'){
scorelineal=40
scorelinealW=20
scorePoly=0
}else{
scorelineal=20
scorelinealW=40
scorePoly=0
}
}else{
calibration='Polynomial'
polynomialcalibration=doTestLin(sylinW,sypol,N)
if(polynomialcalibration=='Linear'){
scorelineal=0
scorelinealW=40
scorePoly=20
}else{
scorelineal=0
scorelinealW=20
scorePoly=40
}
}
Testlist=list("Lineal" = scorelineal, "LinealW" = scorelinealW, "Polynomial" = scorePoly)
return (Testlist)
}
#Function to evaluate critical point
doPuntoCritico=function(x,b,c){
punto=(-b)/(2*c)
if(punto>0){
ptomin=min(x)
ptomax=max(x)
if(punto>ptomin && punto<ptomax){
puntocritico='Yes'
}else{
puntocritico='No'
}
}else{
puntocritico='Error (point not detected)'
}
return (puntocritico)
}
#Function to eval QC lineal
doCalcQClin=function (x,y,a,b){
mediay=mean(y)
for(i in y){
result=((y-(a+(b*x)))/mediay)^2
#interm=((i-ymodel)/mediay)^2
}
result=100*sqrt((sum(result))/(length(y)-2))
result=round(result,1)
return (result)
}
#Function to eval QC polynomical calibration
doCalcQCPol=function (x,y,a,b,c){
mediay=mean(y)
for(i in y){
result=((y-(a+(b*x)+(c*x*x)))/mediay)^2
#interm=((i-ymodel)/mediay)^2
}
result=100*sqrt((sum(result))/(length(y)-3))
result=round(result,1)
return (result)
}
#Function S-chemacal
doSchemaCal=function(all_qc,all_lcrit,sylin,sylinW,sypol,N){
min_qc=min(all_qc)
for(i in all_qc){
factor_qc=round(((min_qc*30)/all_qc),0)
}
min_lcrit=min(all_lcrit)
for(i in all_lcrit){
factor_lcrit=round(((min_lcrit*30)/all_lcrit),0)
}
factor_linealidad=doTestLinT(sylin,sylinW,sypol,N)
SchemaLin=factor_linealidad$Lineal+factor_qc[1]+factor_lcrit[1]
SchemaLinW=factor_linealidad$LinealW+factor_qc[2]+factor_lcrit[2]
SchemaPol=factor_linealidad$Polynomial+factor_qc[3]+factor_lcrit[3]
Schemalist=list("Lineal" = SchemaLin, "LinealW" = SchemaLinW, "Polynomial" = SchemaPol)
return (Schemalist)
}
#Function to evaluation slope lineal
doEvalSlopeL=function(conc,response){
DF<-data.frame(y=response, x=conc)
model=lm(y~x,data=DF)
slope=model$coefficients[2]
return (slope)
}
#Function to evaluation slope linealW
doEvalSlopeLW=function(conc,response){
DF<-data.frame(y=response, x=conc)
modelW=lm(y~x,data=DF,weights=(1/x))
slope=modelW$coefficients[2]
return (slope)
}
#Function to evaluation slope second
doEvalSlopeS=function(conc,response){
DF<-data.frame(y=response, x=conc)
modelP=lm(y~x+I(x^2),data=DF)
coefficient1=modelP$coefficients[2]
coefficient2=modelP$coefficients[3]
slope=coefficient1+(2*coefficient2*(mean(conc)))
return (slope)
}
#Function to evaluation Matrix Effect
doEvalME=function(name,conc1,response1,conc2,response2,type_cal){
if(type_cal==1){
slope_ini=doEvalSlopeL(conc1,response1)
slope_fin=doEvalSlopeL(conc2,response2)
model='Linear'
}else if (type_cal==2){
slope_ini=doEvalSlopeLW(conc1,response1)
slope_fin=doEvalSlopeLW(conc2,response2)
model='Linear-Weighted'
}else{
slope_ini=doEvalSlopeS(conc1,response1)
slope_fin=doEvalSlopeS(conc2,response2)
model='Second-order'
}
ME=100*((slope_fin/slope_ini)-1)
ME=round(ME,1)
if(ME<=20){
interp='low'
}else if(ME>20 && ME<=50){
interp='moderate'
}else{
interp='high'
}
#Results in DF
COMP = c(name)
MODEL = c(model)
MEDF = c(ME)
INTERP = c(interp)
FINALME = data.frame(COMP, MODEL, MEDF, INTERP)
return (FINALME)
}
#Function to evaluation calibration
doEvalCal=function(name,conc,response){
weights=1/conc
DF<-data.frame(y=response, x=conc)
#Linear case
model=lm(y~x,data=DF)
intercept=model$coefficients[1]
slope=model$coefficients[2]
sxy=summary(model)$sigma
r2=summary(model)$r.squared
eq=paste('Y=',round(intercept,3),'+',round(slope,3),'* X')
qcl=doCalcQClin(conc,response,intercept,slope)
coment_linear=''
#Linear-Weighted case
modelW=lm(y~x,data=DF,weights=(1/x))
interceptW=modelW$coefficients[1]
sxyW=summary(modelW)$sigma
slopeW=modelW$coefficients[2]
r2W=summary(modelW)$r.squared
eqW=paste('Y=',round(interceptW,3),'+',round(slopeW,3),'* X')
qclw=doCalcQClin(conc,response,interceptW,slopeW)
coment_linearW=''
#Second order case
modelP=lm(y~x+I(x^2),data=DF)
interceptP=modelP$coefficients[1]
coefficient1=modelP$coefficients[2]
coefficient2=modelP$coefficients[3]
eqP=paste('Y=',round(interceptP,3),'+',round(coefficient1,3),'* X +',round(coefficient2,3),'*X^2')
r2P=summary(modelP)$r.squared
sxyP=summary(modelP)$sigma
slopeP=coefficient1+(2*coefficient2*(mean(conc)))
qcP=doCalcQCPol(conc,response,interceptP,coefficient1,coefficient2)
critical_point=doPuntoCritico(conc,coefficient1,coefficient2)
if(critical_point=='Yes'){
coment_second='Danger: Critical point detected'
}else{
coment_second=''
}
#ISO 11843 linear
gl=length(conc)-2
delta=doCalcDelta(gl)
sxx=sum((conc - mean(conc))^2)
X=(sum(weights*conc))/(sum(1/conc))
lod=((delta*sxy)/slope)*sqrt(1+(1/length(conc))+(X^2/sxx))
#ISO 11843 linear weighted
sxxW=sum(weights*(conc - mean(conc))^2)
XW=(sum(weights*conc))/(sum(1/conc))
lodW=((delta*sxyW)/slopeW)*sqrt(1+(1/(sum(1/conc)))+(XW^2/sxxW))
#ISO 11843 non linear
N=length(conc)
N.hat= 1
Qxx=sum(conc^2) - sum(conc)^2/N
Qx3=sum(conc^3) -(sum(conc) * sum(conc^2)/N)
Qx4=sum(conc^4) - sum(conc^2)^2/N
x.hat=conc[1]
lodP=(sxyP*delta)/(coefficient1 + 2*coefficient2*x.hat)*sqrt(
1/N + 1/N.hat + ((x.hat - mean(conc))^2*Qx4 + (x.hat^2 - sum(conc^2)/N)^2 * Qxx -
2*(x.hat - mean(conc))*(x.hat^2 - sum(conc^2)/N)*Qx3)/
(Qx4 * Qxx - Qx3^2)
)
# Final results
all_qc=c(qcl,qclw,qcP)
all_lcrit=c(lod,lodW,lodP)
chemaresult=doSchemaCal(all_qc,all_lcrit,sxy,sxyW,sxyP,N)
#Results in DF
COMP = c(name,name,name)
MODEL = c('Linear','Linear-Weighted','Second-order')
EQUAT = c(eq,eqW,eqP)
R2 = c(round(r2,4),round(r2W,4),round(r2P,4))
XCRIT = c(signif(lod,3),signif(lodW,3),signif(lodP,3))
QC = c(qcl,qclw,qcP)
min_response=min(response)
COMENTS=c(coment_linear,coment_linearW,coment_second)
MINRESPONSE=c(min_response,min_response,min_response)
SCHEMA = c(chemaresult$Lineal,chemaresult$LinealW,chemaresult$Polynomial)
FINAL = data.frame(COMP, MODEL, R2, EQUAT, XCRIT, QC, SCHEMA, MINRESPONSE, COMENTS)
return (FINAL)
}
#Function to evaluation calibration
doGrafXcrit=function(database,type_cal,xcrit){
long_df=nrow(database)/3
if(type_cal==1){
cali='Linear'
}else if(type_cal==2){
cali='Linear-Weighted'
}else if(type_cal==3){
cali='Second-order'
}
filter=database[(database$XCRIT>xcrit & database$MODEL==cali) ,]
long_filter=nrow(filter)
long_dif=long_df-long_filter
perc_filter=round(100*(long_filter/long_df),0)
perc_resto=round(100*(long_dif/long_df),0)
porcenta=c(perc_filter,perc_resto)
tit_1=paste('Comply=', perc_filter, '%')
tit_2=paste('Non-comply=', perc_resto, '%')
library(ggplot2)
x=filter$COMP
y=round(filter$XCRIT,1)
ggplot(filter, aes(x = x, y = y)) +
geom_segment(aes(x = x, xend = x, y = 0, yend = y),
color = "gray", lwd = 1.5) +
geom_point(size =7, pch = 21, bg = 4, col = 1) +
geom_text(aes(label = y), color = "white", size = 3) +
coord_flip() +
labs(title = paste("Compounds with calibration ",cali, "and Xcritical > ",xcrit ),
subtitle = paste(tit_1," // ",tit_2),
caption = "Data source: doGrafXcrit")
}options(warn = -1)
# 1. Se definen los datos que se quiere ver la calibracion
plag=c('Endosulfan-alfa','Endosulfan-beta','Endosulfan-sulfato')
conc=c(3,5,10,50,100)
alfa=c(7505,11448,24522,115737,227536)
beta=c(4405,10404,15898,88134,164634)
sulfato=c(5073,8039,14720,70672,137259)
# 2. Crear una lista con los vectores para poder iterar sobre ellos
lista_vectores <- list(alfa, beta, sulfato)
# 3. Inicializar un vector para guardar los resultados (si se quieren tener globales)
resultados <- numeric(length(lista_vectores))
# 4. Usar un bucle 'for' para calcular la informacion de cada plaguicida
for (i in 1:length(lista_vectores)) {
# Acceder al vector actual
vector_actual <- lista_vectores[[i]]
# Calcular los datos de la calibracion
parcial=doEvalCal(plag[[i]],conc,vector_actual)
resultados[i] <- parcial
cat("\n--- Resultados Finales de: ",plag[[i]], "\n")
print(parcial)
cat("\n-------------------------------------------\n")
}##
## --- Resultados Finales de: Endosulfan-alfa
## COMP MODEL R2
## 1 Endosulfan-alfa Linear 0.9999
## 2 Endosulfan-alfa Linear-Weighted 0.9997
## 3 Endosulfan-alfa Second-order 1.0000
## EQUAT XCRIT QC SCHEMA MINRESPONSE COMENTS
## 1 Y= 1052.562 + 2270.745 * X 2.060 1.4 70 7505
## 2 Y= 628.172 + 2283.376 * X 0.782 1.6 67 7505
## 3 Y= 402.568 + 2347.847 * X + -0.768 *X^2 1.450 0.9 46 7505
##
## -------------------------------------------
##
## --- Resultados Finales de: Endosulfan-beta
## COMP MODEL R2
## 1 Endosulfan-beta Linear 0.9985
## 2 Endosulfan-beta Linear-Weighted 0.9961
## 3 Endosulfan-beta Second-order 0.9996
## EQUAT XCRIT QC SCHEMA MINRESPONSE COMENTS
## 1 Y= 1075.323 + 1655.348 * X 8.41 5.5 69 4405
## 2 Y= 216.117 + 1680.919 * X 2.71 5.9 67 4405
## 3 Y= -920.79 + 1892.126 * X + -2.359 *X^2 4.96 3.4 46 4405
##
## -------------------------------------------
##
## --- Resultados Finales de: Endosulfan-sulfato
## COMP MODEL R2
## 1 Endosulfan-sulfato Linear 0.9999
## 2 Endosulfan-sulfato Linear-Weighted 0.9998
## 3 Endosulfan-sulfato Second-order 1.0000
## EQUAT XCRIT QC SCHEMA MINRESPONSE COMENTS
## 1 Y= 1287.546 + 1365.031 * X 2.480 1.6 35 5073
## 2 Y= 1049.658 + 1372.111 * X 0.533 1.8 78 5073
## 3 Y= 742.378 + 1429.699 * X + -0.644 *X^2 0.748 0.5 51 5073
##
## -------------------------------------------