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