library(pls)

##################################################################################
#
#                       PCR (Principal Component Regression)
#
##################################################################################

### Carico il dataset

#Description
#A data set with NIR spectra and octane numbers of 60 gasoline samples.
#The NIR spectra were measured using diffuse reflectance as log(1/R)
#from 900 nm to 1700 nm in 2 nm intervals, giving 401 wavelengths.
#http://127.0.0.1:14596/library/pls/html/gasoline.html

originale<-read.table("Dataset gasoline.txt",header=TRUE,sep="\t")
View(originale)

#-- Guardo i dati in grafico:
plot(seq(900,1700,2),originale[1,-1],type="l", xlab="nm",ylab="signal")

Color<-rainbow(nrow(originale))
for (i in c(2:nrow(originale))) {lines(seq(900,1700,2),originale[i,-1],col=Color[i])}

#La risposta è la colonna "octane", le altre colonne sono i predittori

### Bisogna dividere il dataset in Training e Test dataset:

Train<-originale[1:42,] # 70%
Test<-originale[43:60,] # 30%

### Costruisco il modello: prima bisogna dividere i predittori dalla risposta

TrainP<-as.matrix(Train[,c(2:402)])# Predittori
TrainR<-as.matrix(Train[,1]) # Risposta

#-- per costruire il modello uso la cross validation Leave-One-Out (LOO)

Model <- pcr(TrainR ~ TrainP, ncomp = 10, validation = "LOO")

summary(Model)

#-- valuto il numero di componenti di interesse per il modello

plot(RMSEP(Model), legendpos = "topright") # dalla terza componente in poi si abbassa drasticamente
                                           # l'errore del modello calcolato con cross validation LOO

##--Quindi proseguo con la validazione esterna sul Test.
##--divido anche il Test dataset in predittori e risposta:

TestP<-as.matrix(Test[,c(2:402)])# Predittori
TestR<-as.matrix(Test[,1]) # Risposta

Predizione<-predict(Model, ncomp=3,newdata=TestP) # 3 componenti, scelte in precedenza

plot(TestR,Predizione[,,1],pch=16,col="blue",xlab="measured",ylab="predicted",asp=1)
abline(a=0,b=1) # retta y=x

##-- Calcolo R^2:

summary(lm(Predizione[,,1]~TestR)) # R^2 è 0.98


##-- Se avessi scelto solo 2 componenti?(completare)

Predizione2<-predict(Model, ncomp=      ,newdata=TestP)

plot(TestR,          [,,1],pch=16,col="blue",xlab="measured",ylab="predicted",asp=1)
abline(a=0,b=1) # retta y=x
points(TestR,Predizione[,,1],pch=16,col="red")# questi sono i valori di predizione con 3 componenti

summary(lm(           [,,1]~TestR)) # R^2 è 


##################################################################################
#
#                       PLS (Partial Least Square regression)
#
##################################################################################

### Generare un modello PLS con gli stessi dati,
### sempre con cross validation Leave-One-Out (LOO)

Model2 <- plsr(TrainR ~ TrainP, ncomp = 10, validation = "LOO")

summary(Model2)


#-- valuto il numero di componenti di interesse per il modello

plot(RMSEP(Model2), legendpos = "topright") # Quante componenti sono necessarie? 2


##--Quindi proseguo con la validazione esterna sul Test.

PredizionePLS<-predict(Model2, ncomp=  2   ,newdata=TestP)

plot(TestR,PredizionePLS[,,1],pch=16,col="blue",xlab="measured",ylab="predicted",asp=1)
abline(a=0,b=1) # retta y=x

##-- Calcolo R^2:

summary(lm(PredizionePLS[,,1]~TestR)) # R^2 è 0.98 ma, rispetto alla PCR, usando solo 2 componenti