Il documento di analisi spiega tutti i passaggi eseguiti per costruire i vari modelli statistici, che verranno utilizzati per compiti diversi (regressione e classificazione).
Il workflow del progetto si può sintetizzare nel seguente modo:
1. Import e preparazione dei dati
2. Analisi Esplorativa dei dati
3. Analisi di Correlazione tra le variabili
4. Modelli Statistici per compiti di Regressione
5. Modelli Statistici per compiti di Classificazione

DATA MANIPULATION

Il dataset su cui è svolta l’analisi (https://www.kaggle.com/alopez247/pokemon) raccoglie informazioni su tutti i Pokémon usciti fino alla sesta generazione.
Le variabili contenute nel dataset originario sono:

dati = read_excel('pkm.xlsx')
names(dati)
##  [1] "Number"           "Name"             "Type_1"          
##  [4] "Type_2"           "Total"            "HP"              
##  [7] "Attack"           "Defense"          "Sp_Atk"          
## [10] "Sp_Def"           "Speed"            "Generation"      
## [13] "isLegendary"      "Color"            "hasGender"       
## [16] "Pr_Male"          "Egg_Group_1"      "Egg_Group_2"     
## [19] "hasMegaEvolution" "Height_m"         "Weight_kg"       
## [22] "Catch_Rate"       "Body_Style"

Number : è il codice identificativo di ogni Pokémon e si tratta di un numero progressivo, inutile per scopi di analisi.

Name : è il nome del Pokémon, anch’esso inutile per estrarre considerazioni di portata generale.

Type_1 : variabile qualitativa che indica il tipo principale del Pokémon.

Type_2 : alcuni Pokémon sono dotati di 2 tipi, ma come vedremo tra poco la maggior parte di loro ha solo un tipo.

Total : variabile numerica; è la somma di tutte le combat stats che seguono.

HP : hit points; variabile numerica nel range [0,255] che descrive il numero di colpi che un Pokémon può accusare prima di perdere i sensi. Più è alto, più sarà difficile mettere il Pokémon a tappeto.

Attack : variabile numerica nel range [0,255] che descrive la potenza d’attacco del Pokémon. A parità di condizioni e mossa offensiva utilizzata, il Pokémon con valore più alto di Attack toglierà più HP.

Defense : variabile numerica nel range [0,255] che descrive la resistenza del Pokémon. Più è alta, più sarà difficile togliere HP.

Sp_Atk : variabile numerica nel range [0,255]; più è alta, più le mosse speciali del Pokémon considerato saranno dannose per il Pokémon nemico.

Sp_Def : variabile numerica nel range [0,255]; più è alta, meno danno riceverà il Pokémon quando colpito da mosse speciali.

Speed : variabile numerica nel range [0,255]; in duello il Pokémon con valore di Speed maggiore attacca per primo.

Generation : variabile qualitativa con supporto {0,1,..,6} che indica la generazione (in ordine di uscita) a cui appartiene il Pokémon.

isLegendary : variabile qualitativa binaria {0,1} che indica se il Pokémon è Leggendario (=1) oppure no (=0).

Color : il colore del Pokémon; variabile qualitativa con 10 valori diversi, inutile ai fini dell’analisi.

hasGender : variabile binaria; non tutti i Pokémon hanno un genere sessuale, ma come vedremo la maggior parte sì.

Pr_Male : la probabilità che il Pokémon abbia genere maschile. Poiché è una probabilità, si tratta di una variabile numerica con range [0,1].

Egg_Group_1 : variabile qualitativa che indica la razza del Pokémon e determina la possibilità di accoppiamento. Può assumere 15 valori differenti e si è deciso di rimuoverla per non complicare l’analisi.

Egg_Group_2 : come per il Tipo, i Pokémon possono avere due valori di Egg_Group (aumentando il loro spettro di possibili partner).

hasMegaEvolution : variabile binaria che indica se un Pokémon è in grado di mega-evolversi; nel dataset non sono disponibili però informazioni sulle mega-evoluzioni dei Pokémon quindi non verrà considerata nell’analisi.

Height_m : variabile numerica che indica l’altezza (espressa in metri) del Pokémon.

Weight_kg : variabile numerica che indica il peso (espresso in kilogrammi) del Pokémon.

Catch_Rate : variabile numerica nel range [0,255] che indica quanto è facile catturare un Pokémon.

Body_Style : variabile qualitativa che indica l’aspetto anatomico del Pokémon (quadrupede, insettoide, con forma serpentina, bipede con o senza ali ecc.). Può assumere 14 valori differenti e, per non complicare l’analisi, si è deciso di rimuoverla.

Le variabili target sono, rispettivamente per i tasks di regressione e di classificazione, Catch_Rate e isLegendary.

L’obiettivo è quello di rispondere alle seguenti coppie di domande:

1. Quali sono i fattori che spiegano meglio la facilità/difficoltà con cui catturare un Pokémon? Si può misurare l’effetto marginale (ceteris paribus) di questi fattori sul tasso di cattura? Variabile target Catch_Rate 

2. Dato un nuovo Pokémon di cui conosciamo i valori assunti dai predictors considerati, è possibile determinare se esso è un Pokémon leggendario? Variabile target isLegendary

Prima di entrare nella fase di model-building vera e propria, però, si eseguiranno operazioni di manipolazione ed esplorazione dei dati. Sarà durante queste fasi che si potrà rispondere ad altre domande e fare diverse considerazioni, come ad esempio:
- “è comune per un Pokémon avere due tipi e/o due egg_groups?”
- “quanti Pokémon leggendari esistono? Ci sono generazioni con più Pokémon leggendari?”
- “qual è la generazione con più Pokemon in assoluto?” 

# definisco statistiche di combattimento
combat_stats = c('HP','Attack','Defense','Sp_Atk','Sp_Def','Speed')
# trasformo in factor le variabili qualitative e le rinomino
generation = as.factor(dati$Generation)
legendary = as.factor(dati$isLegendary)
type1 = as.factor(dati$Type_1)
type2 = as.factor(dati$Type_2)
# rinomino la variabile target della fase di regressione
catch = dati$Catch_Rate
# dataframe semplificato: solo variabili numeriche
df = cbind(catch,dati[,combat_stats])
# dataframe semplificato esteso con variabili qualitative
df.ext = cbind(generation, legendary, type1, type2, df)

EXPLORATORY DATA ANALYSIS

In questa fase visualizziamo graficamente i dati e facciamo le prime valutazioni di natura qualitativa sulle relazioni tra le variabili.
Vengono fatte le classiche analisi statistiche descrittive univariate: distribuzioni marginali dei valori osservati con barplot per i caratteri discreti ed istogrammi per quelle continue.

# distribuzioni marginali delle variabili type1, type2
t1 <- ggplot(data=df.ext, aes(x=type1, fill=type1)) +
  geom_bar()

t2 <- ggplot(data=df.ext, aes(x=type2, fill=type2)) + 
  geom_bar()

par(mfrow=c(2,1))

t1

t2

# distribuzioni marginali delle variabili generation, legendary
gen <- ggplot(data=df.ext, aes(x=generation, fill=generation)) +
  geom_bar() 

leg <- ggplot(data=df.ext, aes(x=legendary)) +
  geom_bar() 

gen;leg

# quanti leggendari? 
npkm <- data.frame(legendary=sum(df.ext$legendary==1), total=nrow(df.ext),perc=100*sum(df.ext$legendary==1)/nrow(df.ext));
npkm
##   legendary total     perc
## 1        46   721 6.380028

Prime considerazioni:
1. la maggior parte dei Pokémon non ha un tipo 2 (circa la metà delle osservazioni) e questo crea un problema di missing values; per questo motivo si decide di non considerare questa variabile nel resto dell’analisi.
2. sono usciti in media 100 Pokémon per generazione, con il minimo (meno di 75 Pokémon) per la sesta ed il massimo (più di 150 Pokémon) per la quinta.
3. esistono nelle prime sei generazioni 46 Pokémon leggendari, pari a circa il 6% del totale dei Pokémon.

# tbl di contingenza per var. qualitative 
# (generation, type1, legendary)
tab.gen = table(generation); tab.gen
## generation
##   1   2   3   4   5   6 
## 151 100 135 107 156  72
tab.leg = table(legendary); tab.leg
## legendary
##   0   1 
## 675  46
tab.type = table(type1); tab.type
## type1
##      Bug     Dark   Dragon Electric    Fairy Fighting     Fire   Flying 
##       63       28       24       36       17       25       47        3 
##    Ghost    Grass   Ground      Ice   Normal   Poison  Psychic     Rock 
##       23       66       30       23       93       28       47       41 
##    Steel    Water 
##       22      105

Esploriamo quindi le distribuzioni di frequenza dei caratteri quantitativi: istogrammi e curve di densità marginali per le statistiche di combattimento e per la variabile Catch_Rate. Tutte queste variabili assumono valori numerici interi nel set {0,1,…,255} ma la cardinalità del set è sufficientemente alta da trattarle come caratteri quantitativi continui.

# distribuzione marginale di catch 
 print(
    ggplot(data=df, aes(x=df[,1])) +
      geom_density(col='darkblue', fill='lightblue') +
      geom_histogram(aes(y=..density..), alpha=0.4, binwidth = 10) +
      labs(x=names(df)[1]))

# distribuzioni marginali delle statistiche di combattimento
for (i in 2:ncol(df)) {
  print(
    ggplot(data=df, aes(x=df[,i])) +
      geom_density(col='darkblue', fill='lightblue') +
      geom_histogram(aes(y=..density..), alpha=0.4, binwidth = 25) +
      labs(x=names(df)[i])
    )
}

Catch
Graficamente la variabile target per il problema di regressione si presenta con asimmetria positiva : molti valori cadono nella coda destra della curva, con andamento quasi bimodale.

Combat Stats: HP,Attack, Defense, Sp_Atk, Sp_Def, Speed
Tutte le variabili che definiscono l’abilità in combattimento del Pokémon hanno un andamento campanulare e sembrano tutte abbastanza simmetriche. Se ispezioniamo la tabella delle summary statistics, vediamo che - ad eccezione di Catch - i valori di media e mediana (q2) sono molto vicini.

# tbl con summary statistics per var. continue (catch e combat stats)
summary.stats = matrix(0, ncol=6, nrow=7)
colnames(summary.stats) <- c('min','q1','q2','avg','q3','max')
rownames(summary.stats) <- c('catch','hp','atk','def','sp.atk','sp.def','speed')
for (i in 1:7) {
  summary.stats[i,] <- summary(df[,i])
}
summary.stats
##        min q1 q2       avg  q3 max
## catch    3 45 65 100.24688 180 255
## hp       1 50 65  68.38003  80 255
## atk      5 53 74  75.01387  95 165
## def      5 50 65  70.80860  85 230
## sp.atk  10 45 65  68.73786  90 154
## sp.def  20 50 65  69.29126  85 230
## speed    5 45 65  65.71429  85 160

Domanda: Come cambiano le distribuzioni delle variabili considerate, se teniamo in considerazione l’informazione circa la generazione del Pokémon?

Dai grafici che seguono possiamo vedere che le distribuzioni condizionate al gruppo di generazione delle variabili che descrivono l’abilità di combattimento, così come il tasso di cattura, non cambiano significativamente. Non mi aspetto dunque una correlazione significativa tra la variabile Generation e nessuna delle altre variabili.

# distribuzione di catch rate e delle combat stats | generation
print(ggplot(data=df.ext, 
         aes(x=df.ext[,5], col=generation, fill=generation)) +
    geom_density(alpha=0.4) + 
    labs(x=names(df.ext)[5]))

par(mfrow=c(2,3))
for (i in 6:11) {
  print(ggplot(data=df.ext, 
         aes(x=df.ext[,i], col=generation, fill=generation)) +
    geom_density(alpha=0.4) + 
    labs(x=names(df.ext)[i]))
}

Domanda: è vero che i Pokémon Leggendari sono più forti e più difficili da catturare, rispetto a quelli che non sono Leggendari?

Anche per rispondere a questa domanda, si analizzano le distribuzioni delle variabili che raccogliamo nel gruppo combat stats e della variabile catch condizionate alla variabile dicotomica legendary.

# distribuzione di catch rate e delle combat stats | legendary 
print(ggplot(data=df.ext, 
               aes(x=df.ext[,5], col=legendary, fill=legendary)) +
          geom_density(alpha=0.2) + 
          labs(x=names(df.ext)[5]))

for (j in 6:11) {
  print(ggplot(data=df.ext, 
               aes(x=df.ext[,j], col=legendary, fill=legendary)) +
          geom_density(alpha=0.2) + 
          labs(x=names(df.ext)[j]))
} 

# troppi pokemon hanno missing value per la variabile type2,
# perciò la rimuovo dal dataframe esteso con le qualitative
df.ext <- df.ext[,-4]

Chiaramente i Pokémon Leggendari sono più difficili da catturare! La distribuzione azzurra della variabile catch è infatti più schiacciata verso i valori minimi, indicando il gruppo di Pokémon più difficile da catturare con una mera Pokéball.
Per tutte le altre variabili, vediamo che le distribuzioni dei valori di combattimento dei Pokémon leggendari (curve azzurre) sono tutte spostate verso destra rispetto alle analoghe per quelli non leggendari (curve rosse). Possiamo quindi affermare con decisione che i Pokémon leggendari sono sia più difficili da catturare, che più forti in combattimento.

Questo pone le domande centrali per la prossime sezioni:

1. è possibile spiegare il tasso di cattura (catch) data la forza in combattimento del Pokémon? [problema di Regressione; Causation]

2. è possibile prevedere se un Pokémon è leggendario o meno, data la forza in combattimento del Pokémon ed il suo tasso di cattura? [problema di Classificazione; Prediction]

Prima di rispondere a queste domande, però, è necessario studiare la correlazione tra le variabili in esame: se non è vero che il tasso di cattura (Catch) è funzione della forza in combattimento dei Pokémon (combat stats), allora ci aspetteremmo una correlazione molto prossima allo 0.

correlation <- round(cor(df),2)
correlation[upper.tri(correlation)] <- NA

melt.corr <- melt(correlation)

# correlazione tra le variabili numeriche
ggplot(data=melt.corr, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile(color="black") +
  scale_fill_gradient2(midpoint=0, na.value='white', limit=c(-1,1), name='Pearson\nCorrelation') +
  coord_fixed() +
  theme_minimal() +
  geom_text(aes(Var1, Var2, label = value), color = "black", size = 4)

Dalla heatmap dei coefficienti di correlazione, vediamo che la variabile numerica target Catch è modestamente negativamente correlata con tutte le altre variabili che definiscono la forza in combattimento del Pokémon, confermando la nostra intuizione: più il Pokémon in esame è forte, più sarà difficile catturarlo!
Inoltre notiamo come alcune variabili di combattimento siano correlate (positivamente) tra di loro: considerando le stime maggiori o uguali di 0.4, notiamo che le coppie ‘Difesa’ e ‘Difesa Speciale’, ‘Hit Points’ e ‘Attacco’, ‘Attacco Speciale’ e ‘Difesa Speciale’, ‘Attacco Speciale’ e ‘Velocità’ cadono tutte nell’intervallo di modesta correlazione.

spiegare la variabile CATCH [REGRESSION]

Sfruttiamo la correlazione esistente tra le statistiche di combattimento per spiegare la variabilità di Catch Rate. Vengono implementati i seguenti modelli:

1.Modello lineare
Catch ~ Attack, Sp_Atk, Defense, Sp_Defense, HP, Speed
2.Modello lineare con variabile dummy
Catch ~ Attack, … , Speed, legendary
3.Modello polinomiale di 2° grado completo
Catch ~ Attack, Attack^2, AttackSp_Atk, …, AttackSpeed, Sp_Atk, Sp_Atk^2, Sp_AtkDefense, …, HPSpeed, Speed, Speed^2
4.Natural Splines additivi
Catch ~ ns(Attack), … , ns(Speed)
5.Smoothers undimensionali
Catch ~ s(Attack),…,s(Speed)
6.Smoothers bidimensionali con interazione
Catch ~ s(Defense, Sp_Defense), s(Attack,HP), s(Sp_Atk, Speed)
7.Smoothers unidimensionali su Y=log(Catch)
log(Catch) ~ s(Attack),…,s(Speed)
8.Smoother bidimensionali su Y=log(Catch)
log(Catch) ~ s(Defense, Sp_Defense), s(Attack,HP), s(Sp_Atk, Speed)

Il modello lineare multivariato ci servirà come benchmark per confrontare con gli altri modelli e verrà usato per visualizzare graficamente le relazioni tra la variabile target e le variabili esplicative scelte.
Per selezionare il modello migliore si è usato il criterio dell’R quadro aggiustato, in modo da poter confrontare agevolmente modelli con un diverso numero di regressori.

# linear model
mod.a = lm(catch ~ ., data=df)
res.a = summary(mod.a)
r2.a = res.a$adj.r.squared 

# extended linear model
df2 = data.frame(catch=dati$Catch_Rate, atk = dati$Attack,
                 def = dati$Defense, sp_atk = dati$Sp_Atk,
                 sp_def = dati$Sp_Def, hp=dati$HP,
                 speed = dati$Speed, legendary=legendary)


mod.b = lm(catch ~ ., data=df2)
res.b = summary(mod.b)
r2.b = res.b$adj.r.squared

# 2nd degree complete polynomial
df3 = data.frame(catch=df$catch, atk=df$Attack, atk2=df$Attack^2, atk.def=df$Attack*df$Defense,
                 atk.sp_atk=df$Attack*df$Sp_Atk, atk.sp_def=df$Attack*df$Sp_Def, atk.hp=df$Attack*df$HP, atk.spd=df$Attack*df$Speed, sp_atk=df$Sp_Atk, sp_atk2=df$Sp_Atk^2, sp_atk.def=df$Sp_Atk*df$Defense, sp_atk.sp_def=df$Sp_Atk*df$Sp_Def, sp_atk.hp=df$Sp_Atk*df$HP, sp_atk.spd=df$Sp_Atk*df$Speed, def=df$Defense, def2=df$Defense^2, def.sp_def=df$Defense*df$Sp_Def, def.hp=df$Defense*df$HP, def.spd=df$Defense*df$Speed, sp_def=df$Sp_Def, sp_def2=df$Sp_Def^2, sp_def.hp=df$Sp_Def*df$HP, sp_def.spd=df$Sp_Def*df$Speed, hp=df$HP, hp2=df$HP^2, hp.spd=df$HP*df$Speed, spd=df$Speed, spd2=df$Speed^2)

mod.c = lm(catch ~ ., data=df3)
res.c = summary(mod.c)
r2.c = res.c$adj.r.squared

# natural splines 
mod.d = gam(catch ~ ns(HP) + ns(Speed) + ns(Attack) +ns(Sp_Atk) +
              ns(Defense) + ns(Sp_Def), data=df) 
res.d = summary(mod.d)
r2.d = res.d$r.sq 

# unidimensional smoothers
mod.e = gam(catch ~ s(HP) + s(Speed) + s(Attack) + s(Sp_Atk) +
              s(Defense) + s(Sp_Def), data=df) 
res.e = summary(mod.e)
r2.e = res.e$r.sq 

# bidimensional smoothers
mod.f = gam(catch ~ s(HP, Attack) + s(Speed, Sp_Atk) + s(Defense, Sp_Def), data=df) 
res.f = summary(mod.f)
r2.f = res.f$r.sq 

# log catch unidimensional smoothers
mod.g = gam(log(catch) ~ s(HP) + s(Speed) + s(Attack) + s(Sp_Atk) +
              s(Defense) + s(Sp_Def), data=df) 
res.g = summary(mod.g)
r2.g = res.g$r.sq

# log catch bidimensional smoothers
mod.h = gam(log(catch) ~ s(HP, Attack) + s(Speed, Sp_Atk) + s(Defense, Sp_Def), data=df) 
res.h = summary(mod.h)
r2.h = res.h$r.sq 

r2 <- c(r2.a, r2.b, r2.c, r2.d, r2.e, r2.f, r2.g, r2.h); r2
## [1] 0.5415486 0.5426743 0.5654361 0.5415486 0.5604294 0.5767435 0.5919731
## [8] 0.6132918
best.model <- which.max(r2)
result <- list(res.a, res.b, res.c, res.d, res.e, res.f, res.g, res.h)
result[[best.model]]
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(catch) ~ s(HP, Attack) + s(Speed, Sp_Atk) + s(Defense, Sp_Def)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.20993    0.02493   168.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df      F  p-value    
## s(HP,Attack)       9.629  13.33  6.545 7.64e-12 ***
## s(Speed,Sp_Atk)   12.058  16.46 13.605  < 2e-16 ***
## s(Defense,Sp_Def) 13.150  17.66 12.598  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.613   Deviance explained = 63.2%
## GCV = 0.47163  Scale est. = 0.44819   n = 721

Il best model è un modello generalizzato additivo con smoothers unidimensionali e con la variabile target espressa in scala logaritmica; la differenza in termini di devianza spiegata con lo stesso modello con la variabile target espressa nella sua scala di misura è di circa il 3% (differenza significativa). Possiamo quindi concludere che il GAM con smoother unidimensionali e variabile target presa in logaritmi è il miglior modello secondo la misura selezionata, cioè quello in grado di spiegare la maggior parte della variabilità del Catch Rate della popolazione di Pokémon.

Si conclude così l’analisi del problema di regressione, dove l’accento è stato posto sulla stima degli effetti causali e sulla ricerca delle variabili significative (to explain) per spiegare il tasso di cattura di un Pokémon.

Nella parte successiva, invece, il focus sarà la previsione (to predict) della natura Leggendaria o meno del Pokémon a partire da una serie di caratteristiche osservabili: a prescindere dal nesso causale, verranno usati tutti i regressori (e trasformazioni di essi) purché aumentino le performance previsionali.

prevedere la variabile LEGENDARY [CLASSIFICATION]

Nuova variabile target: legendary, variabile qualitativa dicotomica che assume valore 1 se il Pokémon è leggendario.
Obiettivo: sviluppare un algoritmo di classificazione, in grado di prevedere la natura del Pokémon a partire dalle sue features.
Essendo il focus la previsione, si devono creare 2 partizioni del dataset originario ed usare una di queste (training set) per addestrare il modello, l’altra (test set) per testarlo e calcolare la misura di performance. Vengono usati circa 2/3 delle osservazioni come insieme di addestramento del modello (regola di split).

I modelli di classificazione utilizzati sono:
1. Regressione Logistica
2. Analisi Discriminante Lineare (LDA)
3. Analisi Discriminante Quadratica (QDA)
4. Albero di Classificazione
5. Boosting con Stumps
6. Random Forest

Come misura di performance verrà usato il tasso di errore di classificazione, definito come la somma di falsi positivi e falsi negativi (cioè contiamo, rispettivamente e con pari peso, sia un Pokémon erroneamente classificato come Leggendario, che uno erroneamente classificato come Non Leggendario).

set.seed(13)
# training and test sets  
sample <- sample.int(size=0.7*nrow(df), n=nrow(df))
train <- df[sample, ]
test <- df[-sample, ]
train <- data.frame(legendary=legendary[sample], train)
test.y <- legendary[-sample]

# logistic regression
mod_a <- glm(legendary ~ ., data=train, family = 'binomial', control=list(maxit=50))
probs <- plogis(predict(mod_a, newdata=test))
preds.a <- ifelse(probs>=0.5, 1, 0)
error_a <- ifelse(test.y==preds.a,0,1)
error_rate_a <- sum(error_a)/length(error_a); 
error.a <- sum(error_a)
conf_matrix_a <- table(test.y, yhat=preds.a)
conf_matrix_a
##       yhat
## test.y   0   1
##      0 199   2
##      1   1  15
# linear discriminant analysis 
mod_b <- lda(legendary ~ ., data=train, method='mle')
preds.b <- predict(mod_b, newdata=test, type='class')
error_b <- ifelse(test.y==preds.b$class, 0, 1)
error_rate_b <- sum(error_b)/length(error_b)
error.b <- sum(error_b)
conf_matrix_b <- table(test.y, yhat=preds.b$class)
conf_matrix_b
##       yhat
## test.y   0   1
##      0 201   0
##      1  11   5
# quadratic discriminant analysis 
mod_c <- qda(legendary ~ ., data=train, method='mle')
preds.c <- predict(mod_c, newdata=test, type='class')
error_c <- ifelse(test.y==preds.c$class, 0, 1)
error_rate_c <- sum(error_c)/length(error_c)
error.c <- sum(error_c)
conf_matrix_c <- table(test.y, yhat=preds.c$class)
conf_matrix_c
##       yhat
## test.y   0   1
##      0 199   2
##      1   5  11
# classification tree
mod_d <- rpart(legendary ~ ., method="class", data=train)
preds.d <- predict(mod_d, newdata=test, type='class')
error_d <- ifelse(test.y==preds.d, 0, 1)
error.d <- sum(error_d)
error_rate_d <- sum(error_d)/length(error_d)
conf_matrix_d <- table(test.y, yhat=preds.d)
rpart.plot(mod_d)

conf_matrix_d
##       yhat
## test.y   0   1
##      0 200   1
##      1   5  11
# boosting with stumps (d=1)
lambda <- 0.01
d <- 1
cv <- trainControl(method ="repeatedcv",number=5,repeats=2)
grid <- expand.grid(n.trees = c(100, 500, 1000), shrinkage=lambda, interaction.depth=d, n.minobsinnode=c(5,10))
boost <- train(legendary~., train, method = "gbm",
  tuneGrid = grid, verbose = FALSE, trControl=cv)
preds.e <- predict(boost, newdata=test)
error_e <- ifelse(test.y==preds.e,0,1)
conf_matrix_e <- table(test.y, yhat=preds.e)
error.e <- sum(error_e)
error_rate_e <- sum(error_e)/length(error_e)
conf_matrix_e
##       yhat
## test.y   0   1
##      0 199   2
##      1   2  14
# random forest 
rf <- train(legendary ~., train, method="rf", trControl=cv)
preds.f <- predict(rf, newdata=test)
error_f <- ifelse(test.y==preds.f, 0, 1)
error.f <- sum(error_f)
conf_matrix_f <- table(test.y, yhat=preds.f)
error_rate_f <- sum(error_f)/length(error_f)
conf_matrix_f
##       yhat
## test.y   0   1
##      0 199   2
##      1   6  10
# confronto : scelgo il modello con error_rate minimo 
d <- data.frame(error_rates=c(error_rate_a, error_rate_b, error_rate_c, error_rate_d, error_rate_e, error_rate_f), error=c(error.a,error.b,error.c,error.d,error.e,error.f), model=c('logistic', 'lda', 'qda', 'tree','boosting','rf'))
d
##   error_rates error    model
## 1  0.01382488     3 logistic
## 2  0.05069124    11      lda
## 3  0.03225806     7      qda
## 4  0.02764977     6     tree
## 5  0.01843318     4 boosting
## 6  0.03686636     8       rf

Il modello con il minor numero di errori di classificazione è quello di regressione logistica, in particolare dalla matrice di confusione vediamo che solo 1 Pokémon Leggendario è stato erroneamente classificato (falso negativo) e che solo 2 Pokémon non leggendari sono stati scambiati per Leggendari dall’algoritmo (falsi positivi).
La coppia di regressori in magnitudine più incisivi sono le variabili Catch e Attack, e questo risultato si vede anche nell’Albero di Classificazione presentato sopra.

Possiamo quindi concludere affermando che le caratteristiche distintive dei Pokémon Leggendari sono la forza di attacco e l’estrema difficoltà nel catturarli.

mod_a
## 
## Call:  glm(formula = legendary ~ ., family = "binomial", data = train, 
##     control = list(maxit = 50))
## 
## Coefficients:
## (Intercept)        catch           HP       Attack      Defense  
##   -33.13773     -0.14958      0.01630      0.09668      0.05804  
##      Sp_Atk       Sp_Def        Speed  
##     0.06802      0.07229      0.04418  
## 
## Degrees of Freedom: 503 Total (i.e. Null);  496 Residual
## Null Deviance:       227.5 
## Residual Deviance: 26.24     AIC: 42.24
conf_matrix_a
##       yhat
## test.y   0   1
##      0 199   2
##      1   1  15

Bibliografia

  1. https://www.kaggle.com/alopez247/pokemon Dataset
  2. https://it.wikipedia.org/wiki/Pok%C3%A9mon Informazioni generali sui Pokémon
  3. https://it.wikipedia.org/wiki/Analisi_discriminante Modelli di analisi discriminante lineare e quadratica
  4. https://it.wikipedia.org/wiki/Modello_logit Modello di regressione logistica
  5. https://it.wikipedia.org/wiki/Foresta_casuale Modello Random Forest
  6. https://en.wikipedia.org/wiki/AdaBoost Modello di Boosting per variabile dicotomica