Dataa putkeen, käärmettä pyssyyn?

R vastaan Python vol. 2

 

Kirjoitin alkuvuodesta pohdiskelua R:n ja Python:in eroista. Tuolloin vertailin pintapuolisesti molempien kielien datan muokkaustoiminnallisuutta. Molemmat ovat tässä tehtävässä todella hyviä. Vaikka R:ssä on mielestäni viety datan muokkaus, käsittely ja visualisointi toiminnallisesti selvästi pidemmälle tidyverse kokonaisuudella, väitän että Python:issa datan puljaus on laskennallisesti tehokkaampaa. Lisäksi Python:issa on mahdollista hyvin helposti käsitellä isoja datoja iteratiivisesti, mikä ei R:ssä oikein luonnistu.

Ajattelin nyt palata tähän kielikysymykseen vähän syvällisemmin, vertailemalla miten helposti/hyvin näillä kielillä saa toteutettua “automatisoidun” datan esikäsittely-/mallinnusputken. Tällä putkella tarkoitan modulaarista objektia, joka voidaan parametrisoida datalla ja soveltaa tarvittaessa myöhemmin testauksessa/tuotannossa. R:ssä tähän tarkoitukseen on olemassa recipes niminen paketti. Python:issa taas homma pitäisi hoitua sklearn kirjaston pipeline moduulilla.

Tämä kirjoitus, koodeineen ja tuloksineen päivineen on toteutettu R-Studio:ssa, missä on helposti mahdollista käyttää sekä R:ää että Python:ia samassa istunnossa. Tässä kohtaa lukija tyypillisesti vetää herneen nenään, kun helpoksi kehuttua toiminnallisuutta ei avata sen enempää. Siksipä kerronkin tässä heti alkuun miten homma toimii. Olettaen että koneelta löytyy Python/Anaconda asennus, R-Studio:ssa avataan uusi R Notebook, jonka alkuun luodaan ns. setup-blokki; lisätään vain tavallisen {r} blokkiin sana setup, {r setup}. Tämän blokin sisään lisätään seuraavat asiat:

 

# Kutsutaan Pythonin integroinnin mahdollistava "reticulate" paketti:
library(reticulate)

# Lisäksi joku näistä:
use_python("<polku Python asennukseen>")
use_virtualenv("<Python virtuaaliympäristön nimi>")
use_condaenv("<Anaconda-ympäristön nimi>")

 

Muuta ei sitten tarvitakaan (toki reticulate paketti pitää muistaa ensin asentaa).

 

Esimerkkidata

Käytän tässä esimerkkinä klassista koneoppimisdataa Pima intiaanien diabeteksestä. Ainestossa on 21 vuotta täyttäneiden naisten lääketieteellisiä tietoja, mukaan lukien tieto siitä onko henkilöllä todettu diabetes vai ei. Äkkivilkaisulla data näyttää oikein hyvältä; puuttuvia tietoja ei näyttäisi olevan ollenkaan, mikä on itse asiassa aika harvinaista tällaiselle lääketieteelliselle havaintoaineistolle (nimim. kokemusta on). Mutta kuinka ollakaan, datassa puuttuvat arvot on “kätevästi” koodattu nollalla (mikä kuitenkaan tuskin pätee raskauksien määrään). Välttääkseni näistä nollista mahdollisesti aiheutuvia ongelmia olenkin toisaalla korvannut ne NA:lla, mikä tulkitaan puuttuvaksi tiedoksi sekä R:ssä että Python:issa.

 


R-versio: recipes

 

Datan luku

Sitten päästäänkin itse asiaan. Aloitetaan lukemalla datan R:ään ja tulostamalla sen 10 ensimmäistä riviä. Näemme että aineistossa on 8 selittävää muuttujaa vastemuuttujan (Outcome) lisäksi. Valitettavasti erityisesti insuliinitasoissa on paljon puuttuvia arvoja, lähes 50%. Todellisuudessa tällaisen muuttujan käyttäminen mallinnuksessa ei ole kovin järkevää, tai vähintäänkin pitäisi sovittaa eri malli riippuen siitä onko tieto olemassa vai ei.

library(tidyverse)
library(readr)
df = read_csv('diabetes.csv')
head(df,10)
## # A tibble: 10 x 9
##    Pregnancies Glucose BloodPressure SkinThickness Insulin   BMI
##          <int>   <int>         <int>         <int>   <int> <dbl>
##  1           6     148            72            35      NA  33.6
##  2           1      85            66            29      NA  26.6
##  3           8     183            64            NA      NA  23.3
##  4           1      89            66            23      94  28.1
##  5           0     137            40            35     168  43.1
##  6           5     116            74            NA      NA  25.6
##  7           3      78            50            32      88  31  
##  8          10     115            NA            NA      NA  35.3
##  9           2     197            70            45     543  30.5
## 10           8     125            96            NA      NA  NA  
## # … with 3 more variables: DiabetesPedigreeFunction <dbl>, Age <int>,
## #   Outcome <int>

 

Tässä harjoituksessa tulen ihan tarkoituksella sisällyttämään datan muokkausputkeen imputointi-stepin (imputoinnilla tarkoitetaan puuttuvien havaintojen paikkaamista joillakin korvaavilla arvoilla), joten en tee puuttuvien tietojen suurestakaan määrästä ongelmaa.

Koska lopullinen tavoite on tehdä datasta ennustemalli, pitää data vielä pilkkoa opetus- ja testauspaloihin. Tähän tarkoitukseen käytän caret pakettia, jonka createDataPartition() funktio pyrkii pitämään huolen siitä että datan satunnaisotanta on edustava vastemuuttujan suhteen:

 


library(caret)

# Pilkotaan data, 80% käytetään opetukseen:
set.seed(5123)
train = createDataPartition(df$Outcome, p = 0.8)$Resample1

print(paste0('Positiivisten osuus -- koko datassa: ',
             round(table(df$Outcome)[2]/nrow(df),2),
             ', opetus-datassa: ',
             round(table(df[train,'Outcome'])[2]/length(train),2)))
## [1] "Positiivisten osuus -- koko datassa: 0.35, opetus-datassa: 0.35"

 

Reseptiikka

Seuraavaksi muodostan recipes paketilla “reseptin”, jonka mukaan dataa sitten muokataan. Homma alkaa sillä, että recipe() funktiolle kerrotaan miltä data näyttää ja mikä on kunkin muuttujan tulkinta. Koodissa Outcome ~ . formula tarkoittaa, että Outcome käsitellään vasteena ja kaikki muut selittäjinä. Tähän reseptiin voi sitten liittää erinäisen määrän valmistusvaiheita; erilaisia steppejä (step_x()) on tarjolla pyöreästi 70 kappaletta.

Seuraavassa käytän reseptissä neljää vaihetta, joista ensimmäisessä pilkotaan ikä kolmeen luokkaan, toisessa lasketaan dataan uusi muuttuja, kolmannessa täydennetään puuttuvat arvot k-NN menetelmällä (periaate on se, että otetaan kullekin datapisteelle lähimpien naapureiden keskiarvo puuttuvien arvojen tilalle) ja viimeisessä neljännessä vaiheessa skaalataan numeeriset muuttujat nollan ja ykkösen välille:

 


library(recipes)

recip  = recipe(Outcome ~ .,data = df) %>% 
  # Diskretisoidaan ikä:
    step_discretize(Age,options = list(cuts = 3,prefix='AgeClass_')) %>%
  # Muodostetaan uusi muuttuja:
    step_mutate(gluBMI = Glucose / BMI, role = 'predictor') %>%
  # Impoutoidsaan puuttuvat tiedot k-NN menetelmällä:
    step_knnimpute(all_predictors(), neighbors = 10) %>%
  # Skaalataan numeeriset muuttujat välille [0, 1]:
    step_range(all_numeric())

# Printataan resepti
recip
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          8
## 
## Operations:
## 
## Dummy variables from Age
## Variable mutation for  gluBMI
## 10-nearest neighbor imputation for all_predictors()
## Range scaling to [0,1] for all_numeric()

 

Tästä printistä nähdään reseptin tulkinta datan muuttujille, sekä jokaisen “työvaiheen” sanallinen kuvaus. En ole valinnut kyseisiä työvaiheita sen perusteella että nämä olisivat jotenkin optimaalisia datan mallinnuksen kannalta vaan puhtaasti demonstratiivisessa mielessä.

Kun resepti on luotu, se täytyy “kouluttaa”, mikä tarkoittaa kunkin työvaiheen parametrien kalibrointia datan perusteella. Tähän käytetään prep() funktiota, joka valmistelee reseptin käyttäen opetus-dataa. “Treenattusta” reseptistä voidaan sitten puristaa ulos reseptin mukaan käsitelty opetusdata kuvaavasti juice() funktiolla. Vastaavasti saamme reseptillä leivottua testidatan bake() funktiolla. Oleellista tässä on se, että testidata muodostetaan samoilla parametreillä kuin opetusdata. Samalla periaatteella voidaan tuotannossa varmistaa, että data käsitellään aina taatusti samalla tavalla.

 

# Treenataan resepti:
recip_prep = prep(recip, training = df[train,])

# Puristetaan opetusdata ulos:
df_tr = juice(recip_prep)

# Leivotaan testidata:
df_ts = bake(recip_prep, new_data = df[-train,])

 

Ennustaminen

Mielestäni R:n resepti-putki on loistavasti toteutettu. Siinä on kuitenkin aika selvä puute; putkeen ei voi liittää ennustemallin koulutusta (tai se on vähintäänkin melko työlästä). Minulla on tässä tavoitteena optimoida mallin parametrejä riistiinvalidointia hyödyntäen, mikä onnistuu näppärästi caret paketin (joka on muuten todella hyvin dokumentoitu) train() funktiolla.

Käytän mallinnukseen C5.0 mallipuu-menetelmää ja optimoin mallin hyperparametrejä, käyttäen viisinkertaista ristiinvalidointia:

 

require(caret)

# Opetusdata:
X_tr = select(df_tr,-Outcome) %>% as.data.frame()
Y_tr = factor(df_tr$Outcome)

# Sovituksen kontrollointi:
Tcont = trainControl(method='cv', number = 5)
Tgrid = expand.grid(winnow = TRUE, 
                    model = 'tree', 
                    trials = 1:3)

# Mallin sovitus:
model = caret::train(x = X_tr,y = Y_tr,
                     method = 'C5.0',trControl = Tcont,
                     tuneGrid = Tgrid, metric = 'Kappa')

# Tulostetaan malli:
model
## C5.0 
## 
## 615 samples
##   9 predictor
##   2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 491, 492, 493, 492, 492 
## Resampling results across tuning parameters:
## 
##   trials  Accuracy   Kappa    
##   1       0.7171056  0.3645696
##   2       0.7251308  0.3791710
##   3       0.7251439  0.3892372
## 
## Tuning parameter 'model' was held constant at a value of tree
## 
## Tuning parameter 'winnow' was held constant at a value of TRUE
## Kappa was used to select the optimal model using the largest value.
## The final values used for the model were trials = 3, model = tree
##  and winnow = TRUE.

 

Valmista tuli! Lopullisen mallin ennustetarkkuus on (ristiinvalidoituna) pyöreästi 0.73, joka lasketaan vertaamalla ao. konfuusiomatriisin oikeaan osuneiden ennusteiden määrää kaikkiin tapauksiin:

 

 

Parametreistä sen verran että winnow (joka on tai ei ole käytössä) tarkoittaa muuttujien valikointia ja trials tarkoittaa buustaus iteraatioiden määrää (model: puu- malli vai sääntömalli?). Valitsin arviointikriteeriksi Kappa:n, joka on käytännössä korjattu ennustetarkkuus (ottaa huomioon kuinka todennäköisesti saadaan positiivinen lopputulos sattumalta).

Viimeinen homma on laskea eri metriikoita mallista, kun lasketaan mallin ennusteet testidatalle. Tässä kohtaa jätän koodin esittämättä, ihan vain tilaa säästääkseni:

 

Tarkkuus Precision AUC
opetus 0.83 0.75 0.86
testi 0.78 0.71 0.85

 

Tästä taulukosta nähdään, että tarkkuus (eli accuracy) on vähän parempi opetus- kuin testidatalle. Sama pätee mallin precision-arvoon, eli mallin kyky ennustaa positiiviset tapaukset oikein. Tämä arvo kertoo usein paljon paremmin mallin hyvyyden accuracy:yn verrattuna, erityisesti silloin kun positiiviset tapaukset ovat paljon harvinaisempia kuin negatiiviset. Lopuksi, AUC -arvo, joka kertoo mallin kyvystä erotella positiiviset ja negatiiviset tapaukset. Saatu malli on hieman ylioptimistinen, eli opetusdatan antamat metriikat ovat suurempia kuin testidatan. Tässä esimerkkitapauksessa nämä luvut vaihtelevat jonkin verran satunnaislukujen vaikutuksesta, koska dataa on melko vähän.

 


Python-versio: Pipeline

 

Datan luku

Aloitetaan taas lukemalla data, tällä kertaa Python:in ymmärtämässä muodossa (tosin voisin myös käyttää r. notaatiota olemassa olevan df objektin edessä, jolloin Python osaa hyödyntää R:n lukeman datan):

 

import pandas as pd
df = pd.read_csv('diabetes.csv')
print(df.info())
## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 768 entries, 0 to 767
## Data columns (total 9 columns):
## Pregnancies                 768 non-null int64
## Glucose                     763 non-null float64
## BloodPressure               733 non-null float64
## SkinThickness               541 non-null float64
## Insulin                     394 non-null float64
## BMI                         757 non-null float64
## DiabetesPedigreeFunction    768 non-null float64
## Age                         768 non-null int64
## Outcome                     768 non-null int64
## dtypes: float64(6), int64(3)
## memory usage: 54.1 KB
## None

 

Datan esikäsittely

Python:issa caret ja recipe pakettien toiminnallisuus löytyy pitkälti sklearn kirjastosta. Selkeänä rajoitteena on kuitenkin se, että sklearn operoi numpy array datarakenteella, mikä vastaa R:n matriiseja. Näin ollen datan prosessointiputkeen ei voida sisällyttää datan rakenteeseen tapahtuvia muutoksia, vaan ne täytyy hoitaa erikseen pandas kirjastolla:

 

# Luodaan uusi muuttuja kuten yllä:
df['gluBMI'] = df['Glucose'] / df['BMI']
# Diskretisoidaan ikä, kuten yllä:
df['Age'] = pd.cut(df.Age,bins = 3, labels = ['young','middle-aged','old'])
df = pd.get_dummies(df)

 

Datan palastelu

Kuten R:ssä data pilkotaan kahtia. Python:issa tämä hoituu esim. train_test_split() funktiolla, joka myöskin yrittää pitää huolen siitä että satunnaisotanta on edustava.

 

from sklearn.model_selection import train_test_split
# Muutetaan data array-muotoon:
X = df.drop(['Outcome'],axis = 1).values
Y = df['Outcome'].values

# Ositetaan data, 80% opetukseen:
X_tr, X_ts, Y_tr, Y_ts = train_test_split(X,Y, 
                         test_size = 0.2, random_state = 5123)
print('Positiivisten osuus -- koko data: %.2f, opetus-datassa: %.2f' %(Y.mean(), Y_tr.mean()))
## Positiivisten osuus -- koko data: 0.35, opetus-datassa: 0.34

 

Putkityöt

Sitten päästäänkin vihdoin LVI-hommiin, tai siis koostamaan Python data/mallinnusputkea. Tosiaankin, toisin kuin R:ssä, sklearn mahdollistaa mallin kouluttamisen osana putkea, mikä onkin todella näppärää. Valitettavasti k-NN-imputointi ei onnistu yhtä näppärästi kuin R:ssä, joten käytän tässä SimpleImputer() menetelmää, jonka parametrit voidaan optimoida samalla kun putki treenataan. Kuten yllä, MinMaxScaler() skaalaa muuttujat [0,1]-välille ja putken päässä odottaa KNeighborsClassifier(), joka sovittaa dataan ns. lähinaapuriluokittelijan (sklearn kirjasto ei sisällä C5.0 algoritmiä, kun taas k-NN malli löytyy helpommin kuin R:ssä):

 

from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import MinMaxScaler
from sklearn.neighbors import KNeighborsClassifier

# Kootaan tarvittavat välivaiheeta:
steps = [('Imputointi', SimpleImputer()),
         ('Skaalaus', MinMaxScaler()),
         ('Luokittelija',KNeighborsClassifier())]

# Muodostetaan putki välivaiheista:
pipe = Pipeline(steps)

 

Printattuna putken vaiheet näyttävät tältä, mistä näemme kunkin funktion oletusparametrit:

 

print(pipe.steps)
## [('Imputointi', SimpleImputer(copy=True, fill_value=None, missing_values=nan, strategy='mean',
##        verbose=0)), ('Skaalaus', MinMaxScaler(copy=True, feature_range=(0, 1))), ('Luokittelija', KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
##            metric_params=None, n_jobs=None, n_neighbors=5, p=2,
##            weights='uniform'))]

 

Ennustaminen

Python:issa voidaan myös kätevästi syöttää putkitus parametrien optimointialgoritmiin. Seuraavassa paketoin putken säännöllisellä parametrien haarukoinnilla, perustuen viisinkertaiseen ristiinvalidointiin. Lopuksi koko roska treenataan opetusdatalla:

 

import numpy as np
from sklearn.model_selection import GridSearchCV
parameters = {'Imputointi__strategy':['mean','median'],
              'Luokittelija__weights':['uniform','distance'],
              'Luokittelija__p':[1,2],
              'Luokittelija__n_neighbors':np.arange(3,11)}

# Alustetaan GridSearchCV objekti:
cv = GridSearchCV(pipe,parameters,cv=5)

# Treenataan opetus-datalla:
cv.fit(X_tr,Y_tr)

 

Haarukoinnin löytämät “optimaaliset” parametrit:

print(cv.best_params_)
## {'Imputointi__strategy': 'median', 'Luokittelija__n_neighbors': 5, 'Luokittelija__p': 2, 
'Luokittelija__weights': 'uniform'}

 

eli imputointi tehtiin keskiarvokorvauksella ja luokittelussa käytettiin viittä lähintä naapuria luokan määräämiseen, perustuen euklidiseen etäisyyteen (p = 2), eikä etäisyyksiä painotettu mitenkään (‘uniform’).

Lopuksi arvioidaan mallin hyvyyttä samoilla metriikoilla kuin R:ssä. Katsotaan ensin mallin konfuusiomatriisia opetusdatalle:

 

 

Näyttäisi siltä että ainakin väärien positiivisten määrä on tässä mallissa pienempi (taaskin satunnaisluvut voivat vaikuttaa tähän tulokseen).

Metriikka lasketaan hyvin samalla tavalla kuin R:ssä. Alla käytän ns. list ja dict comprehension -rakenteita, jotka vastaavat hyvin pitkälti R:n apply() ja lapply() funktioita:

 

from sklearn.metrics import roc_auc_score, precision_score
splits = {'train':X_tr, 'test':X_ts}
ys = {'train':Y_tr, 'test':Y_ts}

preds = {x:cv.predict(splits[x]) for x in splits.keys()}
predp = {x:cv.predict_proba(splits[x])[:,1] for x in splits.keys()}

out = pd.DataFrame({
 'Tarkkuus':[cv.score(X_tr,Y_tr),cv.score(X_ts,Y_ts)],
 'Precision':[precision_score(ys[ii],preds[ii]) for ii in ys.keys()],
 'AUC':[roc_auc_score(ys[ii],predp[ii]) for ii in ys.keys()]
 })

out.index = ['opetus','testi']

 

Taulokoituna samoin kuin yllä:

 

Tarkkuus Precision AUC
opetus 0.82 0.75 0.90
testi 0.78 0.72 0.81

 

mistä nähdään että malli on hieman ylioptimistinen, koska testidatan metriikka poikkeaa opetusdatan metriikasta (toki on ihan tavallista, että testidatalla saadaan hieman huonommat arvot, koska ko. aineistoa ei ole käytetty mallin opettamiseen ja data oli kokonaisuutena varsin pieni). Varsinkin AUC -arvo on pienempi sekä opetusdataan, että R:n malliversioon verrattuna. Mutta kuten sanottua, reilu vertailu tulisi suorittaa niin, että myös opetus- ja testidatoja arvottaisiin useampi eri setti. Näin päästäisiin kiinni keskimääräiseen performanssiin.

 


Lopputulema

Jaa että kummanko kielen valitsisin tämän harjoituksen perusteella? Molemmissa on puolensa. R:ssä datan esikäsittelyputki on mielestäni paremmin toteutettu (tähän löytyy myös laajennuksi, esim. tekstianalytiikkaan), mutta Python:issa putkeen saa liitettyä myös ennustemallin. R:ssä on paremmin tarjolla erilaisia algoritmejä, mutta Python:issa suurin osa tarpeellisista menetelmistä löytyy samasta paikasta ja ne toimivat moitteettomasti. Mallien (mukaan lukien imputointi, skaalaus, enkoodaus ym. mallien) hyperparametrien optimointi on myös mielestäni Python:issa toteutettu paremmin kuin R:ssä. Toisaalta jos jaksaa itse koodata tämän iteroinnin (mikä ei oikeastaan ole kovin vaikeaa), saa tietty aina mieleisen kielestä riippumatta.

Hmmm… vaikea paikka… Mutta hetkinen, minähän tein koko jutun R-Studio:ssa, missä voin käyttää molempia kieliä rinnan; ei tarvitsekaan valita :D. Näin on mahdollista yhdistää molempien kielien vahvuudet. Rinnakkaista käyttöä helpottaa merkittävästi se, että R-objekteja voi Python blokeissa (r. etuliitteellä) ja Python-objekteja vastaavasti R-blokeissa (py$ etuliitteellä), mikä ei ole mitenkään itsestäänselvää. Eli, ota(n) molemmat. Hyvää koodausta kaikille!

Koulutuksia

Python for Data Science

R-ohjelma tutuksi