# Un altro esempio. Dati da analisi su vini Riesling con tre diverse origini (1=Sudafrica, 2=Germania, 3=Italia, o qualcosa del genere)

wine <- read.table("../dati/wine_tab73_2.dat")
str(wine) # abbiamo 26 osservazioni di 15 variabili (V3 a V17); le prime 2 colonne sono il gruppo ed il numero d'ordine
wine[1:5]
library(MASS)
wine.lda <- lda(wine[,3:17],wine[,1],prior=c(1,1,1)/3) # analisi discriminante lineare con probabilita' a priori uguali
wine.lda
wine.ld <- predict(wine.lda)
wine.ld
plot(wine.ld$x,type="n")
text(wine.ld$x,labels=as.character(wine$V1),col=as.integer(wine$V1))
# una splendida distinzione fra gruppi
# disegno sul grafico anche le medie. Questi sono comandi per farlo
m1 <- mean(subset(wine.ld$x[,1],wine$V1=='1'))
m2 <- mean(subset(wine.ld$x[,2],wine$V1=='1'))
print(c(m1,m2))
points(m1,m2,pch=23,cex=2,col=as.integer(1))# cex indica il tipo di simbolo usato
m1 <- mean(subset(wine.ld$x[,1],wine$V1=='2'))
m2 <- mean(subset(wine.ld$x[,2],wine$V1=='2'))
print(c(m1,m2))
points(m1,m2,pch=23,cex=2,col=as.integer(2))
m1 <- mean(subset(wine.ld$x[,1],wine$V1=='3'))
m2 <- mean(subset(wine.ld$x[,2],wine$V1=='3'))
print(c(m1,m2))
points(m1,m2,pch=23,cex=2,col=as.integer(3))

print(table(wine$V1,wine.ld$class))
# nessun errore come ovvio dal grafico

# passiamo alla cross-validation
wine.ldacv <- update(wine.lda,CV=TRUE) # possiamo usare update anche con questi comandi
wine.ldacv
print(table(wine$V1,wine.ldacv$class))
# 8 errori (su 26 osservazioni). Molti!
# vediamo con un ciclo quali sono le osservazioni errate:
err <- 0
for (i in 1:26) {
	if(wine.ldacv$class[i] == wine$V1[i]) {} else  {err <- err +1;print(i)}
}
err
# guardiamo per esempio il primo errore
wine.ldacv$posterior[5,]
# il metodo da' con quasi certezza una classificazione errata.

# Il problema e' che ci sono troppe variabili (15) e troppo poche osservazioni. 
# Si possono migliorare i risultati, eliminando variabili
# guardiamo la prima variabile discriminante canonica
wine.lda
# rifacciamo l'analisi solo sulla 12,15,16,17

wine.ldarid <- lda (wine$V1 ~ wine$V12+wine$V15+wine$V16+wine$V17, prior = c(1,1,1)/3)
wine.ldarid
wine.ldrid <- predict(wine.ldarid)
print(table(wine$V1,wine.ldrid$class))
# ne sbaglia 3 su 28
wine.ldaridCV <- update(wine.ldarid, CV=TRUE)
print(table(wine$V1,wine.ldaridCV$class))
# anche con la cross-validation abbiamo gli stessi errori.
# usando 4 variabili, il metodo e' molto piu' robusto.


plot(wine.ldrid$x,type="n")
text(wine.ldrid$x,labels=as.character(wine$V1),col=as.integer(wine$V1))