Consideriamo la tabella che contiene i dati (opportunamente elaborati) sulle immagini ricavate dalla scansione di banconote, al fine di riconoscere le banconote false, scaricato da https://archive.ics.uci.edu/ml/datasets/banknote+authentication.
Ai fini dell’implementazione del modello regressivo elementare, rileggiamo la classe in termini di \(\pm1\).
Proviamo una regressione lineare (classica).
Quante osservazioni del training set sono corrette? Confrontiamo i valori della classe con quelli previsti, calcolando l’accuratezza, attraverso una classificazione perentoria,
Può risultare interessante scoprire quali errori sono stati compiuti. Visualizziamo gli errori in una matrice di confusione.
source("s2_cmroc.r")
response=(bn$class>0)
predictor=(predict(bn.lm)>0)
s2_pconfusion(response,predictor)Proviamo a verificare la capacità di predizione.
idx=sample(1372,50)
bncv=bn[-idx,]
bncv.lm=lm(class~.,data=bncv)
predictor=(predict(bncv.lm,bn[idx,])>0)
response=(bn$class[idx]>0)
sum(predictor==response)/50Il risultato non è statisticamente significativo. Ripetiamo più volte l’esperimento su campioni casuali (ovvero stimiamo l’errore con l’autovalutazione).
acc=rep(0,50)
for(i in 1:50){
idx=sample(1372,50)
bncv=bn[-idx,]
bncv.lm=lm(class~.,data=bncv)
predictor=(predict(bncv.lm,bn[idx,])>0)
response=(bn$class[idx]>0)
acc[i]=sum(predictor==response)/50
}
mean(acc)
sd(acc)
hist(acc,10)Proviamo a studiare quanto è robusto il modello, introducendo informazioni false per studiare l’andamento dell’accuratezza.
Consideriamo la regressione logistica sulla stessa tabella.
Valutiamo la diagnostica: accuratezza e matrice di confusione
bn.glm.p=predict(bn.glm,type="response")
sum((bn.glm.p>0.5)==(bn$class>0.5))/length(bn$class)
s2_confusion(bn$class,bn.glm.p)Non è detto che in generale la soglia 1/2 sia la migliore.
In una classificazione {0,1} (dove 1 è positivo, e 0 è negativo) la curva ROC rappresenta la proporzione di risposte corrette rispetto alla soglia di differenziazione. In ascissa c’è la specificità, ovvero il tasso di negativi corretti, in ordinata c’è la sensibilità, ovvero il tasso di positivi corretti, \[ \begin{aligned} \text{specificità} &= \frac{\text{negativi corretti}}{\text{negativi}}\\ \text{sensibilità} &= \frac{\text{positivi corretti}}{\text{positivi}}\\ \end{aligned} \]
Consideriamo la tabella degli indicatori socio-economici di alcuni paesi europei.
Aggiungiamo l’informazione della classe, decisa a-priori tra paesi nord-occidentali e paesi sud-orientali.
La regressione logistica ha problemi di convergenza.
Dobbiamo scegliere un modello più economico. Possiamo utilizzare almeno tre metodi differenti:
Semplifichiamo il modello selezionando i fattori a più alta correlazione con la ``classe’’, ovvero Employment, Export, TaxesTotal e Accidents
Dopo la riduzione del modello, procediamo con la regressione logistica.
I punteggi e le probabilità delle singole nazioni, il numero di assegnazioni errate, l’accuratezza.
round(predict(eu.glm),2)
eu.glm.p=predict(eu.glm,type="response")
sum((eu.glm.p>0.5)==(eu$classe>0.5))
sum((eu.glm.p>0.5)==(eu$classe>0.5))/length(eu$classe)Matrice di confusione e curva ROC.
Usando la regressione multivariata si mantengono i fattori Employment, TaxesTotal, Accidents.
Standardizziamo i dati, in vista dell’analisi delle componenti principali.
Usando l’analisi delle componenti principali, la prima componente è principalmente allineata con Education, HealthyYears e IndustrialProduction Il risultato non è entusiasmante.
loadings(eu.pca)
eu.pca1.glm=glm(classe~Education+HealthyYears+IndustrialProduction,family=binomial,data=eu)
eu.pca1.p=predict(eu.pca1.glm,type="response")
s2_confusion(eu$classe,eu.pca1.p)Oppure, si aggiunge la classe e si vede che classe è allineata principalmente con Employment, Export e TaxesLabour. Il risultato sarà il migliore perché si usa l’informazione relativa alla classe.
eus$classe<-a
eus.pca=princomp(eus)
biplot(eus.pca)
eu.pca2.glm=glm(classe~Employment+Export+TaxesLabour,family=binomial,data=eu)
eu.pca2.p=predict(eu.pca2.glm,type="response")
s2_confusion(eu$classe,eu.pca2.p)Infine, si usano direttamente le (prime tre) componenti principali, ottenendo risultati paragonabili. Il risultato migliora (un po’) aumentando il numero di componenti.
Consideriamo il dataset relativo ai passeggeri del Titanic.
Consideriamo un modello di analisi discriminante lineare allo scopo di prevedere la sopravvivenza dei passeggeri e dell’equipaggio del Titanic.
Valutiamo l’efficacia del modello.
tit.lda.p=predict(tit.lda)
tit.lda.post=tit.lda.p$posterior[,2]
sum((tit.lda.post>0.5)==(tit$survived>0.5))/length(tit$survived)
s2_confusion(tit$survived,tit.lda.post)
tit.lda.roc=s2_roc(tit$survived,tit.lda.post)
s2_roc.plot(tit.lda.roc)Esaminiamo anche il risultato dell’analisi discriminante quadratica.
tit.qda=qda(survived~.,data=tit,CV=F)
tit.qda.p=predict(tit.qda)
tit.qda.post=tit.qda.p$posterior[,2]
sum((tit.qda.post>0.5)==(tit$survived>0.5))/length(tit$survived)
s2_confusion(tit$survived,tit.qda.post)Confrontiamo le due curve ROC.
tit.qda.roc=s2_roc(tit$survived,tit.qda.post)
s2_roc.plot(tit.qda.roc,col="blue")
s2_roc.lines(tit.lda.roc,col="green3")Un indice di confronto interessante è il parametro AUC, ovvero l’area sottostante la curva ROC (migliore se più vicina a 1, peggiore se più vicina a 0).
Aggiungiamo anche il modello di regressione logistica.
I fattori pclass, sex, age sono i più significativi per decidere sulle possibilità di sopravvivenza. Volendo si possono trascurare i fattori che sembrano meno significativi. Non è ovvio se questo dia un beneficio.
summary(glm(survived~.-parch,family=binomial,data=tit))
summary(glm(survived~.-parch-fare,family=binomial,data=tit))Esaminiamo accuratezza e matrice di confusione.
tit.glm.p=predict(tit.glm,type="response")
sum((tit.glm.p>0.5)==(tit$survived>0.5))/length(tit$survived)
s2_confusion(tit$survived,tit.glm.p)Confrontiamo le curve ROC e le aree sotto la curva.
tit.glm.roc=s2_roc(tit$survived,tit.glm.p)
s2_roc.plot(tit.qda.roc,col="blue")
s2_roc.lines(tit.lda.roc,col="green3")
s2_roc.lines(tit.glm.roc,col="red")
s2_auc(tit.glm.roc)Stabiliamo attraverso l’autovalutazione il miglior classificatore.
l=length(tit$survived)
acc=matrix(0,100,3)
for(i in 1:100){
idx=sample(l,100)
titcv=tit[-idx,]
tit.lda=lda(survived~.,data=titcv)
tit.lda.p=predict(tit.lda,tit[idx,])$posterior[,2]
acc[i,1]=sum((tit.lda.p>0.5)==(tit$survived[idx]>0.5))/100
tit.qda=qda(survived~.,data=titcv)
tit.qda.p=predict(tit.qda,tit[idx,])$posterior[,2]
acc[i,2]=sum((tit.qda.p>0.5)==(tit$survived[idx]>0.5))/100
tit.glm=glm(survived~.,family=binomial,data=titcv)
tit.glm.p=predict(tit.glm,tit[idx,],type="response")
acc[i,3]=sum((tit.glm.p>0.5)==(tit$survived[idx]>0.5))/100
}
# lda
mean(acc[,1])
sd(acc[,1])
# qda
mean(acc[,2])
sd(acc[,2])
# reg. logistica
mean(acc[,3])
sd(acc[,3])irisSvolgiamo una analisi discriminante lineare sul dataset iris.
library(MASS)
data(iris)
iris.lda=lda(Species~.,data=iris,prior=c(1/3,1/3,1/3))
plot(iris.lda,col=as.numeric(iris$Species)+1)Una rappresentazione grafica più chiara.
Valutiamo l’accuratezza.
sum(iris$Species!=iris.lda.values$class)
sum(iris$Species!=iris.lda.values$class)/length(iris$Species)Confrontiamo con la capacità distintiva dell’analisi delle componenti principali
iris.pca=princomp(scale(subset(iris,select=-Species)))
iris.pca.pt=predict(iris.pca)
plot(iris.pca.pt,col=as.numeric(iris$Species)+1,pch=20)Proviamo la capacità predittiva, scegliendo tre campioni a caso, uno per specie.
idx=round(runif(3,1,50))+c(0,50,100)
irtrain=iris[-idx,]
irtest=iris[idx,]
irtrain.lda=lda(Species~.,data=irtrain,prior=c(1/3,1/3,1/3))
irtest.pt=predict(irtrain.lda,newdata=irtest)La classe cui appartengono i campioni scelti e le probabilità di appartenenza alle classi di ogni campione del test set.
Una autovalidazione più robusta.
Una modalità per utilizzare la regressione logistica anche nel caso multiclasse è di confrontare ogni classe con tutte le altre, e poi scegliere per ogni osservazione scegliere la classe più probabile.
Regressione logistica per classificare setosa.
ir=iris[,1:4]
ir$class<-as.numeric(iris$Species=="setosa")
set.glm=glm(class~.,family=binomial,data=ir)
set.p=predict(set.glm,type="response")Regressione logistica per classificare versicolor
ir$class<-as.numeric(iris$Species=="versicolor")
ver.glm=glm(class~.,family=binomial,data=ir)
ver.p=predict(ver.glm,type="response")Regressione logistica per classificare virginica
ir$class<-as.numeric(iris$Species=="virginica")
vir.glm=glm(class~.,family=binomial,data=ir)
vir.p=predict(vir.glm,type="response")Calcoliamo la risposta finale e l’accuratezza.
Generare una tabella con 30 righe e 4 colonne con dati casuali. Aggiungere una quinta colonna contenente una classificazione binaria, generata anche essa casualmente. Calcolare la discrepanza tra la classificazione originaria e il risultato di una classificazione logistica. Se si ripete questo esperimento per 1000 volte, calcolare valore medio e deviazione standard della discrepanza.
Generare una tabella A con 170 righe e 7 colonne, in modo che l’ultima colonna contenga valori binari relativi all’appartenenza a una o l’altra di due classi. Generare una seconda tabella B con 5 righe e 6 colonne. Implementare una classificazione degli individui nella tabella B applicando alla tabella A la classificazione logistica e la analisi discriminante lineare. Confrontare i risultati dei due metodi.
Confrontare i metodi di regressione logistica e analisi discriminante lineare e quadratica sul dataset 12titanic usando l’autovalidazione, per stabilire il metodo migliore.
Confrontare il modello di regressione logistica completo e i modellio ridotti, il primo senza parch, il secondo senza parch e fare, per verificare se c’è una sostanziale differenza nei modelli (calcolare accuratezza, AUC, matrici di confusione e confrontare le curve ROC).
Determinare il migliore metodo di classificazione per il dataset Pima.tr. Usare poi i modelli per classificare i nuovi dati contenuti nel dataset Pima.te e verificare se il miglior modello stabilito nella prima parte sia effettivamente il più efficace.