Nell’articolo precedente abbiamo affrontato la regressione lineare calcolata su una proprietà (il numero delle stanze di un’abitazione). Però, realisticamente, ci sono numerose altre caratteristiche di un’abitazione che concorrono a definirne il valore di mercato. Ad esempio la sua distanza dal centro città o dai principali servizi (uffici pubblici, banche, metropolitana, supermercati, eccetera), o la superficie, o il numero di piani su cui è sviluppata, se è dotata di garage autonomo, la presenza di un giardino, l’anno di costruzione, eccetera.
In questo caso è necessario utilizzare la sorella maggiore della regressione lineare già vista: la regressione lineare multipla o multidimensionale.
Non vi spaventate, non vi annoierò con altra teoria di matematica e geometria, perchè quello che è stato detto nell’articolo precedente è ancora valido. Semplicemente, considerate di fare una regressione lineare separata per ognuna delle caratteristiche ed alla fine mettete tutto insieme con una bella somma. Il risultato sarà:
Visto che i bias sono semplici scalari, li posso sommare tra loro ottenendo un unico scalare che chiamo semplicemente b. Già che ci sono, raggruppo tutto con una bella sommatoria, per semplificare la lettura:
che come è evidente, è del tutto simile a quanto visto nella regressione lineare semplice.
Cosa cambia? Se consideriamo due caratteristiche (ad esempio superficie ed anno di costruzione), avremo che il sistema si sviluppa in uno spazio tridimensionale (la terza dimensione è il valore dell’abitazione). Quindi la nostra regressione non sarà più rappresentata da una retta, ma da un piano che minimizzerà le distanze tra lui ed i punti del training set.
Se le caratteristiche su cui fare la stima fossero ancora di più (aggiungendo ad esempio anche il numero di stanze, la distanza dal centro città, eccetera), le dimensioni complessive del sistema sarebbero superiori a 3 e quindi inizieremmo a parlare più in generale di iperspazio del sistema e la regressione sarebbe rappresentata da un iperpiano.
No, non preoccupatevi: non stiamo parlando di StarTrek e di viaggi interplanetari… E come avete intuito, le analogie con la regressione lineare semplice sono tantissime. Se ricordate il codice sorgente dell’articolo precedente, la funzione di regressione lineare di SciKitLearn prevede già i pesi w della regressione sotto forma di array. Questo vi dice qualcosa?
Esempio in Python
Bene, vediamo il nuovo esempio in Python per una regressione lineare multipla.
NOTA: il codice sorgente di questo esempio lo potete trovare qui
In verità non cambia molto dal codice della volta scorsa, per cui evidenzieremo solo i passaggi cardine. Semplicemente andremo a prendere tutte le colonne presenti nel dataset, utilizzando la nomenclatura suggerita dal fornitore dei dati.
# importiamo il dataset direttamente dalla URL dove è archiviato
dataset = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data",
# indichiamo che il file utilizza un numero indefinito di spazi come separatore di colonna
sep="\s+",
# assegniamo i nomi alle colonne (questa volte le utilizziamo tutte!)
# utilizzo la nomenclatura suggerita dallo stesso fornitore
# che si trova a questo indirizzo
# https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names
#
# 0. CRIM per capita crime rate by town
# 1. ZN proportion of residential land zoned for lots over 25,000 sq.ft.
# 2. INDUS proportion of non-retail business acres per town
# 3. CHAS Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
# 4. NOX nitric oxides concentration (parts per 10 million)
# 5. RM average number of rooms per dwelling
# 6. AGE proportion of owner-occupied units built prior to 1940
# 7. DIS weighted distances to five Boston employment centres
# 8. RAD index of accessibility to radial highways
# 9. TAX full-value property-tax rate per $10,000
# 10. PTRATIO pupil-teacher ratio by town
# 11. B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
# 12. LSTAT % lower status of the population
# 13. MEDV Median value of owner-occupied homes in $1000's
names=["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO", "B", "LSTAT", "MEDV"])
Diamo uno sguardo al dataset...
CRIM ZN INDUS CHAS NOX … TAX PRATIO B LSTAT MEDV
0 0.00632 18.0 2.31 0 0.538 … 296.0 15.3 396.90 4.98 24.0
1 0.02731 0.0 7.07 0 0.469 … 242.0 17.8 396.90 9.14 21.6
2 0.02729 0.0 7.07 0 0.469 … 242.0 17.8 392.83 4.03 34.7
3 0.03237 0.0 2.18 0 0.458 … 222.0 18.7 394.63 2.94 33.4
4 0.06905 0.0 2.18 0 0.458 … 222.0 18.7 396.90 5.33 36.2
Dando uno sguardo ai tipi di dato delle varie colonne, ci accorgiamo che tutte le colonne sono valorizzate con numeri reali (float) tranne due che sono valorizzate con numeri interi. Queste due (CHAS e RAD) sono due caratteristiche categoriche.
# diamo uno sguardo ai tipi dati
print("Diamo uno sguardo ai tipi dato...")
print(dataset.info())
Diamo uno sguardo ai tipi dato...
RangeIndex: 506 entries, 0 to 505
Data columns (total 14 columns):
CRIM 506 non-null float64
ZN 506 non-null float64
INDUS 506 non-null float64
CHAS 506 non-null int64
NOX 506 non-null float64
RM 506 non-null float64
AGE 506 non-null float64
DIS 506 non-null float64
RAD 506 non-null int64
TAX 506 non-null float64
PRATIO 506 non-null float64
B 506 non-null float64
LSTAT 506 non-null float64
MEDV 506 non-null float64
- La prima (CHAS) indica se l’abitazione è in prossimità del fiume Charles (i dati sono della città di Boston) e vale 1 in questo caso, 0 altrimenti.
- La seconda (RAD) indica il livello di accessibilità dell’abitazione verso le arterie autostradali della città.
Ora, di tutte queste proprietà, quali sono le più interessanti e che si legano maggiormente al valore di uscita (il valore commerciale medio dell’abitazione)? Se non siete data scientist vi suggerisco di seguire gli articoli tematici del collega Fabio Celli. Nel frattempo possiamo sfruttare le funzionalità di Pandas che possono aiutarci a fare un’analisi empirica del grado di dipendenza delle proprietà.
# valutiamo il grado di correlazione tra le proprietà del dataset
print("Grado di correlazione")
print(dataset.corr())
Grado di correlazione
CRIM ZN INDUS … B LSTAT MEDV
CRIM 1.000000 -0.200469 0.406583 … -0.385064 0.455621 -0.388305
ZN -0.200469 1.000000 -0.533828 … 0.175520 -0.412995 0.360445
INDUS 0.406583 -0.533828 1.000000 … -0.356977 0.603800 -0.483725
CHAS -0.055892 -0.042697 0.062938 … 0.048788 -0.053929 0.175260
NOX 0.420972 -0.516604 0.763651 … -0.380051 0.590879 -0.427321
RM -0.219247 0.311991 -0.391676 … 0.128069 -0.613808 0.695360
AGE 0.352734 -0.569537 0.644779 … -0.273534 0.602339 -0.376955
DIS -0.379670 0.664408 -0.708027 … 0.291512 -0.496996 0.249929
RAD 0.625505 -0.311948 0.595129 … -0.444413 0.488676 -0.381626
TAX 0.582764 -0.314563 0.720760 … -0.441808 0.543993 -0.468536
PRATIO 0.289946 -0.391679 0.383248 … -0.177383 0.374044 -0.507787
B -0.385064 0.175520 -0.356977 … 1.000000 -0.366087 0.333461
LSTAT 0.455621 -0.412995 0.603800 … -0.366087 1.000000 -0.737663
MEDV -0.388305 0.360445 -0.483725 … 0.333461 -0.737663 1.000000
In pratica viene creata una tabella con tutte le proprietà in riga e colonna, e per ogni coppia viene calcolato il grado di correlazione tra i valori. Se il risultato si avvicina a +1, c’è una forte correlazione diretta tra le due proprietà che incrociano quella cella. Se il risultato si avvicina a -1, c’è una forte correlazione inversa. Se il risultato si avvicina a 0, la correlazione sarà nulla o scarsa.
Certo che però leggersi una tabella di decine e decine di valori numerici non è sicuramente facile. Quindi possiamo esplicitare graficamente questi risultati attraverso una cosiddetta heat map, cioè una mappa bidimensionale che associa sfumature di colore ai vari valori delle celle, in modo da rendere più facile identificare le celle che corrispondono a forti correlazioni.
Per fare questo, ovviamente in Python esiste la libreria specifica che ci facilita il compito. Questa è Seaborn, che si appoggia a MatPlotLib per creare i grafici.
import seaborn as sns
# creiamo una heatmap delle correlazioni con seaborn
# i parametri xticklabels e yticklabels servono a fornire il nome di righe e colonne
# forniremo per entrambe l'elenco dei nomi colonna del dataset
sns.heatmap(dataset.corr(), xticklabels=dataset.columns, yticklabels=dataset.columns)
plt.show()
Come si vede, le relazioni dirette hanno colori più chiari tendenti al bianco, mentre le relazioni inverse hanno colori più scuri tendenti al nero. I colori centrali indicano relazioni meno accentuate o nulle. Appare evidente la diagonale bianca in quanto è ovvio che la correlazione tra una proprietà e sè stessa è assolutamente la più alta possibile.
Osservando l’ultima riga della proprietà MEDV, che è il nostro valore di uscita, notiamo che la proprietà RM (numero di stanze) ha una correlazione diretta abbastanza accentuata, così come la proprietà LSTAT (grado di povertà della popolazione della zona) ha una forte correlazione inversa. E in effetti intuitivamente torna: più alto è il numero di stanze di un’abitazione, più il valore commerciale crescerà, mentre più alto è l’incidenza della povertà in zona, meno valore avrà l’abitazione. Dalla heatmap pare anche che le altre proprietà siano correlate in modo più debole con il valore di uscita che stiamo considerando.
Per ulteriore verifica possiamo utilizzare un altro strumento di seaborn che per ogni combinazione di coppie di proprietà, esegue un grafico di dispersione xy. Visto che le proprietà sono molte e quindi il grafico verrebbe enorme, selezioniamo solo quelle più promettenti dall’heatmap di poco fa, aggiungendo anche INDUS, TAX e PRATIO.
# creiamo i grafici delle combinazioni di coppie di proprietà
# NOTA: visto che il grafico diventerebbe enorme, selezioniamo solo le colonne che reputiamo
# più promettenti dopo aver analizzato la heatmap precedente, oltre al nostro valore di uscita
sns.pairplot(dataset[["RM", "LSTAT", "PRATIO", "TAX", "INDUS", "MEDV"]])
plt.show()
Possiamo notare che, a parte la solita diagonale che non ci interessa, sulla riga di MEDV gli unici grafici con una nuvola di valori che si presta ad una regressione lineare sono quelli della proprietà RM e LSTAT, cosa che conferma ulteriormente l’intuizione della heatmap. A dire il vero LSTAT ha una tendenza leggermente curva, ma per ora possiamo approssimare ad un andamento lineare. Nella immagine seguente ho fatto un ingrandimento dei due grafici e disegnato manualmente due righe rosse per chiarire il concetto. Non è ancora una regressione quella!
# associamo ad X i valori di input delle colonne che reputiamo più promettenti
# associamo ad Y i valori di output
X = dataset[["RM", "LSTAT"]].values
Y = dataset["MEDV"].values
# suddividiamo il dataset in due dataset, uno di training ed uno di test
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3)
# istanziamo la classe di calcolo della regressione lineare di SciKitLearn
# la addestriamo e prediciamo i valori con il set di test
lRegr = LinearRegression()
lRegr.fit(X_train, Y_train)
Y_pred = lRegr.predict(X_test)
# calcoliamo l'errore quadratico medio e il coefficiente di determinazione
errore = mean_squared_error(Y_test, Y_pred)
print("Errore:", errore)
punteggio = r2_score(Y_test, Y_pred)
print("Score:", punteggio)
print()
# visualizziamo i valori dei pesi e del bias trovati
print("Valore del peso di RM:", lRegr.coef_[0])
print("Valore del peso di LSTAT:", lRegr.coef_[1])
print("Valore del bias:", lRegr.intercept_)
Errore: 27.48719662493046
Score: 0.6845806026375012
Valore del peso di RM: 5.047448941276316
Valore del peso di LSTAT: -0.6152722462943342
Valore del bias: -1.363513196568526
Questa volta sia l’errore che lo score sono decisamente migliorati, quindi anche la regressione sarà più precisa, proprio in virtù del fatto che abbiamo utilizzato una proprietà ulteriore per aiutarla nella stima.
Allora esageriamo, usiamo tutte le proprietà del dataset, anche quelle che influiscono in modo minore sul valore dell’abitazione. Visto che le proprietà sono molte e hanno valori molto diversi tra loro, conviene procedere ad una standardizzazione del dataset, per evitare che la regressione tenga troppo in considerazione le proprietà con valori molto grandi rispetto a quelle con valori piccoli, col rischio di virare in modo anomalo le correlazioni reali con MEDV. Ci viene ancora in aiuto SciKitLearn con una classe specifica.
from sklearn.preprocessing import StandardScaler
# istanzio la classe di standardizzazione
# standardizzo il train set, creando un modello di standardizzazione in base ai suoi dati
# riutilizzo lo stesso modello di standardizzazione sul test set, in modo da mantenere uniformità tra i due set
ss = StandardScaler()
X_train_std = ss.fit_transform(X_train)
X_test_std = ss.transform(X_test)
Poi ripeto la regressione e visualizzo le metriche, i pesi ed i bias.
Errore: 16.4539169335068
Score: 0.7861973057489404
Valore dei pesi: [('CRIM', -0.9854756667925126), ('ZN', 1.0024283991840004),
('INDUS', 0.7461760403352479), ('CHAS', 0.8420659269345607),
('NOX', -2.5613274503324277), ('RM', 2.5692913135337907),
('AGE', 0.2667217628061168), ('DIS', -2.9071451015959564),
('RAD', 3.035017456004171), ('TAX', -2.3452764870673537),
('PRATIO', -2.013942132020687), ('B', 0.7833809834666242),
('LSTAT', -3.9424253371618834)]
Valore del bias: 22.449435028248587
Come si nota, la qualità della regressione è aumentata ancora, perchè ora ha a disposizione maggiori informazioni per ottimizzare il risultato. Inoltre, avendo standardizzato il dataset, le varie proprietà hanno ora valori equilibrati tra loro, quindi in questo caso i pesi ci dicono anche quali sono le proprietà che la regressione ha trovato essere più interessanti. A valori assoluti più grandi, corrisponde maggiore correlazione trovata. In questo caso i pesi con valore assoluto >2,5 corrispondono alle proprietà NOX, RM, DIS, RAD, LSTAT. Queste proprietà “pesano” di più nel calcolo del risultato.
See you soon!
Bene, nel prossimo articolo vedremo che non esiste solo la regressione lineare, ma anche la regressione polinomiale, utile quando i valori si diffondono con un andamento che non segue una retta ma bensì una curva (indizio…).
Fate sperimentazione sul codice sorgente che trovate nel gitlab interno. E noi ci rivediamo la prossima volta.
Bye!