Verohallinnon avoimen datan tarkastelua

Tässä tekstissä tutkin Verohallinnon avointa dataa, josta löytyvät tiedot julkisten yhteisöjen tuloverotuksesta vuodelta 2014. Datasta löytyvät mm. y-tunnus, verotuskunta, verotettavat tulot, maksuunpannut verot, ennakkoverot, veronpalautukset sekä jäännösverot. Hyödynnän lisäksi verovelkarekisteriä, josta löytyy tieto, onko yrityksellä verovelkaa yli 10000 euroa. Haluan erityisesti selvittää, onko jäännösverolla ja verovelalla yhteyttä, eli toimiiko jäännösvero indikaattorina tulevaisuuden verovelalle. Suuri jäännösvero voi kertoa yrityksen maksuvaikeuksista tai siitä, että yrityksen tulos on ylittänyt odotukset, jolloin ennakkoveron olisi pitänyt olla suurempi. Toisaalta tällöin yritys olisi voinut maksaa puuttuvat ennakkoveronsa ennakontäydennysmaksuna, jolloin korollista jäännösveroa ei olisi syntynyt. Lisätietoa ennakkoverotuksesta löytyy mm. Verohallinnon sivuilta täältä.

Rajataan tarkastelu aluksi yrityksiin, joilla on verotettavaa tuloa yli 10000 euroa. Näin saadaan aineistoa rajattua suurempiin yrityksiin ja niihin, joilla on jäännösveroa. Lasketaan sitten yhteisö- ja jäännösveroprosentti seuraavasti:

R-koodi:

data <- read.csv(file.choose(), header=F,sep=";",dec=",",skip=1)
colnames(data) <- c("Verovuosi","Ytunnus","Verovelvollisennimi","Verotuskunta","Verotettavatulo","Maksuunpannutverotyhteensa","Ennakotyhteensa", "Veronpalautus","Jaannosvero")

data1 <- data[data$Verotettavatulo>10000&data$Maksuunpannutverotyhteensä>0,]

data1$verotusuhdekaikki <- data1$Maksuunpannutverotyhteensa/data1$Verotettavatulo*100

data1$jaannossuhdekaikki <- data1$Jaannosvero/data1$Maksuunpannutverotyhteensä*100

Katsotaan yhteisöjen tuloveroprosenttien persentiilejä ja nähdään, että mukana on myös yrityksiä, joilla yhteisöveroprosentti on reippaasti yli 20.

Katsotaan vielä histogrammia. Rajataan tarkastelu yrityksiin, joilla tuloveroprosentit ovat välillä 15 ja 30. Yritysten määrä tippuu 57816:sta 57113:een.

Lasketaan seuraavaksi kunnille prosenttiosuudet, joista selviävät, kuinka monta prosenttia kunnan yrityksistä on saanut yli 50%:n verojäämän. Tätä varten tehtään funktio, jota käytetään yhdessä lapply-komennon kanssa. Ennen funktiota muodostetaan uusi dataframe, jossa ovat vain kuntien nimet. Tätä varten pitää Verotuskunta-muuttuja ensin jakaa kahteen osaan, koska se sisältää sekä kuntakoodin että nimen. Tällainen kahtiajako on kätevä tehdä reshape-pakkauksen transform-komennolla.

install.packages("reshape")  # Tehdään muuttujan kahtiajako

library(reshape)

data3 = transform(data3, verotuskunta = colsplit(data3$Verotuskunta, split = " ", names = c('numeric', 'nimi','koodi')))

data4 <- data.frame(unique(data3$verotuskunta.nimi)) # Muodostetaan uusi dataframe, jossa kuntien nimet vain kerran

colnames(data4) <- "Kunta"

data4$Kunta <- as.character(data4$Kunta)

osuus <- function(i) {                                   # Tehdään funktio

a <- length(data3$verotuskunta.nimi[data3$verotuskunta.nimi==i])

b <- length(data3$jaannossuhdekaikki[data3$verotuskunta.nimi==i&data3$jaannossuhdekaikki>=50])

c=b/a*100

return(c)

}

data4$Osuus <- lapply(data4$Kunta,osuus) # Lasketaan prosenttiosuudet joka kunnalle

Tehdään karttakuva, josta lasketut prosenttiosuudet selviävät. Kuvaa varten täytyy ensin hakea koordinaatit, jotka löytyvät esim. täältä. Sitten datat yhdistetään kunnan perusteella. 

R-koodi:

koordinaatit<- read.csv(file.choose(),header=TRUE,dec=".",sep=";")

colnames(koordinaatit)[2] <- "Kunta"

datayhd <- merge(koordinaatit,data4, by="Kunta")

datay <- datayhd[!duplicated(datayhd$Kunta),]

datay$Osuus <- as.numeric(datay$Osuus)

g <- sample(nrow(datay),0.05*nrow(datay))    #Otetaan kuvaa varten 5% havainnoista

datas <- datay[g,]

 map <- get_map(location = 'Finland', zoom = 6)  #Ladataan Suomen kartta 

  mapPoints <- ggmap(map)+

   geom_point(aes(x=Longitude,y=Latitude,size=Osuus),data=datas,alpha=0.5,colour="blue")

  mapPoints

mapPointsLegend <-mapPoints +scale_size(range=c(7,25)) +ggtitle("Jäännösveroprosentti >50: Osuudet Suomen kunnissa")

mapPointsLegend

Katsotaan vielä muutamia suurimpien ja pienimpien osuuksien kuntia. Helsingissä Osuus on ~25%. 

head(datay[order(datay$Osuus),])

          Kunta Postal.Code Latitude Longitude Osuus

123    Hailuoto       90480  65.0000   24.7167     0

624   Kivijärvi       43801  63.0667   25.0500     0

705    Kumlinge       22820  60.2667   20.7833     0

896     Luhanka       19950  61.7833   25.7000     0

899  Lumparland       22630  60.1167   20.2500     0

1298  Savukoski       98800  67.2925   28.1581     0

head(datay[order(-datay$Osuus),])

             Kunta Postal.Code Latitude Longitude     Osuus

909      Merijärvi       86220  64.3000   24.4500 100.00000

118           Geta       22340  60.3833   19.8500  75.00000

1205    Rautavaara       73901  63.4833   28.3000  50.00000

1047 Pelkosenniemi       98500  67.1108   27.5106  41.66667

117          Föglö       22710  60.0167   20.4167  37.50000

1536         Vaala       91700  64.4333   26.8000  36.84211

Seuraavaksi tutkitaan, onko suuren jäännösveroprosentin yrityksillä suurempi todennäköisyys saada verovelkaa. Rajaan tutkittavan joukon verotettavan tulon mukaan välille 50000-80000. Tästä joukosta valitsen ne, joilla on yli 99%:n jäännösvero ja ne, joilla on alle 1%:n jäännösvero. Näistä kummastakin joukosta otan sitten 30 kappaleen satunnaisotannan. 

data33 <- data3[data3$Verotettavatulo>=50000 & data3$Verotettavatulo<=80000,]

datajaannos99 <- data33[data33$jaannossuhdekaikki>99,]

datajaannos1 <- data33[data33$jaannossuhdekaikki <1,]

a <- sample(nrow(datajaannos99),30)

b <- sample(nrow(datajaannos1),30)

datajaannos99t <- datajaannos99[a,]

datajaannos1t <- datajaannos1[b,]

Yli 99%:n jäännösveron otoksessa verovelkaisia oli 7 kappaletta, eli 23.3%:a kaikista. 

Alle 1%:n jäännösveron otoksessa verovelkaisia ei ollut yhtään kappaletta. 

Tehdään vielä Fisherin exact testi. Testillä voidaan testata, onko populaatioiden dikotomisten muuttujien arvojen suhteellisissa osuuksissa eroja.

data = matrix(c(7,23,0,30), ncol=2)
fisher.test(data,alternative = "two.sided")

Fisher's Exact Test for Count Data
data:  data
p-value = 0.01054
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.649705      Inf
sample estimates:
odds ratio 
 Inf 

Ratkaistaan myös luottamusväli yli 99%:n ryhmän verovelallisten osuudelle.

prop.test(x=7,n=30,conf.level = 0.95)

1-sample proportions test with continuity correction

data:  7 out of 30, null probability 0.5

X-squared = 7.5, df = 1, p-value = 0.00617

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

 0.1063502 0.4270023

sample estimates:

        p 

0.2333333 

Laskin verovelallisten osuuden lopuksi myös 45-55 %:n jäännösveroryhmälle. 30:n yrityksen otoksessa oli verovelallisia 3 kappaletta, eli 10%:a kaikista.

Verovelkarekisterissä näkyvät tiedot vain viimeiseltä 6 kuukaudelta, joten käytetyn datan (2014) ja verovelkadatan välillä on yli vuoden ero. Tämä antaa viitteitä siitä, että suuri jäännösveroprosentti tänään voi merkitä suurempaa verovelkaisuuden todennäköisyyttä tulevaisuudessa. Ei ole tietenkään suuri yllätys, että alle 1%:n jäännösverolaisilla ei ollut verovelkaa, koska hehän ovat jo maksaneet veronsa ennakkoveroina, jolloin verovelkaa ei voi sillä hetkellä syntyä. Tämä tulos on kuitenkin hyödyllinen, koska se ennakoi verovelvollisen käytöstä tulevaisuudessa. 

R ja suuret aineistot

Esittelen nyt suurten aineistojen käsittelyä R:llä. Käytän harjoitusaineistona Englannin ja Walesin asuntojen hintadataa, joka löytyy täältä. Datassa on tietoja asuntojen hinnoista Englannissa ja Walesissa vuodesta 1995 lähtien. Koko data sisältää yli 22 miljoonaa riviä ja 16 muuttujaa. Muuttujina löytyvät mm. asunnon hinta, uusi/vanha asunto, tontti oma/vuokralla, myyntipäivämäärä ja asuntotyyppi. Asuntojen koosta ei ole kuitenkaan dataa. Olen kiinnostunut tutkimaan, onko kalleimpien asuinpaikkojen keskihinnoissa ollut eroja vuosina 2014 ja 2015. Rajaan tarkastelun omakotitaloihin, joilla on oma tontti.

Aloitetaan lataamalla data R:ään. Tähän kannattaa käyttää data.table-pakettia, josta löytyy fread-komento. Tämä on huomattavasti nopeampi tapa ladata suuri määrä dataa kuin esim. read.csv. Tarvitsen suuresta datamäärästä dataa vuosille 2014 ja 2015, joten jätän osan datasta lataamatta käyttämällä skip-komentoa. Dataa latautuu 2.5 miljoonaa riviä. Tällaisen määrän dataa R pyörittää vielä hyvin 6GB:n keskusmuistilla.

library(data.table)
data <- fread(file.choose(),header=F,skip=19000000)
colnames(data) <- c("ID","Price","Date","Postcode","Propertytype","Oldnew","Duration","PAON","SAON","Street","Locality","Towncity","District","County","PPD","RS")
data$Price <- as.numeric(data$Price)
data$Date <- as.Date(data$Date, format="%Y-%m-%d")
data$Year <- as.numeric(format(data$Date,"%Y")) # Teen vuosi-muuttujan, jotta vuodet on helppo valita jatkossa.

Nyt kun dataan on tehty sopivat muunnokset, lasketaan ensin uusien omakotitalojen keskihinnat vuosille 2014 ja 2015.

mean(data$Price[data$Year=="2014"  & data$Oldnew=="Y"&data$Duration=="F"&data$Propertytype=="D"]) #Keskihinnat 2014 ja 2015
mean(data$Price[data$Year=="2014"  & data$Oldnew=="Y"&data$Duration=="F"&data$Propertytype=="D"])

Keskihinnat ovat 2014: ~358000 ja 2015: ~373000

Tarkastellaan seuraavaksi keskihintoja maakuntatasolla.

data2015<- do.call(data.frame,aggregate(data$Price[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"],list(data$County[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"],data$Oldnew[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"]),function(x) c(mean=mean(x),length=length(x))))

data2014<- do.call(data.frame,aggregate(data$Price[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"],list(data$County[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"],data$Oldnew[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"]),function(x) c(mean=mean(x),length=length(x))))

Otetaan tarkastelusta pois alueet, joista on alle 30 havaintoa.

data2015l <- subset(data2015, data2015$Length>=30,)
data2014l <- subset(data2014, data2014$Length>=30,)

Järjestetään data keskiarvon mukaan suurimmasta pienimpään.

data2015l <- data2015l[order(-data2015l$Mean),]
data2014l <- data2014l[order(-data2014l$Mean),]

Katsotaan heatmappia, jossa on alhaisen ja korkean keskihinnan maakuntia vuodelta 2015.

datax <- data2015l[c(1:6),]
datay <- data2015l[c(184,189,185,200,191,195.192),]
data5 <- rbind(datax,datay)

library(ggplot2)
ggplot(data = data5, aes(x = Oldnew, y=County, fill=mean)) + 
geom_raster() + scale_fill_gradient2(low="red", mid="yellow", high="green",na.value="grey") +
labs(fill='Keskihinta',title="Asuntojen keskihinnat")

Tehdään seuraavaksi  kuva, jossa näkyvät neljä keskihinnaltaan kalleinta uusien omakotitalojen aluetta vuodelta 2015.

data2015 <- subset(data2015l,data2015l$Oldnew=="Y",)
data2015r <- data2015[,"Mean",drop=FALSE]  #drop=FALSE pitää data.frame -muodon,
row.names(data2015r) <- data2015$County  
data2015tm <- t(as.matrix(data2015r[1:4,,drop=FALSE]))


b <- barplot(data2015tm,,col=c("3"),las=1,yaxt="n",ylab="£",cex.names=0.62,,main="Englannin ja Walesin kalleimmat asuinalueet vuonna 2015 \n (Vertailu: Uusien omakotitalojen keskihinnat)")
axis(2,seq(0,1600000,by=200000))
legend("topright",legend=c("Keskihinta","Kaikkien uusien omakotitalojen keskihintataso"),col=c("3","1"),pch=c(19,NA),lty=c(NA,1))
text(b,60000,round(data2015tm,0))
abline(h=372630,lwd=2)

Katsotaan lopuksi, miten näiden alueiden keskihinnat ovat kehittyneet aiemmasta vuodesta.

data20157<- data2015[1:4,c("County","Mean")]
data2014 <- subset(data2014l,data2014l$Oldnew=="Y",)
yhd <- merge(data20157,data2014,by="County",all.y=F)
yhd$erotus20152014 <- yhd$Mean.x -yhd$Mean.y
yhd$muutosprosentti <- (yhd$erotus20152014/yhd$Mean.y)*100
colnames(yhd) <-c("County","mean2015","Oldnew","mean2014","length","erotus20142014","muutosprosentti")
yy <- yhd[,"muutosprosentti",drop=FALSE]
row.names(yy) <- yhd$County

muutosprosentti <- t(as.matrix(yy))
b <- barplot(muutosprosentti,col=c("5"),las=1,yaxt="n",ylab="%",cex.names=0.62,main="Keskihintojen muutos vuosina 2014-2015")
axis(2,seq(-20,25,by=1))
legend("topright",legend=c("Muutosprosentti"),col=5,pch=19)

text(b,1,round(muutosprosentti,2))

Vaihtelua vuosien välillä on siis melkoisen paljon. Ongelmana tässä datassa on, ettemme tiedä asuntojen kokoa. Asuntojen hintojen kehitystä olisi luontevinta tutkia vertailemalla neliöhintoja alueittain. Suuret hintavaihtelut selittyvät luultavasti osittain sillä, että tutkimme uusia omakotitaloja. Tällaisia taloja myytiin esim. Suur-Lontoossa 223 kappaletta vuonna 2015. Windsorissa vastaava luku oli vain 35. Tällaisilla kalliilla ja harvinaisilla taloilla hintavaihtelut ovat luultavasti suurempia kuin edullisemmilla taloilla, joita on markkinoilla tarjolla koko ajan. Alla on taulukko, jossa on muutosprosentteja vanhoille omakotitaloille kalleimmilla alueilla. Näemme, että luvut ovat Suur-Lontoolle, Poolelle, Surreylle ja Windsorille huomattavasti maltillisemmat. Lisäksi Count-sarakkeiden myydyistä kappalemääristä näemme, että Suur-Lontoossa ja Windsorissa myydään runsaasti enemmän vanhoja kuin uusia omakotitaloja, mikä ei tietenkään ole yllätys. Windsorista käytetyn omakotitalon sai vuonna 2015 halvimmillaan 94000 punnan hintaan, kun taas uudesta joutui maksamaan 325000 puntaa. Kallein myyty omakotitalo maksoi siellä 6.4 miljoonaa puntaa.

2015-2014 muutosprosentit vanhoille omakotitaloille

County	Mean2015	Count2015	Mean2014	Count2014	Muutosprosentti
GREATER LONDON	947019	5418	917501	5651	3.22
WINDSOR AND MAIDENHEAD	851027	607	812147	640	4.79
SURREY	793702	5742	784542	6032	1.17
HERTFORDSHIRE	723683	3622	677935	3670	6.75
BUCKINGHAMSHIRE	700343	2778	648867	2882	7.93
BRIGHTON AND HOVE	588432	516	588451	544	-0.00
OXFORDSHIRE	562062	2584	532447	2760	5.56
WOKINGHAM	554932	1009	521463	1016	6.42
READING	539202	246	472896	280	14.02
BRACKNELL FOREST	531927	482	466597	495	14.00
WEST BERKSHIRE	527765	782	503797	831	4.76
WEST SUSSEX	491486	4233	468644	4466	4.87
BATH AND NORTH EAST SOMERSET	465811	614	490501	625	-5.03
HAMPSHIRE	463830	7483	440626	7603	5.27
POOLE	459666	1213	426080	1172	7.88

Alla on kuva Windsorin ja kaikkien alueiden hintajakaumista. Tarkastelu on tässä rajattu alle kahden miljoonan punnan vanhoihin omakotitaloihin. Jakaumat ovat positiivisesti vinoja, joten ei ole yllätys, että kummankin jakauman keskiarvo ylittää niiden mediaaniarvon.

Alla on vielä karttakuva, josta selviävät asuntojen keskihinnat alueittain. Tätä kuvaa varten on ladattava Englannin koordinaatit ja yhdistettävä ne alkuperäiseen dataan. Tämän jälkeen kuva voidaan muodostaa esim. qplotilla tai ggplotilla. Alla olevan kuvan tein qplotilla.

R toimii suurilla aineistolla suhteellisen hyvin, mutta selkeä yläraja alkaa tulla vastaan, kun datan koko ylittää 3 gigaa. Ongelmia suurten aineistojen kanssa tulee varsinkin, jos tavoitteena on mallintaa dataa esim. regressiomalleilla. Tällöin laskemisajat voivat olla jopa tunteja. R:lle on kehitetty pakkaus nimeltä SparkR, joka pyrkii parantamaan R:n toimintaa Big datan kanssa. Pakkaus yhdistää Sparkin ja R:n, eli Sparkia käytetään R:n kautta. Kirjoittelen SparkR:sta jatkossa lisää.

Ristiintaulukointi ja CLT-simulointi

Tässä tekstissä käsittelen ristiintaulukointia sekä muita hyödyllisiä datan tarkastelumenetelmiä.
Ristiintaulukointi saadaan tehtyä table-komennolla. Komento laskee kahdelle luokitellulle muuttujalle niiden väliset frekvenssit. Aloitetaan simuloimalla yritysdataa, jossa on mukana toimialaluokitus, tieto mahdollisesta konkurssista sekä nettotulosprosentti.

R-koodi:

toimiala <- sample(letters[1:4],1000,replace=TRUE)
konkurssi <- rbinom(1000,1,0.02)
nettotulos <- rnorm(500,mean=0,sd=5)+runif(500,-5,25)
data <- data.frame(toimiala,konkurssi,nettotulo

Aloitetaan datan tutkiminen toimialaluokitus-muuttujasta katsomalla eri toimialaluokkien osuudet.

table(data$toimiala)

toimiala
  a   b   c   d
241 267 244 248

Ratkaistaan myös prosenttiosuudet

prop.table(table(data$toimiala))

toimiala
    a     b     c     d
0.241 0.267 0.244 0.248

Visualisoidaan barplotilla ja piechartilla:

aa <-barplot(prop.table(table(toimiala)),col=c("5","6","2","3"),ylim=c(0,1),main="Eri toimialaluokkien osuudet kaikista")
text(aa,0.15,round(prop.table(table(toimiala)),2),cex=1.3,pos=3)

library(plotrix)
tabl <- prop.table(table(toimiala))
lbl <- paste(names(tabl), "\n", tabl, sep="")
pie3D(tabl,labels=lbl,explode=0.05,main="Eri toimialojen osuudet kaikista")

Eri toimialoista on siis lähes saman verran havaintoja.

Katsotaan seuraavaksi barplot konkurssi-muuttujasta.

a <- barplot(prop.table(table(konkurssi)),col=c("10","1"),ylim=c(0,1),main="Konkurssit")
text(a,0.2,round(prop.table(table(konkurssi)),3),cex=1.1,pos=3)
legend("topright",legend=c("Konkurssi","Ei konkurssia"),pch=c(19,19),col=c(10,1))

Konkurssiin on mennyt datan tuhannesta yrityksestä 2.5 prosenttia. Aloitetaan nyt ristiintaulukointi ja tarkastellaan konkursseja toimialoittain. Ensin table-komennolla:

table(konkurssi,toimiala)

         toimiala
konkurssi   a   b   c   d
        0 233 258 240 244
        1   8   9   4   4

prop.table(table(konkurssi,toimiala),margin=2)

        toimiala
konkurssi          a          b          c
        0 0.96680498 0.96629213 0.98360656
        1 0.03319502 0.03370787 0.01639344
         toimiala
konkurssi          d
        0 0.98387097
        1 0.01612903

Tehdään vielä barplot kuva konkurssien osuuksista eri toimialoilla

a <- barplot(prop.table(table(konkurssi,toimiala),2),beside=T,yaxt="n",ylab="%",xlab="Toimiala") # Tehdään pylväskuvio
legend("topright",legend=c("Ei konkurssia","Konkurssi"),bg="white",pch=c(19,19),col=c(1,"grey"))
axis(2,at=seq(0,1,by=0.1))
text(a, 0.1, round(prop.table(table(konkurssi,toimiala)),2),cex=1,pos=3)

Testataan prop.testilla, onko eri toimialojen välillä tilastollisesti merkittävää eroa konkurssien osuuksissa. Menetelmällä testataan nollahypoteesia, jonka mukaan osuudet ovat kaikissa toimialoissa samat. Nollahypoteesia ei kuitenkaan hylätty.

• H₀: The proportion of cases is the same in each age group: p₁ = p₂ = p₃ = p₄ = p₅ = p₆
• H1: The proportion of cases is not the same in each age group: at least one p_i is different from the others

yhteensävektori = tapply(data$konkurssi, data$toimiala, length)
konkurssejavektori = tapply(data$konkurssi[data$konkurssi==1],data$toimiala[data$konkurssi==1], length)
prop.test(konkurssejavektori,yhteensävektori)


Katsotaan seuraavaksi keskimääräiset nettotulokset eri toimialoilla ja eri konkurssi-muuttujan arvoilla.

tapply(data$nettotulos,list(data$konkurssi,data$toimiala),mean)

      a        b         c         d
0 9.963086 9.310432 10.066982  9.532678
1 4.142009 5.814025  1.901524 18.058030

Toimialoilla a,b ja c konkurssiin menneillä yrityksillä oli alempi nettulosprosentti.

z <- tapply(data$nettotulos,list(data$konkurssi,data$toimiala),mean)
a=barplot(z,beside=T,main="Nettotulosprosentti toimialan ja konkurssi-statuksen mukaan",col=c(12,14),ylim=c(0,20),args.legend=list(bg="white",x=12,y=21),legend.text=c("Ei konkurssia","Konkurssi"))
text(a,0.4, round(z,1),cex=1.3,pos=3)

Central Limit Theorem -simulointi

Lopuksi näytän vielä, kuinka Central Limit Theoremin takia voimme yleistää otoksesta saadut keskiarvot koko populaatioon riippumatta populaation jakaumasta. Saamamme luottamusvälit ovat siis valideja vaikka esim. nettotulosprosentin jakauma olisi vinoutunut. Testataan tätä R:n For-loopin avulla. Muodostetaan ensin Poisson-jakautunutta dataa, joka edustaa populaation jakaumaa ja on tarkoituksella kaukana normaalijakaumasta. Tästä populaatiosta otamme sitten otoksia ja testaamme, osuuko populaation keskiarvo otoksesta saadun keskiarvon luottamusvälille. Lasketaan otoksen keskiarvoille luottamusväli 95%:n mukaan, eli simuloinnista voimme odottaa, että n.5%  kerroista populaation keskiarvo ei osu luottamusvälillemme. Alla kuva jakaumasta.

R-koodi:

a <- rpois(100000,1)
hist(a, main="Poisson-jakauma",col="5",xlab="")
a <- rpois(100000,1)


hist(a)
a <- data.frame(a)
summary(a)


tulosvektori <- vector("numeric",10000)  #Vektori,johon loopin tulokset tulevat


se <- function(x){                  #Standard errorin funktio
sd(x)/sqrt(length(x))
}



for (i in 1:10000){                                      
set.seed(i)
d <- sample(nrow(a),0.02*nrow(a))
b <- a[d,]
tulosvektori[i] <- ifelse(mean(a[,1])>=mean(b)-1.96*se(b)& mean(a[,1]) <=mean(b)+1.96*se(b),1,0)
}

prop.table(table(tulosvektori))

tulosvektori
    0     1
0.047 0.953

Populaation oikea keskiarvo osui 10000:sta lasketusta luottamusvälistä 95.3%:iin, eli 95%:n luottamusväli piti hyvin paikkansa.

Vastaava tulos saatiin myös alla olevalla beta-jakaumalle.

Nämä tulokset ja luottamusvälin toiminta perustuvat siihen, että otoksista lasketut keskiarvot noudattavat normaalijakaumaa. Tämä voidaan nähdä katsomalla otoskeskiarvojen jakaumaa.

R-koodi:

for (i in 1:10000){
set.seed(i)
d <- sample(nrow(a),0.02*nrow(a))
b <- a[d,]
vektori[i] <- mean(b)
}

MySQL:n integroiminen R:ään

Tässä tekstissä esittelen, kuinka R:n ja MySQL:n saa yhdistettyä dbConnect-pakkauksen avulla. Yhdistämisen jälkeen voimme hakea tietoa R:llä suoraan MySQL:sta SQL-komennoin. Ensimmäinen vaihe on MySQL:n asentaminen koneelle, jonka jälkeen voidaan R:stä käsin ottaa yhteys ohjelmaan.

Alla on scripti, jolla yhteys muodostetaan. Huom. sqltk-kohdan käyttäjätiedot luodaan MySQL:n asennusvaiheessa.

R-koodi:

install.packages("dbConnect")  # Asennetaan tarvittava pakkaus

library(dbConnect)

library(RMySQL)

sqltk = dbConnect(MySQL(), user='käyttäjänimi', password='salasana', dbname='tietokannan nimi', host='host-nimi') # Haetaan MySQL:sta haluttu tietokanta

dbListTables(sqltk) # Näyttää haetussa tietokannassa olevat taulut

Seuraavaksi luon MySQL:ssa pienen tietokannan, josta voidaan hakea tietoa luodun yhteyden avulla. 

MySQL-koodi: 

create database datamyynnit;

use datamyynnit;

create table myyja(MyyjaID int(10) unsigned auto_increment primary key not null, 

Kaupunki varchar(30) not null, Puhelinnumero varchar(20) not null);

create table myynnit(MyyjaID int(10) not NULL, tilausnumero int(10) unsigned primary key not null, 

maara int(10) not null, tuote varchar(30) not null,hinta int(10) not null, pvm date not null);

insert into myyja (Kaupunki,Puhelinnumero)  values("Espoo","034343434");

insert into myyja (Kaupunki,Puhelinnumero)  values("Helsinki","042004332");

insert into myyja (Kaupunki, Puhelinnumero) values("Espoo","034363555");

insert into myyja (Kaupunki, Puhelinnumero) values("Helsinki","045545555");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("2","15","1","tuote1","900","2015-10-10");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("1","16","1","tuote2","800","2015-10-15");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("1","17","1","tuote1","900","2015-10-15");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("3","18","1","tuote2","800","2015-10-28");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("3","19","2","tuote1","900","2015-10-26");

insert into myynnit (MyyjaID,tilausnumero,maara,tuote,hinta,pvm) values("4","20","2","tuote1","900","2015-10-25");

Nyt otetaan yhteys tähän tietokantaan R:llä ja haetaan dataa SQL-komennolla.

R-koodi: 

sqltk = dbConnect(MySQL(), user='root', password='1234', dbname='datamyynnit', host='host')

myynnit <- dbGetQuery(sqltk,"select myyja.myyjaID,myyja.kaupunki,myynnit.maara,myynnit.tuote,myynnit.hinta,myynnit.pvm from myyja,myynnit where myyja.myyjaid=myynnit.myyjaid order by myyja.myyjaid")

print(myynnit)

myyjaID kaupunki maara  tuote hinta  pvm

       1      Espoo     1 tuote1   900 2015-10-15

       1      Espoo     1 tuote2   800 2015-10-15

       2  Helsinki     1 tuote1   900 2015-10-10

       3      Espoo     2 tuote1   900 2015-10-26

       3      Espoo     1 tuote2   800 2015-10-28

       4  Helsinki     2 tuote1   900 2015-10-25

Yllä on siis poimittu data. Tätä dataa voi nyt käsitellä normaalisti R:ssä, koska se on perinteisessä dataframe-muodossa. Tarkastellaan seuraavaksi dataa hieman: Luodaan uusi sarake, jossa lasketaan kullekin riville määrä*hinta. Tämä  sarake on siis myyjä-/päivä-/ ja tuotekohtainen kokonaismyynti. Esim. myyjä 1 on myynyt 15.10.2015 yhden kappaleen tuotetta 1, joten hänen sarakearvonsa on 900. 

myynnit$kokonaismyynti <-myynnit$maara*myynnit$hinta

Ratkaistaan  myyjäkohtaiset kokonaismyynnit käyttäen äsken luotua saraketta.

tapply(myynnit$kokonaismyynti,myynnit$myyjaID,sum)

  1    2    3    4 

1700  900 2600 1800 

Ratkaistaan myyjäkohtaiset kokonaismyynnit tuotteelle 1

tapply(myynnit$kokonaismyynti[myynnit$tuote=="tuote1"],myynnit$myyjaID[myynnit$tuote=="tuote1"],sum)

  1    2    3    4 

 900  900 1800 1800 

Ratkaistaan myyjäkohtaiset tuotteen 1 kokonaismyynnit Espoossa ennen 28.10.2015

myynnit$pvm <- as.Date(myynnit$pvm,format="%Y-%m-%d")

tapply(myynnit$kokonaismyynti[myynnit$tuote=="tuote1" & myynnit$kaupunki=="Espoo"& myynnit$pvm <= "2015-10-28" ],myynnit$myyjaID[myynnit$tuote=="tuote1" & myynnit$kaupunki=="Espoo "& myynnit$pvm <= "2015-10-28" ],sum)

 1    3 

 900 1800 

Katsotaan lopuksi vielä myyjien kokonaismyynnit jokaiselle datassa esiintyvälle päivälle

for (i in factor(unique(myynnit$pvm))){

print(paste(c("Myyjä1","Myyjä2","Myyjä3","Myyjä4"),"pvm",i,tapply(myynnit$kokonaismyynti[myynnit$pvm==i],myynnit$myyja[myynnit$pvm==i],sum),"€"

))

}

"Myyjä1 pvm 2015-10-10 NA €"  "Myyjä2 pvm 2015-10-10 900 €"
[3] "Myyjä3 pvm 2015-10-10 NA €"  "Myyjä4 pvm 2015-10-10 NA €" 
[1] "Myyjä1 pvm 2015-10-15 1700 €" "Myyjä2 pvm 2015-10-15 NA €"  
[3] "Myyjä3 pvm 2015-10-15 NA €"   "Myyjä4 pvm 2015-10-15 NA €"  
[1] "Myyjä1 pvm 2015-10-25 NA €"   "Myyjä2 pvm 2015-10-25 NA €"  
[3] "Myyjä3 pvm 2015-10-25 NA €"   "Myyjä4 pvm 2015-10-25 1800 €"
[1] "Myyjä1 pvm 2015-10-26 NA €"   "Myyjä2 pvm 2015-10-26 NA €"  
[3] "Myyjä3 pvm 2015-10-26 1800 €" "Myyjä4 pvm 2015-10-26 NA €"  
[1] "Myyjä1 pvm 2015-10-28 NA €"  "Myyjä2 pvm 2015-10-28 NA €" 
[3] "Myyjä3 pvm 2015-10-28 800 €" "Myyjä4 pvm 2015-10-28 NA €" 

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.