PLS-DA Origano

library(mdatools)
library(kableExtra)
library(tidyverse)
library(pracma) #per detrend
library(plotly)
library(neuralnet)

dftr = read.csv('KS_oregano_TR.csv', sep = ';')
dftr = dftr[,-ncol(dftr)]
dftrdata = dftr[,-c(1:9)]

variabilix = colnames(dftrdata)
variabili = gsub("X", "", colnames(dftrdata))

colnames(dftrdata) = variabili

dfks = read.csv('Origano_KS.csv')
dfks = dfks[-1,]

dftable = table(dftr$Info.contaminante, dftr$Tipologia.reale) %>% data.frame()
dftable %>% pivot_wider(names_from = Var2, values_from = Freq)

## # A tibble: 6 × 3
##   Var1            Adulterant Oregano
##   <fct>                <int>   <int>
## 1 Inflorescences           0     200
## 2 Myrtle pure           1200       0
## 3 Olives pure           1200       0
## 4 Oregano pure             0    4600
## 5 Strawberry pure       1200       0
## 6 Sumac pure            1200       0

dftr %>% group_by('T' = Contaminante) %>% summarise(Count = n())

## # A tibble: 5 × 2
##   T          Count
##   <chr>      <int>
## 1 Myrtle      1200
## 2 Olives      1200
## 3 Oregano     4800
## 4 Strawberry  1200
## 5 Sumac       1200

attr(dftrdata, "xaxis.name") = "Wavelength"
attr(dftrdata, "xaxis.values") = colnames(dftrdata) %>% as.numeric()

mdaplot(dftrdata, type = 'l', main = 'Raw data', cgroup = as.factor(dftr$Contaminante))

dftrdatadetrend = detrend(dftrdata %>% t(), tt = 'linear')

attr(dftrdatadetrend, "xaxis.name") = "Wavelength"
attr(dftrdatadetrend, "xaxis.values") = colnames(dftrdata) %>% as.numeric()


mdaplot(t(dftrdatadetrend), type = 'l', main = 'Detrend', cgroup = as.factor(dftr$Contaminante))

dftrdata2 = scale(t(dftrdatadetrend), scale = F)

attr(dftrdata2, "xaxis.name") = "Wavelength"
attr(dftrdata2, "xaxis.values") = colnames(dftrdata) %>% as.numeric()

mdaplot(dftrdata2, type = 'l', main = 'Detrend + Mean center', cgroup = as.factor(dftr$Contaminante))

PCA

par(mfcol= c(1,1))
pc = pca(dftrdata2, ncomp = 10)

plotVariance(pc)

scorevec = pc$res$cal$scores[,2]
indexvec = c(1:length(scorevec))

dfscore = data.frame('PC2' = scorevec, 'Index' = indexvec)
dfscore = cbind(dftr[,1:8],dfscore)

ggplot(dfscore, aes(x=Index, y=PC2, color = Contaminante))+ geom_point(size = 1)+
  theme_bw()+
  geom_hline(yintercept = 0, linetype ='dashed', alpha = 0.5)

classes = as.numeric(dftr$Tipologia.reale == 'Oregano') %>% factor()

PLS-DA

model = plsda(dftrdata2, dftr$Tipologia.reale, cv = dftr$CV_TR, center = F, lim.type = 'jm')

par(mfrow = c(1,1))
plotRMSE(model)

model = selectCompNum(model , 10)
plotPredictions(model)

predvec = model$res$cv$y.pred[,10,2]+0.5
dfpred = data.frame(index = c(1:length(predvec)), Origano = predvec, Contaminante = dftr$Contaminante, Tipologia = dftr$Tipologia.reale )

ggplot(dfpred, aes(index, Origano, color = Tipologia)) + geom_point()+
  scale_color_manual(values = c('brown3', 'forestgreen'))+
  theme_bw()+
  labs(color = 'Contaminante')+
  scale_y_continuous(breaks=seq(-2,4,0.5))+
  geom_hline(yintercept = 0.5, linetype = 'dashed', alpha = 0.6)

meanorigano = mean(predvec[4801:9600])
meanadulterante = mean(predvec[0:4800])

devorigano = sd(predvec[4801:9600])
devadulterante = sd(predvec[0:4800])

normorigano = rnorm(4800, meanorigano, devorigano)
normadult = rnorm(4800, meanadulterante, devadulterante)
normvec = append(normadult, normorigano)

dfnorm = data.frame(normvec, Tipologia = dftr$Tipologia.reale )

ggplot(dfnorm, aes(normvec)) + 
  #geom_histogram(position = position_dodge2())+
  scale_fill_manual(values = c('brown3', 'forestgreen'))+
  #geom_density(alpha = 0.4)+
  stat_function(fun= dnorm,args = list(mean = meanorigano, sd = devorigano),geom = 'polygon', fill = 'green4', alpha = 0.5)+
  stat_function(fun= dnorm,args = list(mean = meanadulterante, sd = devadulterante),geom = 'polygon', fill = 'brown3', alpha = 0.5)+
  scale_x_continuous(breaks=seq(-2,3.5,0.5))+
  theme_bw()+
  coord_flip(xlim = c(-2,3.5), ylim = c(0,1))+
  xlab('')+
  ylab('')

model$cvres$c.pred[,8,2] %>% table()

## .
##   -1    1 
## 4801 4799

SVM

library('e1071')

## 
## Caricamento pacchetto: 'e1071'

## Il seguente oggetto è mascherato da 'package:pracma':
## 
##     sigmoid

classes = as.numeric(dftr$Tipologia.reale == 'Oregano')

classifier = svm(formula = classes ~ .,
                 data = scale(dftrdata),
                 type = 'C-classification',
                 kernel = 'linear')

y_cv = predict(classifier, newdata = scale(dftrdata))

dfpred = data.frame(index = c(1:length(y_cv)), Origano = y_cv, Contaminante = dftr$Contaminante, Tipologia = dftr$Tipologia.reale )

ggplot(dfpred, aes(index, Origano, color = Tipologia)) + geom_point()+
  scale_color_manual(values = c('brown3', 'forestgreen'))+
  theme_bw()+
  labs(color = 'Contaminante')+
  geom_hline(yintercept = 0, linetype = 'dashed', alpha = 0.4)+
  ggtitle('SVM')

table(classes, y_cv)

##        y_cv
## classes    0    1
##       0 4578  222
##       1  202 4598

Artificial neural network

classes = as.numeric(dftr$Tipologia.reale == 'Oregano')
colnames(dftrdata2) = variabilix

dftrdata3 = cbind('Tipologia' = classes, dftrdata2) %>% as.data.frame()



formula <- as.formula(paste("Tipologia ~ ", paste(variabilix, collapse = " + ")))