Yksisuuntainen varianssianalyysi

Varianssianalyysin avulla voimme tutkia, eroavatko eri ryhmien määrällisesti mitattavat ominaisuudet toisistaan. Voimme olla esim. kiinnostuneita siitä, onko miesten ja naisten keskipituuksissa eroja. Varianssianalyysissa on neljä perusoletusta:

1. Ryhmien jäsenet eivät esiinny muissa kuin omassa ryhmässään (ns. Independence-oletus)
2. Ryhmien arvot on poimittu normaalijakaumasta (tosin mallin residuaalien normaalisuus riittää)
3. Ryhmien varianssit ovat yhtäsuuret
4. Otoskoot ovat samansuuruiset.

Tärkein näistä oletuksista on ensimmäinen, koska jos joidenkin ryhmien jäsenet esiintyvät myös muissa ryhmissä, niin tulokset ovat vääristyneitä, eikä tälle ole muuta korjausta kuin palauttaa oletus toteen. Oletukset 2 ja 3 eivät usein täysin toteudu, mutta esim. varianssin yhtäsuuruus -oletukseen on keksitty ratkaisuja (esim.Welch-correction), joilla tuloksista saadaan luotettavia, vaikka oletus ei pätisi. Mikäli arvoja ei ole poimittu normaalijaukamasta, niin voi käyttää "ei-parametrista" testiä, kuten Kruskall-Wallis -testiä, jossa tätä oletusta ei ole. Toisaalta oleellisinta on, että mallin residuaalit ovat normaalisti jakautuneita, joten ennen Kruskall-Wallis -testin käyttöä on hyvä tutkia residuaalit. Viimeinen ehto on helppo toteuttaa, eli jos jostain ryhmästä on enemmän dataa kuin toisesta, niin käyttää tästä ryhmästä vain toisen ryhmän verran dataa, jolloin otoskoot ovat ryhmille samansuuruiset.

Teen seuraavaksi esimerkin sukupuolten keskipituuksien tutkimisesta. Muodostan datan, jossa kummankin sukupuolen pituudet ovat normaalijakautuneita. Ajatellaan, että data "populaatio" kertoo koko populaation pituudet ja että meillä on mahdollisuus ottaa tuosta populaatiosta 400:n havainnon otos. Tiedämme siis naisten ja miesten oikeat keskipituudet, joten voimme tutkita, pääsemmekö oikeisiin tuloksiin otoksen avulla.

R-koodi:

pituudet <- rnorm(1000,mean=170,sd=20) # Muodostetaan normaalijakautunut populaatio-data
sukupuoli <- factor(rep(c("N","M"),500))
populaatio <- data.frame(pituudet,sukupuoli)
f <- sample(nrow(populaatio),0.2*nrow(populaatio))
otos <- populaatio[f,]
tapply(otos$pituudet,otos$sukupuoli,mean) # Saadaan pituuden keskiarvo ja varianssi sukupuolen mukaan. Keskiarvot: M=170.54 & N= 170.92. Varianssit: M=349.47 & N=432.53
tapply(otos$pituudet,otos$sukupuoli,var)
boxplot(otos$pituudet~otos$sukupuoli) # Boxplot-kuva

library(car) # Seuraavaksi tehdään Levene-testi, jossa tutkitaan varianssien yhtäsuurutta. Nollahypoteesi varianssien yhtäsuuruudesta jää voimaan
leveneTest(otos$pituudet~otos$sukupuoli)

shapiro.test(otos$pituudet[which(otos$sukupuoli=="N")]) # Tehdään Shapiro-testi, joka kertoo, onko "pituudet" jakautunut normaalisti. Testin mukaan "pituudet" on jakautunut normaalisti
# Tehdään histogrammi, josta näemmä naisten ja miesten pituusjakaumat
hist(otos$pituudet[which(otos$sukupuoli=="N")],col="blue",main="Naisten pituudet",xlab="")
hist(otos$pituudet[which(otos$sukupuoli=="M")],col="red",main="Miesten pituudet",xlab="")

Nyt voimme tehdä yksisuuntaisen varianssianalyysi-testin

R-koodi:

oneway.test(otos$pituudet~otos$sukupuoli,var.equal=T) # P-arvo on 0.85, joten nollahypoteesi keskiarvojen yhtäsuuruudesta jää voimaan. Tämä ei ole yllätys, koska kumpikin sukupuoli sai arvonsa samasta normaalijakaumasta, jossa keskiarvo oli 170 ja keskihajonta 20.

Katsotaan vielä, etteivät mallin residuaalit poikkea hirveästi normaalijakautuneista.

m1 <- lm(otos$pituudet~otos$sukupuoli)
par(mfrow=c(2,1))
res <- resid(m1)
x <- seq(min(res),max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
hist(resid(m1),freq=F)
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)

Kuvista nähdään, että residuaalit ovat normaalijakautuneet. Residuaalien ei tarvitse kuitenkaan olla täysin normaalijakautuneita, koska otoskoon ollessa riittävän suuri Central limit theorem pätee, jolloin luottamusvälit saadaan laskettua luotettavasti. Katso tästä lisää loppukappaleesta.

Yllä olevassa kvantiili-kvantiili -käyrän kuvassa (QQ-plot) residuaalit laitetaan suuruusjärjestykseen ja niitä verrataan normaalijakauman kvantiileihin. Jos havainnot osuvat lähelle suoraa, niin havaintoja voidaan pitää normaalijakautuneina. Katsotaan vielä, miltä kuva näyttää, kun havainnot ovat positiivisesti ja negatiivisesti vinoista jakaumista.

Muodostetaan seuraavaksi uusi data, jossa sukupuolten pituudet tulevat eri normaalijakaumista ja eroavat.

pituudetN <- rnorm(500,mean=165,sd=5)
pituudetM <- rnorm(500,mean=185,sd=5)
sukupuoliN <- factor(rep("N",500))
sukupuoliM <- factor(rep("M",500))
populaatio2 <- data.frame(c(pituudetN,pituudetM),factor(c(as.character(sukupuoliN),as.character(sukupuoliM))))
colnames(populaatio2) <- c("pituudet","sukupuoli")
f <- sample(nrow(populaatio2),0.2*nrow(populaatio2))

otos <- populaatio[f,]
boxplot(otos$pituudet~otos$sukupuoli)

leveneTest(otos$pituudet~otos$sukupuoli) # Nollahypoteesi jää voimaan
hist(otos$pituudet[which(otos$sukupuoli=="N")],col="blue",main="Naisten ja miesten pituudet",xlab="",xlim=c(150,210))
hist(otos$pituudet[which(otos$sukupuoli=="M")],col="red",xlab="",add=TRUE)
legend("topright",c("Naiset","Miehet"),col=c("blue","red"),lty=c(1,1),lwd=c(5,5))

oneway.test(otos$pituudet~otos$sukupuoli,var.equal=T) # P-arvo alle 0.01, joten nollahypoteesi hylätään, eli eri sukupuolten keskipituuksissa on eroja.

Saimme siis tietoa, että eri sukupuolten keskipituuksissa on eroa. Emme kuitenkaan aina tiedä, millä ryhmillä nämä erot ovat, koska varianssianalyysin vaihtoehtoinen hypoteesi on, että ainakin kaksi keskiarvoa eroaa toisistaan tilastollisesti merkitsevästi. Esimerkissämme oli vain kaksi ryhmää, joten on helppoa nähdä, että miesten pituudet ovat suurempia. Mikäli ryhmiä on enemmän ja haluamme tietää, mitkä ryhmät eroavat toisistaan ja kuinka paljon, voimme tehdä Tukey:n HSD -testin. Huom. Testiä voi käyttää vain, jos varianssianalyysin nollahypoteesi keskiarvojen yhtäsuuruudesta on ensin hylätty. Lisäksi ryhmien varianssien tulee olla samansuuruisia. Jos näin ei ole, niin voi käyttää Games-Howell -testiä, joka löytyy R:n "userfriendlyscience" -pakkauksesta. Tehdään nyt Tukey:n HSD -testi esimerkkidatalle seuraavasti:

TukeyHSD(aov(otos$pituudet~otos$sukupuoli))

Testi antaa 95%:n luottamusvälin, jolle ryhmien keskiarvojen erotus osuu. Esimerkissämme tämä luottamusväli on naisten ja miesten pituuden erotukselle (N-M) :  -21.35, -18.50.

Tarkastetaan vielä, että mallin residuaalit ovat approksimatiivisesti normaalisti jakautuneet. Alla olevasta kuvasta nähdään, että näin myös on.

m2 <- lm(otos$pituudet~otos$sukupuoli)
res <- resid(m2)
x <- seq(-15,max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
par(mfrow=c(2,1),mar=c(3,3,2,2))
hist(res,freq=F,ylim=c(0,0.074))
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)

Muodostan seuraavaksi dataa kolmesta ryhmästä, jotta näemme, kuinka yksisuuntainen varianssianalyysi toimii, kun ryhmiä on yli kaksi.

R-koodi:

pituusR1 <- rnorm(500,mean=165,sd=5)
pituusR2 <- rnorm(500,mean=185,sd=5)
pituusR3 <- rnorm(500,mean=175,sd=5)

ryhmä1 <- factor(rep("R1",500))
ryhmä2 <- factor(rep("R2",500))
ryhmä3 <- factor(rep("R3",500))
populaatio3 <- data.frame(c(pituusR1,pituusR2,pituusR3),factor(c(as.character(ryhmä1),as.character(ryhmä2),
as.character(ryhmä3))))
colnames(populaatio3) <- c("pituudet","ryhmä")
f <- sample(nrow(populaatio3),0.2*nrow(populaatio3))

otos <- populaatio3[f,]
boxplot(otos$pituudet~otos$ryhmä)

library(ggplot2) # Tehdään tiheyskuva

ggplot(otos, aes(pituudet, fill = ryhmä)) + geom_density(alpha = 0.7)

oneway.test(otos$pituudet~otos$ryhmä,var.equal=T) # H0 hylätään
TukeyHSD(aov(otos$pituudet~otos$ryhmä),conf.level=0.95) # Saadaan alla oleva taulukko

Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = otos$pituudet ~ otos$ryhmä)

$`otos$ryhmä`
           diff        lwr       upr p adj
R2-R1  20.53805  18.720077 22.356033     0
R3-R1  10.38782   8.537313 12.238317     0
R3-R2 -10.15024 -11.991966 -8.308513     0

Tuloksista näemme, että kaikki erot ovat tilastollisesti merkitseviä 95%:n tasolla. Lisäksi luottamusvälejä tarkastelemalla näemme että:

R2:n keskipituus on suurempi kuin R1:n
R3:n keskipituus on suurempi kuin R1:n
R2:n keskipituus on suurempi kuin R3:n

Ryhmät ovat keskipituuden mukaisessa järjestyksessä siis R2, R3 ja R1.

Tarkastetaan lopuksi mallin residuaalit ja näemmä, että ne ovat approksimatiivisesti normaalisti jakautuneet.

m3 <- lm(otos$pituudet~otos$ryhmä)
res <- resid(m3)
par(mfrow=c(2,1),mar=c(3,3,2,2))
x <- seq(-15,max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
hist(res,freq=F)
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)

Tarvitseeko residuaalien olla normaalisti jakautuneita? 

Seuraavaksi katsotaan, kuinka Central limit theorem:n "lait" pätevät simuloidulle datalle. Simuloin regressiomalleja, joissa residuaalit noudattavat erilaisia jakaumia ja näytän, että poikkeustapauksia lukuunottamatta regressiomallin muuttujan regressiokerroin saa arvonsa normaalijakaumasta, riippumatta virhetermien jakaumasta. On tärkeää, että regressiokertoimet tulevat normaalijakaumasta, koska luottamusvälit ja p-arvot perustuvat vahvasti tähän oletukseen. Kun regressiokertoimet tulevat normaalijakaumasta, voimme esim. sanoa, että n.68% kerroista oikea regressiokertoimen arvo sijaitsee kerroinestimaattimme ja +-1:n keskivirheen sisällä. Jos jakauma on joku toinen, niin luottamusvälimme ja p-arvomme ovat todennäköisesti vääriä.

Regressiomalli on tässä simuloinnissa muotoa y=2+5*x+"virhetermi", eli oikea regressiokerroin on 5, joten katsotaan, kuinka lähelle tuota arvoa päästään, kun simuloidaan 1000 mallia erilaisilla virhetermin jakaumilla.

Ensimmäinen virhetermin jakauma on positiivisesti vino jakauma.

kertoimet = matrix(0,1000,2) # Luodaan matriisi, johon regressiomallin kertoimet tulevat
for(i in 1:1000) # Toistetaan alla oleva 1000 kertaa
{
    x = rnorm(100) 
    y = 2 + 5*x + rbeta(100,1,5)
    m = lm(y~x)
    kertoimet[i,] = coef(m)
}

par(mfrow=c(2,2))
hist(rbeta(100,1,5),main="Positiivisesti vino jakauma (n=100)")
qqnorm(resid(m),main="Regressiomallin residuaalit QQ Plot")
qqline(resid(m))
qqnorm(kertoimet[,2],main="Kertoimet QQ Plot")
qqline(kertoimet[,2])
hist(kertoimet[,2], main="Kertoimet Histogram")
range(kertoimet[,2]) # Kertoimet vaihtelevat välillä 4.952 - 5.044

Seuraavaksi vastaavat kuvat, kun virhetermi seuraa negatiivisesti vinoa jakaumaa. Kertoimen arvo vaihtelee välillä 4.959-5.043

Yllä näimme, että virhetermien vinous ei vääristä regressiokertoimia. Seuraavaksi katsotaan ongelmallista T-jakaumaa pienillä vapausasteilla. Tällainen virhetermin jakauma hankaloittaa regressiokertoimen estimointia. Alla kuvassa diagnostiikkakuvat T-jakautuneille virhetermeille (df=2.01). Nyt regressiokerroin saa arvoja väliltä: 1.106 - 9.787, eli hännät ovat paksuja.

Seuraavaksi lisätään hieman vapausasteiden määrää (df=3), jolloin kertoimien jakauma alkaa jo muistuttamaan enemmän normaalijakaumaa ja kerroin vaihtelee lähempänä oikeaa: 4.359 - 6.158

Seuraavaksi kuva, kun df=5. Nyt kertoimet vaihtelevat välillä 4.525 - 5.437.

Katsotaan vielä kuvat, kun virhetermit tulevat normaalijakaumasta (mean=0,sd=1). Kertoimien vaihteluväli on 4.656-5.277.

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))

Machine learning: Bagging ja puumallit

Esittelen koneoppimiseen kuuluvan bagging-menetelmän, jolla puumallien ennustusten varianssia pystytään vähentämään. Puumallien yleinen ongelma on ylisovittaminen (Overfitting), jonka takia mallit ennustavat hyvin harjoitusdataan, mutta huonosti uuteen dataan. Ne siis löytävät harjoitusdatasta rakenteita, jotka eivät kuitenkaan ole yleistettävissä muuhun dataan. Alla kuvassa malli on ylisovittanut ja hakenut punaisia pisteitä sinisten joukosta. Musta viiva antaa kuitenkin parhaan ennustustuloksen, kun mallia kokeillaan uuteen dataan. Punaiset pisteet sinisten joukossa eivät siis tässä muodossa toistu usein, kun otetaan uusia otoksia perusjoukosta.

Lähde: https://kaggle2.blob.core.windows.net/competitions/kaggle/2489/logos/front_page.png

Puumalli on mahdollista sovittaa täydellisesti harjoitusdataan, jolloin malli ennustaisi virheettömästi harjoitusdatan. Tämä johtuu siitä, että puumalli voi aina tehdä uuden oksan, kunnes kussakin oksan lehdessä on enää 1 havainto jäljellä. Esim. luottoriskiluokituksiin liittyen tämä tarkoittaisi sitä, että jokainen yksilö luokiteltaisiin loppuun asti, eli kukin yksilö kulkisi puussa omaan lehteensä, jossa ei olisi ketään muuta. Esimerkiksi kaksi lainanhakijaa, joista toisella olisi todettu maksukyvyttömyys, voitaisiin puumallilla jakaa eri lehtiin, vaikka heidän ainoa eronsa olisi yksi euro säästötilillä. Tällöin puumalli veisi hakijat samaa reittiä viimeiseen oksaan, jossa se tekisi jaon säästötilin rahasumman mukaan. Jos tällaista puuta käyttäisi luotonantoon, niin näistä kahdesta käytännössä identtisestä hakijasta vain toinen saisi lainan. Ylisovittamisen estämiseksi puumalleihin voi asettaa erilaisia ehtoja, kuten havaintojen minimimäärät lehdissä ja oksien maksimimäärät mallissa. Lisäksi ns. boostrap aggregating (bagging) on tehokas keino välttää ylisovittaminen.

Bagging

Sen sijaan, että käyttäisimme yhden puumallin harjoitusdatasta muodostamaa ennustetta, voimme ottaa harjoitusdatasta monta otosta ja muodostaa niiden pohjalta useita puumalleja. Kun näillä puumalleilla sitten ennustetaan testidataa, niin käytämmekin yksittäisten mallien sijaan niiden ennusteiden keskiarvoa. Tällä tavoin saamme luotua vakaamman ennusteen, jolloin yksittäisten mallien ylisovittuminen ei ole enää yhtä suuri ongelma (Baggingin avulla pienet muutokset datassa eivät aiheuta suuria muutoksia ennusteissa). Havainnollistan baggingin vaikutusta simuloimalla datan, jossa on kaksi muuttujaa y ja x. Lisään y:n vähän "kohinaa", jotta näemme, kuinka lähelle mallin ennusteet menevät, kun kohina häiritsee puumalleja. Muodostan sitten datasta harjoitusdatan (60% havainnoista) sekä testidatan (40% havainnoista). Alla on kuva testsetistä, johon ennustukset tehdään. Kuvassa oleva sininen viiva kuvaa parasta mahdollista ennustetta, koska se on simuloidun datan funktio ilman hajontaa aiheuttavaa satunnaista kohinaa.
R-koodi:

x <- sort(runif(300, 0, 10))
y <- sin(2.5*x)/x + rnorm(300,0,sd=0.10) 
data <- data.frame(y,x)
t <- sort(sample(nrow(data),0.6*nrow(data)))
trainingset <- dataa[t,]
testset <- dataa[-t,]
plot(testset$x, testset$y)
realvalue <- sin(2.5*x)/x
lines(x, realvalue,col="blue",lwd=2)
legend("topright",c("sin(2.5*x)/x",
"sin(2.5*x)/x + rnorm(300,0,sd=0.10)"),col=c("blue","black"),lty=c(1,NA),pch=c(NA,1),lwd=c(2,1))

Seuraavaksi muodostetaan puumalli käyttäen koko trainingset-dataa ja lasketaan sen antamien ennusteiden perusteella RMSE (Root Mean Squared Error, lisää infoa: RMSE). RMSE:lla verrataan ennusteiden osuvuutta oikeisiin arvoihin. Mitä pienempi arvo, sitä paremmin malli on ennustuksissaan onnistunut.
R-koodi:

library(rpart)
m0 <- rpart(y~.,data=trainingset,control=rpart.control(minsplit=5))
testset$pre0 <- predict(m0,testset) 
RMSE <- sqrt(mean((testset$y - testset$pre0)^2))
print(RMSE)
plot(testset$x, testset$y) # Uusi kuva, jossa m0:n ennusteet punaisella
realvalue <- sin(2.5*x)/x
lines(x, realvalue,col="blue",lwd=2)
lines(testset$x,testset$m1,col="red")
legend("topright",c("sin(2.5*x)/x",
"sin(2.5*x)/x + rnorm(300,0,sd=0.10)","m0 (RMSE=0.151)"),col=c("blue","black","red"),lty=c(1,NA,1),pch=c(NA,1,NA),lwd=c(2,1,1))

Ensimmäisen puumallin RMSE arvo on 0.151. Seuraavaksi otan traininigset-datasta viisi 20%:n satunnaisotosta (palautuksella), joista jokaisella muodostan oman puumallin. Ennustan näillä malleilla sitten testsettiin ja lasken RMSE:t. Sitten lasken mallien ennusteista bagging-menetelmän mukaisesti keskiarvot ja lasken uuden RMSE:n.
R-koodi:

a <- sample(nrow(trainingset),0.2*nrow(trainingset))
b <- sample(nrow(trainingset),0.2*nrow(trainingset))
c <- sample(nrow(trainingset),0.2*nrow(trainingset))
d <- sample(nrow(trainingset),0.2*nrow(trainingset))
e <- sample(nrow(trainingset),0.2*nrow(trainingset))
data1 <- trainingset[a,]
data2 <- trainingset[b,]
data3 <- trainingset[c,]
data4 <- trainingset[d,]
data5 <- trainingset[e,]
m1 <- rpart(y~., data=data1, control=rpart.control(minsplit=5))
m2 <- rpart(y~., data=data2, control=rpart.control(minsplit=5))
m3 <- rpart(y~., data=data3, control=rpart.control(minsplit=5))
m4 <- rpart(y~., data=data4, control=rpart.control(minsplit=5))
m5 <- rpart(y~., data=data5, control=rpart.control(minsplit=5))
pre1 <- predict(m1,testset)
pre2 <- predict(m2,testset)
pre3 <- predict(m3,testset)
pre4 <- predict(m4,testset)
pre5 <- predict(m5,testset)
RMSE1 <- sqrt(mean((testset$y - pre1)^2)) #lasketaan kullekin mallille RMSE-arvot
print(RMSE1)

Yllä olevassa kuvassa näkyvät m1-m5:n ennusteet sekä niiden yhteiset mAvg:n ennusteet. Pienimmän RMSE:n antaa malli mAvg, eli bagging on parantanut ennusteita verrattuna yksittäisiin malleihin m1-m5. Alla on vielä yllä oleva kuva ilman sinistä viivaa, jotta bagging:n merkitys näkyy selkeämmin. Kuvasta näkee, kuinka vaihtelevia yksittäisten puumallien ennustukset ovat, vaikka ne ovat kaikki peräisin saman trainingsetin arvoista. Baggingin avulla näistäkin ailahtelevista puumalleista saadaan luotua vakaa ennustus (ennusteet ovat vakaita, koska baggingin takia tapahtuu konvergoitumista, eli lisäämällä otosten määrää trainingsetistä mAvg:n tulokset eivät radikaalisti poikkea nykyisestä), koska keskiarvotus jättää yksittäisten puumallien rajut ylisovittamiset vähemmälle huomiolle. Lisäksi, jos enemmistö puumalleista on ennusteessaan melko oikeassa, niin nämä ennusteet saavat enemmän painoarvoa (esim. katso kuvaa, kun x=3). Tietysti, jos enemmistö puumalleista olisi väärässä, niin myös huonot ennusteet painottuisivat. Kuitenkin Ensemble learning -menetelmät perustuvat ajatukseen, että keskinkertaisista malleista saadaan luotua hyviä yhdistämällä, ja pahimmassakin tapauksessa yhdistetty malli on yksittäisen mallin veroinen, eli ei muutu sitä huonommaksi. Kahden mallin yhdistäminen voi olla hyvä ajatus erityisesti, jos kummankin mallin ennusteet poikkeavat toisistaan ja mallit ovat melko varmoja oikeissa ennusteissaan, kun taas epävarmoja ollessaan väärässä. Tällaiset mallit yhdistämällä niiden oikeaksi menneet ennusteet pysyvät sellaisina myös keskiarvotusten jälkeen ja lisäksi osa toisen mallin epäonnistumisista korjaantuu, koska toinen malli ennusti ne oikein. Näin väärien ennustusten määrä onkin yhdistetyssä mallissa pienempi kuin yksittäisissä. Lisää tieto löytyy mm. täältä

"A necessary and sufficient condition for an ensemble of classifiers to be more accurate than any of its individual members is if the classifiers are accurate and diverse. An accurate classifier is one that has an error rate of better than random guessing on new x values Two classifiers are diverse if they make different errors on new data points"

Lähde: http://web.engr.oregonstate.edu/~tgd/publications/mcs-ensembles.pdf

Yhdistetty malli mAvg on tällä hetkellä m0:n veroinen ennustuksissaan, joten baggingista ei ole ollut käytännössä vielä huomattavaa hyötyä. Seuraavaksi otan viisi uutta satunnaisotosta trainingset-datasta ja muodostan uudet puumallit, joiden ennusteet yhdistän m1-m5:n kanssa. R-koodi on vastaava kuin aiempi. Nyt RMSE on 0.141. Alla kuva, jossa vanha ja uusi malli.

Alla vielä kuva koneoppimiseen liittyvän Random forest -menetelmän ennusteista. RMSE on tälle mallille 0.119. Random forest -menetelmässä algoritmi muodostaa useita puumalleja (tässä 5000 kpl) ja keskiarvottaa niiden ennusteet (luokittelevassa ennustamisessa se käyttää ennusteiden moodia). Puumallit tehdään bagging-menetelmällä, mutta myös muuttujista käytetään vain osaa kussakin puumallissa (muuttujat kuhunkin jakoon (split) valitaan sattumanvaraisesti). Oletuksena R:ssa muuttujia "n" käytetään kussakin jaossa neliöjuuri(n) (luokittelu) tai n/3 (regressio) määrä. Random forest on kätevä menetelmä, jolla saa parempia ennusteita kuin yksittäisillä puumalleilla. Random forest ei kuitenkaan juuri paranna ennusteita verrattuna logistiseen regressioon, kun ennustetaan binääristä selitettävää muuttujaa. Random forest löytää kuitenkin muuttujien välisiä yhteisvaikutuksia (Interactions) ja on siten avuksi logistisen regression rinnalla. Näiden mallien yhdistämisestä kirjoittelen jatkossa.

Lineaarinen regressio ja Interaction-termit

 Lineaarinen regressio sopii tilanteisiin, joissa haluamme tietää, kuinka eri muuttujat vaikuttavat selitettävään muuttujaan määrällisesti. Lineaarisella regressiolla voi esimerkiksi tutkia, miten ihmisten tulotaso vaikuttaa heidän ostamansa auton hintaan. Lineaarisella regressiolla voi laskea myös todennäköisyyksiä jollekin binääriselle tapahtumalle, mutta logistinen regressio sopii tällaiseen huomattavasti paremmin mm. koska logistinen regressio pakottaa todennäköisyyden 0:n ja 1:n välille, kun taas lineaarinen regressio ennustaa myös ulkopuolen (< 0 & >1) arvoja. Mikäli ratkaisee tämän ongelman asettamalla alle nollan arvot nollaksi ja yli yhden arvot yhdeksi, niin syntyy äärimmäisiä 0:n ja 1:n todennäköisyyksiä, eli malli ennustaisi tapahtumia täysin varmoiksi. Lineaarinen regressiomalli rikkoo myös esim. virhetermien normaalisuusehdon, joten logistinen malli on parempi vaihtoehto. Logistisen regressiomallin ennusteet  käyttäytyvät pehmeämmin kuin lineaarisen regressiomallin ("ei lineaarisesti"), kun lähestytään todennäköisyyksien ääripäitä. Alla olevasta kuvasta näkyy, kuinka logistisen mallin todennäköisyydet pehmenevät ääripäissä.

Lähde: http://www.appstate.edu/~whiteheadjc/service/logit/intro.htm

Vaikka lineaarinen regressio ei toimi logistisen mallin veroisesti todennäköisyyksien mallintamiseen, on se hyödyllinen muihin tarkoituksiin. Sen avulla voi saada intuitiivista tietoa nopeasti eri muuttujien välisistä mahdollisista yhteyksistä. Lineaarista mallia käytettäessä on tärkeää, että selittävän ja selitettävän muuttujan välillä vallitsee lineaarinen yhteys. Mikäli yhteys ei ole lineaarinen, on mahdollista käyttää esim. neliöjuuri-muunnettua muuttujaa normaalin muuttujan kanssa, jolloin epälineaarisia yhteyksiä voidaan mallintaa. Tällaisessa tapauksessa täytyy kuitenkin muistaa tarkastella, onko tämä teoreettisesti oikein. Esim. tulotason kasvu voi vaikuttaa ostetun auton hintaan nostavasti tiettyyn pisteeseen "x" saakka, kunnes kasvu hiipuu tai jopa laskee. Tällöin suoran viivan vetäminen ei kuvaa yhteyttä oikein, koska mallin mukaan auton hinta kasvaa vielä, kun menemme kriittisen tulotason "x" yli. Tällaisessa tapauksessa voi harkita selittävän muuttujan luokittelua, jos haluaa pitää mallin helpommin tulkittavana, eli ei lisäile epälineaarisia yhteyksiä löytäviä muuttuja-muunnoksia. Yksi regressiomallien hienouksia on, että voimme tutkia, kuinka muuttuja "x vaikuttaa selittettävään muuttujaan "y" kun kontroillaan muuttujaa "z". Eli voimme esim. tutkia, kuinka tulotaso vaikuttaa auton hintaan, kun pidämme sukupuolen vakiona. Kontrollointi tapahtuu lineaarisessa mallissa seuraavasti:

Haluamme tietää mallin (y=x+z) kertoimen x:lle. Ensin ratkaistaan x:n vaikutus y:hyn regressoimalla x z:lla ja tallentamalla residuaalit. Tämän jälkeen y regressoidaan näillä residuaaleilla. Poistamme siis ensin x:sta z:n vaikutuksen ja sitten vasta regressoimme y:n tämän x:n, mutta z:sta riippumattoman osan kanssa. Vastaavalla tavalla saamme myös z:n kertoimen. Yksi regressiomallin yleisistä ongelmista (Omitted variable bias) aiheutuu juurikin tämän kontrolloinnin puutteesta, eli jos tärkeä muuttuja jätetään mallista pois, niin muuttujien regressiokertoimet voivat olla vääristyneitä. Esim. x:n kerroin voi olla suurempi kuin mitä se olisi jos olisimme ottaneet z:n vaikutuksen huomioon.

Alla vielä vastaava esimerkki R-koodilla:

sukupuoli <- (rbinom(100,1,0.5)) #luodaan sukupuoli-muuttuja Bernoullin jakaumasta
tulot <- rnorm(100,15000,100)+150*sukupuoli # luodaan toinen muuttuja, joka korreloi sukupuolen kanssa
y=10000+(0.1*tulot)+(500*sukupuoli) # y kuvaa auton hintaa. Tästä näemme oikeat regressiokertoimet.

#Ratkaistaan ensin tulojen regressiokerroin

regTs <- lm(tulot~sukupuoli)
residuaalitT <- residuals(regTs)
m1 <- lm(y~residuaalitT)
round(coef(m1),1) # Saadaan kerroin 0.1. Round-komento pyöristää tässä luvut yhteen desimaaliin.
regSt <- lm(sukupuoli~tulot) # Ratkaistaan vastaavasti sukupuolen regressiokerroin 
residuaalitS <- residuals(regSt)
m2 <- lm(y~residuaalitS)
round(coef(m2),2)# Saadaan 500

sum(residuaalitT*y)/(sum(residuaalitT^2)) #Regressiokertoimet saadaan laskettua myös näin
sum(residuaalitS*y)/(sum(residuaalitS^2)) 

sum(y-0.1*tulot-500*sukupuoli)/100 # Intercept saadaan lopuksi tästä kaavasta
#Lopputulos on siis y=10000+(0.1*tulot)+(500*sukupuoli)

Omitted variable bias tulee ilmi, kun jätämme sukupuolen pois tarkastelusta:

m3 <- lm(y~tulot)
round(coef(m3),1)
#saadaan tulojen kertoimeksi 0.3

Koska tulot ja sukupuoli korreloivat positiivisesti (cor=0.63) ja sukupuoli vaikuttaa auton hintaan positiivisesti, niin tulojen kerroin on liian suuri, jos jätämme sukupuolen pois mallista. Vastaavasti, jos sukupuoli vaikuttaisi auton hintaan negatiivisesti (y=10000+(0.1*tulot)-(500*sukupuoli)), niin tulojen kerroin olisi alaspäin vääristynyt -0.1.

Interaction-termit

Yllä olevassa esimerkissä tulojen kasvu vaikuttaa auton hintaan kullakin sukupuolella samoin . On kuitenkin mahdollista, että eri sukupuolten väliset preferenssit eroavat ja toinen sukupuoli käyttää tulojen noususta suhteellisesti enemmän  rahaa autoonsa. Tällainen yhteys voidaan mallintaa lisäämällä malliin ns. Interaction-muuttuja, jolloin tulojen muutos voi vaikuttaa eri tavoin auton hintaan eri sukupuolilla. Interaction-muuttuja lisätään R:lla seuraavasti: 

m4 <- lm(y~tulot+sukupuoli+tulot:sukupuoli)

Nyt tulojen kerroin riippuu sukupuolesta. Esim. sukupuolella "1" tulojen kerroin voi olla 0.5 ja sukupuolella "0" 0.7. 

Alla on kuva datasta, jossa tulot vaikuttavat ostetun auton hintaan (y) eri tavoin. Selvästi sukupuoli "0" käyttää autoonsa enemmän rahaa tulojen noustessa. 

Jos emme ota sukupuolten välistä eroa huomioon, saamme seuraavanlaiset regressioviivat

Vastaava kuva, kun lisäämme interaction-termin malliin:

Interaction-termin lisääminen parantaa mallin ennusteita selkeästi.

R taloustieteellisen tehtävän tukena

Esittelen nyt R:n käyttöä erään taloustieteellisen tehtävän kanssa. R:lla voi tehdä laskutoimituksia ja erityisen hyödyllinen se on kuvien tekemisessä. Tehtävä on seuraavanlainen: 

Hyödykkeen  X

  kysyntä D  on        p = 7 - q/2

  ja tarjonta S on      p = 3/4 + q/2

  määrää hyödykkeen  tasapainohinta ja kysytty  määrä.

Määrää myös suurin mahdollinen veron tuotto, jonka julkinen valta voi 

saada asettamalla hyödykkeelle  t:n suuruisen yksikköveron.

Huom. Käytän * merkkiä tarkoittamaan kertomista, eli 2*2 on kaksi kertaa kaksi. Vastaavasti ^ on potenssimerkki, eli 2^2 tarkoittaa kaksi potenssiin kaksi.

Tässä tehtävässä tulee asettaa p (kysyntä)=p (tarjonta), eli kun piirtää kysyntäkäyrän ja tarjontakäyrän, niin kohdassa p(kysyntä)=p(tarjonta) käyrät leikkaavat (kyseisellä hinnalla kysyntä on sama kuin tarjonta, joten markkinat ovat tasapainossa.

Asetetaan p(kysyntä)=p(tarjonta) ja ratkaistaan määrä q → 7-q/2=3/4+q/2 → q=25/4

Markkinoiden tasapainomäärä on siis 25/4 kappaletta hyödykettä. Markkinahinta saadaan nyt sijoittamalla saatu arvo q=25/4 kysyntäkäyrän- tai tarjontakäyrän yhtälöön (kumpikin tuottaa saman hinnan, koska ratkaisimme määrän q asettamalla ensin hinnat samaksi):

Sijoitetaan tasapainomäärä kysyntäkäyrän yhtälöön p=7-(25/4)*1/2 → p=(56/8)-(25/8) → p=31/8

Vastaus on siis: Hyödykkeen X tasapainohinta: p=31/8 ja Kysytty määrä: q=25/4

R:lla voidaan ratkaista tasapainohinta ja -määrä sekä luoda kuva markkinatasapainosta seuraavasti:

R-Koodi:

kysyntä <- function(q) { 7-(1/2)*q } #Kysyntäkäyrän funktio

plot(kysyntä,from=1,to=10,xlab="Määrä q",ylab="Hinta p",lwd=2) # Kuva kysyntäkäyrästä

tarjonta <- function(q) {3/4+(1/2)*q } # Tarjontakäyrän funktio

curve(tarjonta,add=TRUE,col="blue",xlab="Määrä",lwd=2) # Lisätään kuvaan tarjontakäyrän

legend("topright",c("Kysyntäkäyrä","Tarjontakäyrä"),lty=c(1,1),col=c("black","blue"),lwd=c(2,2))

g <- function(q) { 7-(1/2)*q-3/4-(1/2)*q } # Kysyntä=Tarjonta

uniroot(g,lower=1,upper=10) # Ratkaistaan edellinen funktio.

kysyntä(6.25) # Ratkaistaan tasapainohinta sijoittamalla tasapainomäärä kysyntäkäyrän funktioon

text(7.5,3.9,"(6.25,3.875)") # Lisätään kuvaan tasapainohinnan ja -määrän koordinaatit

Seuraavaksi ratkaistaan suurin mahdollinen veron tuotto, jonka julkinen valta voi saada asettamalla hyödykkeelle t:n suuruisen yksikköveron. Yksikkövero tarkoittaa, että jokaisesta myydystä tuotteesta maksetaan vero. Julkinen valta tietää kysyntäkäyrän (p=7-q/2) sekä tarjontakäyrän (p=3/4+q/2) ja on päättänyt asettaa veron yrityksille. Yritykset laittavat yksikköveron T tarjontakäyrän yhtälöön, eli yhtälö on nyt muotoa: p=3/4+q/2+T

Veron seurauksena tarjontakäyrä nousee siis veron verran ylös, jolloin kysyntä vähenee. Tästä aiheutuu se, ettei markkinahinta nouse veron verran. 

Ratkaistaan ensin tasapainomäärä, kun p(kysyntä)=p(tarjonta+vero). Tämän jälkeen muodostamme verotulojen funktion, jonka derivoimme veron T suhteen. Derivaatta asetetaan nollaksi, jolloin voimme ratkaista verotulot maksimoivan veron määrän. Verotulojen funktio on määrä*vero, eli jos vaihdettu määrä on 8 ja vero 1 niin valtio saa verotuloja 1*8=8. Haluamme siis maksimoida nuo verotulot. Valtio voi nostaa veroa esim. tuosta yhdestä eurosta kahteen euroon, mutta tällöin kysytty määrä vähenee, eli se ei olekaan enää tuo kahdeksan. Se, kuinka paljon kysyntä vähenee, riippuu kysyntäkäyrän muodosta. Kysynnän hintajousto liittyy tähän, eli jos kysyntä laskee hinnan noustessa suhteessa enemmän kuin hinta (joustava kysyntä), niin käytetty rahamäärä vähenee, ja täten myös verotulot tippuvat. Toistaalta, jos kysyntä vähenee hinnan nousteessa suhteessa vähemmän kuin hinta, niin käytetty rahamäärä kasvaa ja verotulot kasvavat (joustamaton kysyntä).

Eli nyt ratkaistaan tehtävä asettamalla p(kysyntä)=p(tarjonta+vero) ja ratkaistaan q →

7-q/2=3/4+q/2+T → q=25/4-T

Verotulojen funktio on qT eli kerrotaan yhtälön kummatkin puolet T:lla ja saadaan → qT=(25/4)T-T^2

Nyt meillä on verotulojen funktio, mutta emme tiedä, mikä veron määrä T maksimoi tämän funktion. Meidän täytyy derivoida verotulojen funktio veron T suhteen. Derivointi siis kertoo meille funktion käyttäytymisestä, eli jos derivaatta on positiivinen, niin itse funktio kasvaa. Derivaatan nollakohdassa derivaatan merkki vaihtuu, joten nämä ovat funktion maksimeja tai minimejä. Jos derivoimme verotulojen funktion T:n suhteen ja emme aseta derivaattaa nollaksi, niin emme tiedä, missä kohtaa funktio saa “ääriarvonsa” eli verotulot maksimoituvat.

Derivoidaan siis qt ja asetetaan derivaatta nollaksi. d(qT)=(25/4)-2T=0

Nyt voimme ratkaista yhtälöstä T:n → 2T=25/4 → T=25/8

Saamamme T on siis maksimiarvo, koska verotulojen funktion toinen derivaatta T:n suhteen on negatiivinen (-2)

Voimme katsoa verotulojen funktiosta verotulojen määrän tällä T:n arvolla. Se on 9,77

(qt=(25/4)*(25/8)-(25/8)^2=9,77)

Voimme myös ratkaista uuden markkinatasapainon hinnan ja määrän. Ratkaistaan ensin määrä q sijoittamalla T=25/8 yhtälöön q=25/4-T

eli yhtälö jonka saimme, kun asetimme p=p. Ratkaistaan: q=25/4-25/8 → q=25/8

Ja nyt voimme tarkistaa tasapainohinnan kysyntä- tai tarjontakäyrän yhtälöstä. Esim. kysyntäkäyrästä: p= 7-q/2 → p=7-(25/8)(1/2) → p=7-25/16 =5,44

Näemme siis, että veron takia uusi tasapainomäärä (25/4 > 25/8) on pienempi ja hinta suurempi 31/8<5,44) kuin ennen veroa.

On myös hyvä havaita, että tuottajien saama hinta on markkinahintaa pienempi, koska he maksavat veron. Tuottajien saama hinta on siis P-T. Hinnan voimme ratkaista tuottajien tarjontakäyrästä → P=3/4+q/2+T → P-T=3/4+q/2

Tasapainossa määrä oli 25/8, joten tuottajien saama hinta on P-T=3/4+25/16=37/16

Alla vielä Kuva uudesta tasapainosta ja kuvan R-koodi.

R-koodi:

kysyntä <- function(q) { 7-(1/2)*q }

plot(kysyntä,from=1,to=10,xlab="Määrä q",ylab="Hinta p",lwd=2)

tarjonta <- function(q) {3/4+(1/2)*q }

curve(tarjonta,add=TRUE,col="blue",xlab="Määrä",lwd=2)

legend(6,6.55,c("Kysyntäkäyrä","Tarjontakäyrä","Tarjontakäyrä(T)"),lty=c(1,1,1),col=c("black",

"blue","red"),lwd=c(2,2,2))

g <- function(q) { 7-(1/2)*q-3/4-(1/2)*q }

uniroot(g,lower=1,upper=10)

kysyntä(6.25)

text(7.5,3.9,"(6.25,3.875)")

g <- function(q) { 7-(1/2)*q-3/4-(1/2)*q-25/8 } #Kysyntä=Tarjonta

uniroot(g,lower=1,upper=10) #Uusi markkinatasapainomäärä

kysyntä(3.125) # Uusi markkinatasapainohinta

tarjontaverolla <- function (q) {3/4+(1/2)*q+25/8 }

curve(tarjontaverolla,add=TRUE, col="red",lwd=2) # Lisätään uusi tarjontakäyrä(T)

text(4.45,5.45,"(3.125,5.438)")