Mallin ennustuskyvyn arviointi: Ristiinvalidointi

Kun teemme mallin, jolla haluamme ennustaa, niin on tärkeää arvioida mallin ennustuskykyä etukäteen. Tämä johtuu siitä, että mallit yleensä ylisovittuvat dataan, jolla ne luodaan, eli niiden ennustuskyky ei olekaan enää hyvä uudella datalla. Tämän ongelman ratkaisuksi on keksitty erilaisia ristiinvalidointimenetelmiä, joista esittelen K-ositus -ristiinvalidoinnin (K-fold cross-validation). Menetelmän avulla saamme mallin ennustuskyvystä vakaamman kuvan kuin, jos käyttäisimme perinteistä trainingset ja testset -jakoa. Ristiinvalidointi on erityisen tärkeää, kun datan määrä on pieni, koska tällöin saamme äärimmäisiä tuloksia todennäköisemmin.

K-ositus ristiinvalidoinnin vaiheet:

1. Jaa data k:hon samansuuruiseen osaan. Esim. kun k=3, niin data jaetaan kolmeen samansuureiseen osaan.
2. Muodosta malli käyttäen k-1 osia (trainingset) ja ennusta sillä puuttuvan osan (testset) arvot. Toista tämä, kunnes jokainen osa on ollut kerran testsettinä.
3. Laske kunkin mallin testsetille tuottamien AUC-arvojen keskiarvo.


Esimerkki, kun k=3 ja osat ovat osa1, osa2 ja osa3. Muodostetaan ensin mallit:

Malli 1: Trainingset: osa1 ja osa2. Testset: osa3.
Malli 2: Trainingset: osa1 ja osa3. Testset: osa2.
Malli 3: Trainingset: osa2 ja osa3. Testset: Osa1.

Seuraavaksi malleilla ennustetaan testsetteihin ja saaduista AUC-arvoista lasketaan keskiarvot. Jos saamme esim. AUC-arvot 0.82, 0.78 ja 0.76, niin keskiarvo on ~0.787. Ristiinvalidoinnin tärkeys tulee ilmi tässä kohtaa, koska ilman sitä olisimme voineet päätyä Malli 1:n antamaan AUC-arvoon (0.82), joka antaa liian positiivisen kuvan mallin ennustustarkkuudesta. Ristiinvalidointi antaa siis vakaamman kuvan mallin ennustustarkkuudesta. Esimerkissä k oli kolme, mutta se voi olla myös esim. kymmenen. Lisäämällä k:n arvoa voimme pienentää ristiinvalidoinnista saamaamme AUC-arvon varianssia. Varianssia saadaan myös pienemmäksi, kun ristiinvalidointi suoritetaan useampaan kertaan ja kunkin ristiinvalidoinnin saadut AUC-arvot keskiarvoistetaan. Voimme esim. suorittaa k=10 ristiinvalidoinnin viisi kertaa, jolloin saamamme lopputulos olisi viiden ristiinvalidoinnin tuottaman AUC-arvon keskiarvo.

Tarkastellaan seuraavaksi aiemmin käyttämäämme luottadataa (data) ristiinvalidoinnin näkökulmasta. Lasken mallille ensin AUC-arvon, kun käytämme 60% trainingset ja 40% testset jakoa.


data <- read.delim(file.choose(), header=FALSE,sep="")
names(data) <- c("check_acc", "duration", "history", "purpose", "amount", "sav_acc", "employment", "install_rate", "ptatus", "other_debtor", "resid", "property","age", "other_install", "housing", "other_credits", "job", "num_depend", "teleph", "foreign", "default")
data$install_rate <- as.factor(data$install_rate)
data$resid <- as.factor(data$resid)
data$other_credits <- as.factor(data$other_credits)
data$num_depend <- as.factor(data$num_depend)
data$default <- as.factor(ifelse(data$default==1,"0","1"))

trainingdata$age <- as.factor(ifelse(trainingdata$age <= 30,"0-30",ifelse(trainingdata$age <=40,"31-40","40+")))
trainingdata$amount <-as.factor(ifelse(trainingdata$amount<=2500,"2500",ifelse(trainingdata$amount<=5000,'2600-5000','5000+')))
trainingdata$duration <- as.factor(ifelse(trainingdata$duration <= 12,"1-12",ifelse(trainingdata$duration <=24,"13-24","24+")))

d <- sort(sample(nrow(data),0.6*nrow(data)))
trainingdata <- data[d,]
testdata <- data[-d,]

model1 <- glm(default~.,family=binomial(link=logit),data=trainingdata)
testdata$log <- predict(model1, testdata, type = "response")
f<- prediction(testdata$log,testdata$default)
performance(f,"auc") # AUC-arvo: 0.791

library(ggplot2)
library(plotROC)
plot <- ggplot(testdata,aes(d=default,m=log))+geom_roc()+style_roc()
plot

Saimme AUC-arvoksi 0.791, mutta onko tämä AUC-arvo sellainen, jonka voimme odottaa saavamme myös jatkossa? Tietysti on mahdollista saada kyseinen arvo myös jatkossa, mutta ristiinvalidoinnin avulla saamme parempaa tietoa AUC-arvon odotusarvosta, eli voimme karsia epävarmuutta. Seuraavaksi teen ristiinvalidoinnin, kun k=5.
R-koodi:

library(caret)
d <- createFolds(data$default, k = 5, list = TRUE, returnTrain = FALSE) # Jaetaan data viiteen osaan
names(d)[c(1,2,3,4,5)] <- c("t1","t2","t3","t4","t5") 

m1 <- glm(default~., family=binomial, data=data[c(d$t1,d$t2,d$t3,d$t4),])
g <- predict(m1,data[d$t5,],type="response")
pre1 <- prediction(g,data[d$t5,"default"])
a1 <- performance(pre1,"auc") # AUC: 0.824

m2 <- glm(default~., family=binomial, data[c(d$t1,d$t2,d$t3,d$t5),])
g <- predict(m2,data[d$t4,],type="response")
pre2 <- prediction(g,data[d$t4,"default"])
b1 <- performance(pre2,"auc") # AUC: 0.721

m3 <- glm(default~., family=binomial, data[c(d$t1,d$t2,d$t4,d$t5),])
g <- predict(m3,data[d$t3,],type="response")
pre3 <- prediction(g,data[d$t3,"default"])
c1 <- performance(pre3,"auc") # AUC: 0.801

m4 <- glm(default~., family=binomial, data[c(d$t1,d$t3,d$t4,d$t5),])
g <- predict(m4,data[d$t2,],type="response")
pre4 <- prediction(g,data[d$t2,"default"])
d1 <- performance(pre4,"auc") # AUC: 0.808


m5 <- glm(default~., family=binomial, data[c(d$t2,d$t3,d$t4,d$t5),])
g <- predict(m5,data[d$t1,],type="response")
pre5 <- prediction(g,data[d$t1,"default"])
e1 <- performance(pre5,"auc") # AUC: 0.733


(as.numeric(a1@y.values)+as.numeric(b1@y.values) # Lasketaan keskiarvo
+as.numeric(c1@y.values)+as.numeric(d1@y.values)+as.numeric(e1@y.values))/5

Keskiarvoksi saadaan 0.778, joten näyttää siltä, että ristiinvalidointi paljasti, ettei malli ollutkaan niin hyvä kuin aluksi vaikutti. Kun katsomme saatuja AUC-arvoja, niin huomaamme, kuinka paljon ne vaihtelevat välillä 0.72-0.80. Tämä vaihtelu johtuu datan pienestä koosta (1000 havaintoa). Toistan seuraavaksi ristiinvalidoinnin viisi kertaa ja lasken saatujen arvojen keskiarvon. Arvojen keskiarvoksi saadaan 0.773 ja arvojen vaihteluväli on 0.764-0.777, eli arvojen hajonta on huomattavasti pienempi. Voimme siis päätellä, että mallin ennustukset uudella datalla tuottavat noin 0.77 suuruisen AUC-arvon.

Machine learning: Blending

Esittelen nyt tehokkaan koneoppimisen menetelmän, jolla mallien ennusteita voidaan parantaan. Menetelmää kutsutaan nimellä "Blending". Ideana on yhdistää usean mallin ennusteet ja näin saavuttaa parempi ennustustarkkuus. Voimme siis yhdistää esim. logistisesta regressiosta ja random forestista saamamme ennusteet ja näin saavuttaa paremmat ennusteet. Blending toimii parhaiten, jos eri mallit pystyvät edes kohtalaisiin ennusteisiin (AUC>0.5) ja niiden välinen korrelaatio ei ole suuri. Jos yhdistämme kaksi mallia, jotka korreloivat vahvasti, niin emme saa uutta tietoa, vaan uusi ennuste vastaa lähes kokonaan yksittäisen mallin ennustetta.  Mikäli mallit eivät korreloi vahvasti, niin parhaimmillaan yhden mallin oikeat ennusteet dominoivat, kun toinen malli on väärässä ja näin saadaan aikaan parempia ennusteita. Blending ei aina paranna mallien ennusteita, mutta se ei kuitenkaan yleensä huononna ennusteita verrattuna yksittäisiin malleihin, joten sitä kannattaa kokeilla. Lisätietoa löytyy mm. täältä

Esimerkki, kun blending parantaa ennustustarkkuutta: Logistisen mallin ennusteet ovat kolmen lainanhakijan maksuhäiriöiden todennäköisyyksistä seuraavat:

Hakija1: 0.09, Hakija2: 0.05, Hakija3: 0.045.

Vastaavat random forestista saadut todennäköisyydet ovat:

Hakija1: 0.08, Hakija2: 0.09,  Hakija3: 0.1

Kuvitellaan, että todellisuudessa Hakijat 1 ja 3 saivat maksuhäiriön. Logistinen malli ennusti siis hakijat maksuhäiriön todennäköisyyden mukaisesti suuruusjärjestykseen seuraavasti: Hakija1, Hakija2, Hakija3. Random forestilla vastaava järjestys oli: Hakija3, Hakija2, Hakija1. Oikea järjestys olisi ollut Hakija1, Hakija3, Hakija2, eli hakijoiden 1 ja 3 olisi pitänyt saada suuremmat maksuhäiriön todennäköisyydet kuin toisen hakijan. Yhdistämällä mallit seuraavasti saamme optimaalisen mallin: 0.5* Logistisen mallin ennusteet + 0.5* Random forestin ennusteet. Nyt hakijoiden maksuhäiriöiden todennäköisyydet ovat: Hakija1: 0.085, Hakija2: 0.07, Hakija3: 0.0725. Saimme yhdistämällä mallin, joka on parempi kuin yksittäiset mallit.

Kokeilen seuraavaksi kolmen eri mallin yhdistämistä. Muodostan mallit jo tutusta luottodatasta (GermanCreditData). Malleina käytän logistista regressiota, random forestia ja support vector machinea.

R-koodi:

data <- read.delim(file.choose(), header=FALSE,sep="") # Avataan data ja nimetään muuttujat
names(data) <- c("check_acc", "duration", "history", "purpose", "amount", "sav_acc", "employment", "install_rate", "ptatus", "other_debtor", "resid", "property","age", "other_install", "housing", "other_credits", "job", "num_depend", "teleph", "foreign", "default")

data$install_rate <- as.factor(data$install_rate)

data$resid <- as.factor(data$resid)

data$other_credits <- as.factor(data$other_credits)

data$num_depend <- as.factor(data$num_depend)

data$default <- as.factor(ifelse(data$default==1,"0","1"))
data$amount <- as.factor(ifelse(trainingdata$amount<=2500,"2500",ifelse(data$amount<=5000,'2600-5000','5000+')))
data$duration <- as.factor(ifelse(trainingdata$duration <= 12,"1-12",ifelse(datadata$duration <=24,"13-24","24+")))

data$age <- as.factor(ifelse(testdata$age <= 30,"0-30",ifelse(data$age <=40,"31-40","40+")))
d <- sort(sample(nrow(data),0.6*nrow(data)))  # Muodostetaan training- ja testdata

trainingdata <- data[d,]

testdata <- data[-d,]

model1 <- glm(default~check_acc+duration+history+purpose+  # Tehdään logistinen malli

amount+

sav_acc+

employment+

install_rate+

ptatus+

other_debtor+

resid+

property+

age+

other_install+

housing+

other_credits+

job+

num_depend+

teleph+

foreign,family=binomial(link=logit),data=trainingdata) 

testdata$log <- predict(model1, testdata,type="response") # Muodostetaan mallin ennusteet

m1.scores <- prediction(testdata$log, testdata$default)

performance(m1.scores,"auc") # AUC 0.773

library(randomForest)  # Muodostetaan random forest -malli

model2 <- randomForest(default~check_acc+duration+history+purpose+

amount+

sav_acc+

employment+

install_rate+

ptatus+

other_debtor+

resid+

property+

age+

other_install+

housing+

other_credits+

job+

num_depend+

teleph+

foreign, ntree=1000, data=trainingdata)

testdata$rf <- predict(model2, testdata,type="prob")[,2]  # Muodostetaan mallin ennusteet

m2.scores <- prediction(testdata$rf, testdata$default)

performance(m2.scores,"auc") # AUC 0.794

library(kernlab)

model3 <- ksvm(default~check_acc+duration+history+purpose+ # SVM-malli

amount+

sav_acc+

employment+

install_rate+

ptatus+

other_debtor+

resid+

property+

age+

other_install+

housing+

other_credits+

job+

num_depend+

teleph+

foreign, data=trainingdata,

cross=5,type='eps-svr' )

svm <- predict(mode13, testdata,type="response") # Muodostetaan mallin ennusteet

logsvm <- glm(testdata$default~svm,family=binomial) #Huom. SVM ei anna todennäköisyyksiä suoraan, joten ne ratkaistaan erikseen logistisella mallilla.

testdata$svm <- predict(logsvm, testdata,type="response")

m3.scores <- prediction(testdata$svm, testdata$default)

performance(m3.scores,"auc") # AUC 0.778

Random forest ennusti parhaiten testidataan (AUC: 0.794). Kokeillaan parantaako mallien keskiarvottaminen ennusteita.

R-koodi:

preyhd <- (testdata$m1.log+testdata$rf+testdata$svm)/3
m4.scores <- prediction(preyhd, testdata$default)
performance(m4.scores,"auc") # AUC 0.796

Mallien yhdistäminen paransi hieman AUC-arvoa. Tehtävänä on seuraavaksi löytää optimaaliset painot ennusteille, eli voimme esim. painottaa random forestin ennusteita enemmän. Painojen löytämiseen voi käyttää esim. meta-luokittelijaa, jolloin puhutaan ns. stacking:sta. Tällöin käytämme esim. logistista regressiomallia apuna oikeiden painojen löytämiseen. Idea on laittaa mallien ennusteet mallin selittäviksi tekijöiksi ja muodostaa näiden kertoimista painot.

Esittelen nyt menetelmän, jolla testidatan AUC-arvo voidaan maksimoida. Täytyy tietysti muistaa, ettei tälle testidatalle optimaalisen AUC:n antavat painot välttämättä toimi parhaiten uudelle datalle. Saamme kuitenkin arvokasta tietoa oikeiden painojen suunnasta.


install.packages("Metrics") 
library("Metrics")
cols <- c("m1.log","rf","svm")
fn.opt.pred <- function(pars, data) {
    pars.m <- matrix(rep(pars,each=nrow(data)),nrow=nrow(data))
    rowSums(data*pars.m)
}
fn.opt <- function(pars) {
    -auc(testdata$default, fn.opt.pred(pars, testdata[,cols]))
}
pars <- rep(1/length(cols),length(cols))
opt.result <- optim(pars, fn.opt, control = list(trace = T))

test.pred <- fn.opt.pred(opt.result$par, testdata[,cols])

print(opt.result$par)
# Saadaan optimaaliset painot: 0.3890623 0.8060392 0.1767558

Näillä painoilla saamme AUC-arvoksi 0.799. On hyvä huomata, että random forest saa suurimman painon. Random forest oli alunperin malleista ennusteissaan paras, joten parhaalle mallille on luontevaa antaa suurin paino malleja yhdistellessä.  

Dynaaminen regressiomalli

Kun haluamme sovittaa regressiomallin aikasarjadataan, niin on tärkeää huomioida aikasarjoille tyypillinen epästationaarinen rakenne, eli rakenne, jossa aikasarjan aiempi arvo ennustaa sen tulevaa arvoa. Alla on esimerkki epästationaarisesta aikasarjasta, jossa jälkimmäisen kuukauden arvo on lähes aina aiemman kuukauden arvoa suurempi.

Lähde: http://www.uow.edu.au/student/qualities/statlit/module3/9time/index.html

Jos tällaiseen dataan sijoittaa regressiomallin suoraan, niin sekä regressiokertoimet että keskivirheet ovat vääristyneitä. Tällöin on kyse ns. Spurious Regression:sta, jolloin saisimme pienet p-arvot laittamalla muuttujaksi esim. maapallon keskilämpötilan kehityksen samana ajankohtana. Spurious regression siis altistaa mallin vääränlaisille tulkinnoille. Aikasarjadataa käsitellessä on siis tärkeää tutkia ensin, onko data epästationaarisessa muodossa vai ei. Datan saa stationaariseen muotoon differensoimalla sen tarpeeksi monta kertaa. Yleensä 1-2 kertaa riittää. Jos datan muodostaneen AR(1)-termin kerroin on yli yksi, niin data voi tarvita kolme differensointia. Alla on esimerkki datasta, joka vaatii kolme differensointia:
R-koodi:

z <-rnorm(100,1,0)
g <- rnorm(100,mean=0,sd=1)
for (t in 2:100) z[t]=1.1*z[t-1]+g[t]

Dynaaminen regressiomalli

Käsittelen nyt tilannetta, jossa regressiomallin virhetermi seuraa ARIMA-mallia, eli regression tehokkuusehdot eivät täyty, koska sen virhetermit korreloivat keskenään. Tällöin emme hyödynnä kaikkea informaatiota mallissamme ja lisäksi keskivirheet ovat vääristyneitä. On siis tärkeää hyödyntää ARIMA-mallia yhdessä regressiomallin kanssa. Muodostan ensin datan, jossa y on riippuvainen x:sta sekä ARIMA-mallia seuraavasta virhetermistä.
R-koodi

y  <- rnorm(100,mean=0,sd=1)
g <- rnorm(100,mean=0,sd=1)
x <- rnorm(100,mean=0,sd=5)
e <- rnorm(100,mean=0,sd=1)
for (t in 2:100) e[t]=e[t-1]+g[t] # Muodostan yksikköjuurellisen (seuraa ARIMA-mallia) virhetermin
for (t in 2:100) y[t]=x[t]+e[t]
library(tseries)
adf.test(y) # Dickey-Fuller -testin mukaan y on stationaarinen.

Seuraavaksi muodostan regressiomallin, jolla saamme x:n kertoimeksi 1.040 (Keskivirhe on 0.038). Testaan sitten virhetermien autokorrelaation Breusch-Godfrey -testillä, jotta näen, ovatko mallin tulokset luotettavia. Testin p-arvo on pienempi kuin 0.05, mikä viittaa virhetermeissä esiintyvään autokorrelaatioon. Alla olevasta kuvasta näemme myös selvästi tämän autokorrelaation.
R-koodi:

model1<- lm(y[2:100]~x[2:100])
summary(f)
library(lmtest)
bgtest(f,5) 

Model1:n autokorrelaatiokuvat.

Osittaisautokorrelaation piikki ensimmäisessä viiveessä on tyypillistä AR(1)-prosessille, eli dynaamisen regressiomallin muodostaminen on hyvä aloittaa mallilla, jossa on x sekä AR(1)-termi.
R-koodi:

library(forecast)
data <- data.frame(y[2:100],x[2:100])
c <- ts(data,start=c(1999,1),frequency=12) # Muunnan datan aikasarja-muotoon.
arimamodel1 <- Arima(c[,1],xreg=c[,2],order=c(1,0,0)) # Dynaaminen regressiomalli, jossa mukana x ja AR(1)-termi.
Box.test(resid(arimamodel1),type="Ljung",lag=5,fitdf=2) # testaa mallin residuaalien autokorrelaatiota.
res <-resid(arimamodel1)
hist(res, breaks=10, density=10, col="blue", main="Dynaamisen regressiomallin residuaalit", xlab="Residuaalit",freq=FALSE)  # Tekee histogrammin residuaaleista.
x<-seq(min(res),max(res),length=40) 
y<-dnorm(x,mean=mean(res),sd=sd(res)) 
lines(x, y, col="black", lwd=2)

Box-Ljung -testin arvo 0.877 tarkoittaa, että nollahypoteesi (H0: ei autokorrelaatiota virhetermeissä) jää voimaan. Yllä olevasta histogrammista näemme, että residuaalit ovat approksimatiivisesti normaalijakautuneita. Mallin kerroin x:lle on 1.000 ja keskivirhe 0.013, eli pääsimme erittäin lähelle x:n oikeaa arvoa. Testasin vielä muita ARIMA-malleja, mutta tulokset eivät juuri parantuneet. Alla on vielä kuva dynaamisen mallin ennustamista arvoista ja oikeista arvoista.
R-koodi:

plot(arimamodel1$x,col="blue",ylab="y")
lines(fitted(arimamodel1),col="red")
legend("topright",c("Mallin ennusteet","Oikeat y:n arvot"),col=c("red","blue"),lty=c(1,1))