# Supplementary Material
# Fiber Features, a dedicated software for the automatic
# characterization of Anthropogenic Microfibers from environmental
# samples

# Di Mauro R.[a,b], Russo R. E.[c, *] , Alves N.M.[a] , Rodriguez J.S.[a] ,
# Rodriguez J.[a] , Diaz M.V.[a,b,d] , Pastore J.I.[b,c]
# [a] Instituto Nacional de Investigación y Desarrollo Pesquero (INIDEP), Paseo Victoria
# Ocampo 1, Mar del Plata, B7602HSA , Buenos Aires, Argentina
# [b] Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Godoy Cruz
# 2290, CABA, 1425, Buenos Aires, Argentina
# [c] Laboratorio de Instrumentación Virtual y Robótica Aplicada (LIVRA-ICYTE), Juan B.
# Justo, Mar del Plata, 7600, Buenos Aires, Argentina
# [d] Instituto de Investigaciones Marinas y Costeras (IIMyC), Paseo Victoria Ocampo 1,
# Mar del Plata, B7602HSA , Buenos Aires, Argentina
# [*] Corresponding author. Email address: rodrigo.e.russo@fi.mdp.edu.ar (Russo R. E.)

rm(list = ls())
#read data
getwd() # show working directory

data_CV<-read.csv("CV.csv",header=T)
str (data_CV)

# install.packages("ggplot2", type = "binary")
library(ggplot2)

p <- ggplot(data_CV, aes(x = CV, fill = ind)) + 
  geom_density(aes(cv_hw), fill = "blue", alpha=0.3, col = rgb (0,0,1,0.2), linewidth=1.5)+
  geom_density(aes(cv_lw), fill = "red",  alpha=0.3, col = rgb (1,0,0,0.2), linewidth=1.5)+
  scale_fill_discrete(name = "CV density values", labels = c("HW", "LW")) 
p

#Stats metrics 
#CV_Highly Weathered (HW)
mean(na.omit(data_CV$cv_hw))
quantile(na.omit(data_CV$cv_hw))

#CV_Low Weathered
mean(na.omit(data_CV$cv_lw))
quantile(na.omit(data_CV$cv_lw))

#CV threshold value
#Calculate Youden threshold for the intersection of empiric distribution:

simple_roc2 <- function(ref, test){
  tab = table(test, ref) # head(tab)
  data.frame(threshold=paste('>=',rownames(tab)), 
             ref0 = tab[,1], 
             ref1 = tab[,2],  
             FPR = rev(cumsum(rev(tab[,1])/sum(tab[,1]))),
             TPR = rev(cumsum(rev(tab[,2])/sum(tab[,2]))),
             row.names=1:nrow(tab))
}
a <- na.omit(data_CV$cv_lw) 
b <- na.omit(data_CV$cv_hw)

test=c(a,b)
ref=c(rep(0, length(a)), rep(1, length(b)))
res = simple_roc2(ref, test)
res$FNR = 1-res$TPR # 1-Se
pos.optimal.threshold = which.min(res$FPR+res$FNR)
optimal.threshold=row.names(table(test, ref))[pos.optimal.threshold] 

# CRITICAl RAI:
rai = optimal.threshold
rai

###########################################################
#Estimation of Confidence Interval by simulation
###########################################################
rai_subc<-array()

for (i in 1:100000) {
  l<-runif(1,20,length(na.omit(data_CV$cv_hw)))
  subc_CV_HW<-sample(na.omit(data_CV$cv_hw),size =l, replace=F)
  l<-runif(1,20,length(data_CV$cv_lw))
  subc_CV_LW<-sample(data_CV$cv_lw,size =l, replace=F)
  
  dato1<-subc_CV_LW
  dato2<-subc_CV_HW
  
  simple_roc2 <- function(ref, test){
    tab = table(test, ref) # head(tab)
    data.frame(threshold=paste('>=',rownames(tab)), 
               ref0 = tab[,1], 
               ref1 = tab[,2],  
               FPR = rev(cumsum(rev(tab[,1])/sum(tab[,1]))),
               TPR = rev(cumsum(rev(tab[,2])/sum(tab[,2]))),
               row.names=1:nrow(tab))
  }
  a <- na.omit(dato1)  
  b <- na.omit(dato2) 
  
  test=c(a,b)
  ref=c(rep(0, length(a)), rep(1, length(b)))
  res = simple_roc2(ref, test)
  res$FNR = 1-res$TPR # 1-Se
  pos.optimal.threshold = which.min(res$FPR+res$FNR)
  optimal.threshold=row.names(table(test, ref))[pos.optimal.threshold] # Youden threshold  

  # CRITICAL RAI:
  rai_subc[i] = as.numeric(optimal.threshold)
}

rai_subc

rai_subc <- as.numeric(rai_subc)
hist(rai_subc)

#Stats metrics
mean(rai_subc)
median(rai_subc)
min(rai_subc)
max(rai_subc)
quantile(rai_subc)
quantile(rai_subc, probs = 0.975)
quantile(rai_subc, probs = 0.025)

boxplot(rai_subc)