Pähkinä - numeroiden lukumäärät
Sain jälleen kerran ratkottavakseni pähkinän, joka kuuluu jokseenkin seuraavasti: “Kyltintekijän pitää numeroida kadun varren talot numeroilla 1-100. Kuinka monta yhdeksikköä hän tarvitsee?” Kyseisellä tavalla muotoiltuna pähkinä on riittävän helppo käsipelillä ratkottavaksi, mutta entäs jos onkin tarpeen laskea kaikkien tarvittavien numeroiden määrät tai numeroimiseen käytettävät numerot muodostavat laajemman sarjan, esimerkiksi kaikki luvut yhdestä miljoonaan?
Ongelman voi ratkoa nopeahkosti monillakin ohjelmointikielillä, mutta tässä eräitä “yksirivisiä” ehdotelmia R:ää käyttäen.
Yhdeksikköjen lukumäärä:
sum(unlist(strsplit(as.character(1:100), ""))==9)
Kaikkien lukujen lukumäärät:
table(unlist(strsplit(as.character(1:100), "")))
Kirjainten lukumäärät, esimerkiksi laatikoiden nimeämiseen:
teksti<-c("puutarhatontut", "pajavasarat", "valokuvat") table(unlist(strsplit(teksti, "")))
Riskin visualisointi
Erilaisten terveysriskien hahmottaminen voi välillä olla vaikeaa. Useinhan riskeistä puhutaan esimerkiksi prosentteina tyyliin: “Sinulla on 1% mahdollisuus kuolla ennen 30. syntymäpäivääsi”. Yhden prosentin riski voi olla helpommin käsitettävissä, jos siitä puhutaan luonnollisilla frekvensseillä eli tyyliin “1 henkilö sadasta kuolee ennen 30. syntymäpäiväänsä”. Luonnollisten taajuuksien visualisointi on myös helppoa.
Kirjan The Illusion of Certainty kirjoittajat ehdottavat, että riskien havainnollistamiseen käytettäisiin Risk Characterization Theatre -esitystapaa. Siinä teatterin istumakaaviota käyttäen voidaan eri väreillä merkitä esimerkiksi sairastumisen riskissä olevien henkilöiden määrää osuutena kaikista henkilöistä. Stubborn Mule on tehnyt näistä R-pohjaisen visualisoinnin, jolla hän havainnollistaa kirjassa esitettyä esimerkkiä ja tupakoinnin aiheuttamaa riskiä.
Tuossa teatterimuotoisessa esityksessä minua häiritsi suuresti sen muoto ja ihmisjoukon suuruus. Kuviossa on nimittäin tuhat henkilöä, mutta kuinka moni on nähnyt jossakin salissa 1 000 yhdellä kertaa? Lisäksi kuvio on useampien riskien yhtäaikaista visualisointia ajatellen aavistuksen hankalan muotoinen. Niinpä kirjoitin oman funktion, jossa henkilöiden määrä perustuu Finnairinkin käyttämän Boeing 757-200 -koneen istumapaikkakaavioon. Koneeseen mahtuu 227 matkustajaa.
Itse funktio rct() näyttää R-koodina seuraavalta:
rct<-function(plans="Finnair Boeing 757-200", cases, main) { if(plans=="Finnair Boeing 757-200") { seats<-227 available<-"#009900" occupied<-"grey75" m<-matrix(nrow=6, ncol=39, data=available) m[1:3,1]<-"white" m[1,10]<-"white" m[1,11]<-"white" m[6,11]<-"white" m[6,30]<-"white" x11(width=7, height=0.89) par(mar=c(0.1,2,2,0.1), xpd=NA) plot(x=1, y=1, xlim=c(1, ncol(m)), ylim=c(1, nrow(m)), col="white", type="n", axes=F, xlab="", ylab="", main=main, family="sans") lines(x=c(0.75,-2.5), y=c(0.5,3)) lines(x=c(-2.5,-2.5), y=c(3,4)) lines(x=c(-2.5,0.75), y=c(4,6.5)) lines(x=c(39.5, 40.5), y=c(0.5,2.5)) lines(x=c(40.5, 40.5), y=c(2.5,4.5)) lines(x=c(40.5,39.5), y=c(4.5,6.5)) } if(cases>seats) { stop(paste("The plane can only accommodate ", seats, " passangers. You have asked for more cases than available seats!", sep="")) } m2<-m m2[m2!="white"]<-1:seats m[m2 %in% as.character(sample(1:seats, cases))]<-occupied for(i in 1:nrow(m)) { for(j in 1:ncol(m)) { points(x=j, y=i, pch=15, col=m[i,j]) } } }
Funktio ottaa parametreina kuvioon harmaalla merkittävien kuolleiden henkilöiden lukumäärän, joka voi siis olla korkeintaan 227. Lisäksi kuviolle voi antaa otsikon. Esimerkiksi komento:
rct(cases=100, main="Testi")
piirtää kuvion, jossa 100 penkkiä on väritetty harmaalla (kuolleet), ja jäljelle jäävät 127 penkkiä (elossa olevat) vihreällä.
Miltä esimerkiksi Suomen vuoden 2010 kuolleisluvut (Tilastokeskus) näyttävät tätä visualisointi käyttäen? Tämä on esitetty seuraavassa kuviossa, joka siis vastaa kysymykseen: “Jos vauvoja syntyi yhteensä 227, kuinka moni heistä on elossa täytettyään X vuotta?”
Saman tiedon voi toki visualisoida paljon pienemmässä tilassa esimerkiksi viiva- tai pylväskaaviota käyttäen, mutta tällainen esitystapa voi toisinaan olla vaikuttavampi.
(g)lmer on epäluotettava?
Zhang et al. ovat julkaisseet vastikään artikkelin, jossa he vertaavat eräiden tilasto-ohjelmistojen sekamallien sovittamiseen kehitettyjä ominaisuuksia. Artikkelissa esitetyt tulokset eivät ylipäätään ole kovin mairittelevia, sillä parhaiksi havaittujenkin ohjelmistojen (SAS:n NLMIXED) tyypin I -virhetaso on noin kaksinkertainen ilmoitettuun verrattuna. Erityisen huonosti tuntuu kuitenkin pärjäävän R:n laajennuspaketin lme4 funktion glmer(). Lisäksi R sai pyyhkeitä myös epäluotettavasta mallin paramatrien arvioinnista.
Yritetäänpä toistaa artikkelin tulokset ajantasaisella 32-bittisellä R:llä (R 2.14.2., lme4 0.999902344-0, Matrix 1.0-4, splines 2.14.2). Asiasta on keskusteltu myös ainakin kahdessa viestiketjussa (1, 2) R-sig-ME -listalla.
1. Davisin aineiston analyysi
Davisin aineiston analyysi voidaan toistaa R:ssä seuraavalla koodilla:
# HSAUR contains respiratory data data("respiratory", package = "HSAUR") # Correct an error in the data respiratory$status[rownames(respiratory)=="428"]<-"poor" # Convert pretest to covariate "baseline" resp <- subset(respiratory, month > "0" ) resp$baseline <- rep(subset(respiratory, month == "0")$status, each=4) # Invert the factor levels to correspond to Davis et al. resp$treatment2<-factor(resp$treatment, levels=c("treatment","placebo")) resp$sex2<-factor(resp$sex, levels=c("male", "female")) # Fit the models library(lme4) # Laplace approximation print(m_glmer_4.L <- glmer(status ~ centre + treatment2 + sex2 + age + baseline + (1|subject), family=binomial,data=resp), cor=FALSE) # Gauss-Hermite approximation print(m_glmer_4.L <- glmer(status ~ centre + treatment2 + sex2 + age + baseline + (1|subject), family=binomial,data=resp, nAGQ=16), cor=FALSE)
Artikkelissa esitetyt SAS:n NLMIXED-ohjelmalla saadut mallin parametrien estimaattein arvot ja tässä R:llä lasketut ovat hyvin lähellä toisiaan, vakiotermiä lukuunottamatta. Tässä lasketut R:llä tuotetut parametrien estimaattien arvot eroavat selvästi artikkelissa esitetyistä, eikä ole vielä selvillä, miten ne oikeastaan on laskettu. On toki myönnettävä, että esimerkiksi iän vaikutuksen estimaatti on R:llä laskettuna noin 33% pienempi kuin SAS:n NLMIXED-ohjelmalla laskettu estimaatti, joten suhteellinen virhe voi pienillä estimaatin arvoilla olla suurikin.
SAS R Intercept 2.25 0.70 Center 0.97 0.95 Treatment -2.10 -2.18 Sex -0.24 -0.24 Age -0.03 -0.02 Baseline 2.94 2.94
R:ssä funktio glmer() tuottaa hyvin samanlaiset estimaatit, käyttipä mallin sovittamiseen Laplace- tai Gauss-Hermian -menetelmää.
2. Simulaatiot
Otetaan tähän esimerkiksi yksi simulaatio, jossa mallin molempien parametrien estimaattien (Beta) pitäisi olla 1, koska simulaatio on toteutettu siten. Oletetaan kovarianssimatriisin (Cor) luomisessa käytettävän parametrin (Tau) arvoksi simulaatioissa 2, ja ryhmäkohtaisen havaintojen (N) lukumääräksi 500. Toistetaan aineiston simulointi ja mallin sovitus 1000 kertaa (R), ja vedetään toistoista saadut tulokset lopuksi yhteen. Simulaatioissa käytetty koodi on seuraava (muokattu aiemmin mainituista R-sig-ME:n viestiketjuista, erityisesti Ben Bolker ja Dave Fournier):
# Parameters of the simulation N<-500 Beta<-c(1,1) Tau<-2 Cor<-matrix(c(1,0.25,0.25,1),nrow=2) R<-1000 # Simulation function for the data library("lme4") simfun <- function(n,beta=c(1,1),tau=2.0, Cor=matrix(c(1,0.25,0.25,1),nrow=2)) { D <- expand.grid(id=factor(1:n),t=1:3) D$x <- rnorm(3*n) # x_{it} ~ N(0,1) X <- model.matrix(~x,data=D) Z <- model.matrix(~id-1+id:x,data=D) b <- MASS::mvrnorm(n,mu=c(0,0),Sigma=tau^2*Cor) D$Y <- rbinom(nrow(D), prob=plogis(X %*% beta + Z %*% c(b)), size=1) D } # Simulations m1<-matrix(ncol=2, nrow=R, data=NA) for(i in 1:R) { set.seed(i) dat=simfun(N) fm <- glmer(Y~ x + ( x|id), data=dat,family=binomial) m1[i,] <- cbind(as.vector(fixef(fm))[2], as.vector(sqrt(vcov(fm,useScale = FALSE)[2,2]))) } # Calculate the confidence intervals based on normal approximation cilower<-m1[,1]-m1[,2]*1.96 ciupper<-m1[,1]+m1[,2]*1.96 # How many times is the known parameter outside of the confidence interval? sum(cilower>1 | ciupper<1)/R
Artikkelissa SAS:n NLMIXED-ohjelmalle ilmoitettu tyypin I-virheen frekvenssi oli noin 5-10%, ja liki samaan tulokseen on ADMB-työkalun osalta päätynyt Dave Fournier. Yllä suoritetun simulaation perusteella R:n glmer()-funktion tyypin I-virheen frekvenssi sen sijaan todella näyttää olevan paljon korkeampi, noin 25%.
Mitä tästä sitten pitäisi päätellä? Ainakin se voitaneen varmuudella sanoa, että glmer()-funktion antamat luottamusvälit ja p-arvot ovat liian optimistisia, mikä johtaa paljon oletettua korkeampaan väärien positiivisten tulosten määrään. Funktion glmer() antamien p-arvojen tulkinnassa kannattanee siis olla hieman varovainen.
Muste-workshop
Survo-käyttäjien yhdistys järjestää huhtikuussa nk. Muste-workshopin, jonka tarkoituksena on kerätä yhteen R- ja Survo-asiantuntijat keskustelemaan Muste-projektista. Ohessa lainaus workshopin alustuksesta, josta käy ilmi tapaamisen tarkoitus:
Muste on R-ohjelmistolle kehitetty, avoimeen lähdekoodiin perustuva ja vapaasti
saataville tuleva toteutus professori Seppo Mustosen elämäntyönään luomasta Survoohjelmistosta, jonka kehityskaari on jo peräti 50 vuoden mittainen.
Musteen päävastuullinen kehittäjä on Reijo Sund, joka 2009 esitti idean Musteesta ja
sen toteuttamisesta R-ohjelmiston puitteissa. R on alun perin Robert Gentlemanin ja
Ross Ihakan alullepanema toteutus professori John Chambersin johdolla kehitetystä
tilastollisesta S-ohjelmointikielestä. Nykyisin R on laajan kansainvälisen käyttäjäyhteisön kehittämä, vapaasti saatava ohjelmisto, joka on noussut 2000-luvulla suosioon etenkin tilastotieteen, matematiikan ja tietojenkäsittelytieteen piirissä, mutta myös yhä useammalla näitä menetelmätieteitä hyödyntävillä aloilla.
Muste liittää Survon pitkän kehityshistorian aikana syntyneet ja hioutuneet ominaisuudet ja innovaatiot aidoksi osaksi R-ohjelmistoa. Muste on yksi ominaisuuksiltaan monipuolisimmista R:n laajennuspaketeista. Se tarjoaa editoriaalisen käyttöliittymän lisäksi lukuisia vaihtoehtoisia työtapoja R-käyttäjille mm. aineistojen hallinnointiin, esikäsittelyyn ja muokkaukseen, matriisilaskentaan, kuvien piirtoon ja raportointiin.
Juuri nyt ollaan erittäin mielenkiintoisessa kehitysvaiheessa, sillä Muste toimii jo
luotettavasti monenlaisissa ympäristöissä ja sisältää yli 90 % Survon toiminnoista.
Erilaisten luovien käyttötapojen kekseliäät yhdistelmät hakevat ja löytävät uusia uomia.
Lisätietoja asiasta löytyy workshop-ilmoituksesta.
Tyhjien korvaaminen edellisellä arvolla
R on täynnä toinen toistaan hyödyllisempia valmiita funktioita. Usein on tilanne, että data ei ole rakenteeltaan eheää ja sitä joutuu käsittelemään paljonkin ennen varsinaista analyysiä. Törmäsin tänään tilanteeseen, jossa datassa oli mukana väliotsikot yhdessä sarakkeessa ja väliotsikoiden alla tyhjää (NA). Tavoitteenani oli täyttää nämä tyhjät rivit väliotsikoilla, ja tämän jälkeen suodattaa datasta pois tietyt rivit väliotsikoiden perusteella. Mietin, että miten tämän asian voisin toteuttaa, kunnes keksin, että joku on varmasti tämän jo miettinyt valmiiksi. Nopea googlaus, ja vastaus oli löytynyt.
Ohessa esimerkki:
library(zoo) DF <- data.frame(var1 = c("a", NA, NA, NA, "b", NA, NA, NA, "c", NA), var2 = 1:10)
> DF var1 var2 1 a 1 2 <NA> 2 3 <NA> 3 4 <NA> 4 5 b 5 6 <NA> 6 7 <NA> 7 8 <NA> 8 9 c 9 10 <NA> 10
DF$var1 <- na.locf(DF$var1)
>DF var1 var2 1 a 1 2 a 2 3 a 3 4 a 4 5 b 5 6 b 6 7 b 7 8 b 8 9 c 9 10 c 10 >
Suomen sisäinen muuttoliike
Louhos-blogissa oli taannoin visualisointi suomalaisten maailmanlaajuisista muuttoliikkeistä. Samantyyppisen kuvion voi laatia tarkemmalla tasolla Suomen sisälläkin. Sopivaa aineistoa löytyy suoraan Tilastokeskuksen muuttoliike-tilastoista ja kuntien koordinaatit saa selvitettyä, ainakin sinnepäin, Wikipedian sivulta. Tilastokeskuksen tarjoamassa aineistossa on vain se ongelma, ettei kaikkien kuntien välisiä muuttoliikkeitä ole mahdollista hakea yhtä aikaa palvelun asettaman rivirajoituksen vuoksi. Niinpä seuraavassa onkin esimerkki Kangasalan tulo- ja lähtömuutosta vuosina 1987-2010.
Päätin kuvata muuttolikkeen kartalla viivoilla, joiden väri kertoo muuton suunnan ja paksuus muuttajien lukumäärän. Karttapohjan hakemiseen käytettiin sorvi-pakettia.
Seuraava R-koodi tuottaa alla olevan kaltaisen kuvan. Koodi olettaa tarvittavien tietojen löytyvän datakehikosta nimeltä d2. Sen tuottamiseen tarvittava R-koodi on kuvan alapuolella.
library(sorvi) gadm<-GetGADM(resolution = "FIN_adm", taso = 4) plot(gadm, col="black", border="grey25", bg="black", lwd=0.5) for(i in 1:nrow(d2)) { lines(y=c(d2$lat.tulo[i], 61.463), x=c(d2$lon.tulo[i], 24.076), lwd=log10(abs(d2$nettomuutto[i]))*3, col=d2$col[i]) } title("Kangasala - nettomuutto") legend(x="topright", lwd=c(6), col=c("#0066CC", "#CC0000"), lty=1, legend=c("Tulomuutto","Lähtömuutto"), text.col="grey75") legend(x="bottomright", lwd=c(1,2,3)*3, col="grey75", lty=1, legend=c("10", "100", "1000"), title="Henkilöä", text.col="grey75")
Suurin muuttovoitto Kangasalle suuntautuu Tampereelta. Suurin muuttotappio puolestaan suuntautuu lähialueille, kuten Nokialle, Pälkäneelle ja Orivedelle, sekä suuriin kaupunkeihin, kuten Helsinkiin ja Turkuun.
d2<-structure(list(lat.tulo = c(60.633, 63.923, 61.35, 60.801, 63.1, 63, 63.752, 60.803, 62.3, 61.983, 61.316, 60.633, 62.183, 63.544, 64.083, 65.317, 60.867, 61.567, 65.967, 64.677, 63.1, 65.817, 60.713, 61.067, 60.417, 65.733, 60.46, 64.125, 61.433, 60.8, 66.5, 60.819, 62.9, 60.117, 61.867, 64.233, 60.407, 63.844, 60.851, 62.27, 61.633, 60.285, 61.483, 60.944, 60.4, 62.317, 61.772, 60.209, 65.017, 60.981, 61.317, 62.033, 61.169, 61, 60.459, 61.317, 61.267, 61.467, 61.677, 60.17, 61.333, 61.5), lon.tulo = c(25.317, 24.958, 25.283, 21.411, 21.6, 23.817, 25.329, 23.483, 27.133, 24.083, 22.135, 24.876, 27.833, 29.145, 24.55, 25.367, 26.7, 25.183, 29.183, 24.459, 23.083, 24.533, 24.438, 24.383, 24.333, 24.567, 26.946, 29.511, 21.883, 25.736, 25.717, 23.621, 27.683, 24.433, 25.2, 27.683, 25.036, 23.128, 23.043, 23.7, 23.2, 25.05, 21.783, 25.937, 25.117, 27.917, 23.067, 24.658, 25.467, 25.655, 23.617, 24.633, 23.858, 24.45, 22.25, 23.75, 24.033, 23.5, 24.363, 24.931, 24.28, 23.75), nettomuutto = c(11L, 11L, 11L, -11L, 11L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, -13L, 13L, 13L, -14L, 14L, 14L, -14L, 14L, 15L, 16L, 16L, -17L, -17L, 18L, -18L, 18L, 18L, -20L, 20L, 21L, 22L, 25L, 26L, 26L, 26L, 27L, 27L, 28L, -29L, 30L, 31L, 33L, 34L, 35L, -36L, 43L, 54L, 56L, -71L, 85L, -104L, -110L, -118L, -136L, -141L, -166L, -211L, -219L, -373L, 4000L ), col = c("#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#CC000066", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#0066CC66", "#0066CC66", "#CC000066", "#0066CC66", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#CC000066", "#0066CC66")), .Names = c("lat.tulo", "lon.tulo", "nettomuutto", "col"), class = "data.frame", row.names = c(177L, 184L, 197L, 301L, 303L, 2L, 20L, 60L, 207L, 245L, 28L, 42L, 62L, 187L, 322L, 46L, 116L, 120L, 126L, 231L, 88L, 94L, 153L, 30L, 314L, 92L, 115L, 119L, 295L, 190L, 243L, 17L, 122L, 104L, 70L, 74L, 292L, 110L, 152L, 318L, 44L, 306L, 214L, 182L, 96L, 307L, 49L, 12L, 194L, 132L, 310L, 178L, 1L, 45L, 290L, 144L, 304L, 185L, 192L, 35L, 229L, 281L))
Millimetripaperia automaattisesti
R-Help -postilistalla oli vastikään mielenkiintoinen kysymys millimetripaperin tuottamisesta R:llä. Tähän tarvitaan laajennuspaketti grid, mutta sen asentamisen jälkeen työ käy helposti (koodi kopioitu suoraan ylle linkatusta postista):
library(grid) pdf("grid.pdf", width=21/2.54,height=29.7/2.54) grid.grill(h = unit(seq(0, 297, by=1), "mm"), v = unit(seq(0, 210, by=1), "mm"), gp=gpar(col="grey",lwd=0.1)) grid.grill(h = unit(seq(0, 297, by=5), "mm"), v = unit(seq(0, 210, by=5), "mm"), gp=gpar(col="grey20",lwd=0.2)) grid.grill(h = unit(seq(0, 29.7, by=1), "cm"), v = unit(seq(0, 21, by=1), "cm"), gp=gpar(col="black",lwd=1)) dev.off()
Tuloksena on seuraavanlainen ruudukko, jonka voi sitten tulostaa paperille:
Samassa keskustelussa mainittiin myös komento grid(), joka tuottaa automaattisesti vaalean harmaan ruudukon esimerkiksi hajontakuvion taa. Perinteinen tapa on ollut käyttää komentoa abline() esimerkiksi seuraavasti:
plot(0,0,type="n") abline(h=c(-1,0,1), col="grey75") abline(v=c(-1,0,1), col="grey75") points(x=seq(-1,1,0.05), rnorm(41), pch=16)
Komentoa grid() käyttäen sama koodi yksinkertaistuu muotoon:
plot(0,0,type="n") grid() points(x=seq(-1,1,0.05), rnorm(41), pch=16)
Alla olevassa kuvassa on esitetty molempien tapojen tuottamat hajontakuviot alekkain:
R Inferno - dropping dimensions
Patrick Burnsin kirja R Inferno esittelee erittäin suuren määrän mahdollisia sudenkuoppia, jotka liittyvät R-kielen käyttöön. Kirja on monellakin tapaa hyvin laadittu, mutta omien taitojen etenemistä voi hyvin mitata sen avulla, mihin kirjan lukuun useimmat kohdatut ongelmat sijoittuvat, sillä kirja noudattelee kevyesti Danten Inferno-runoa.
“If you are using R and you think you’re in hell, this is a map for you.”
Taulukkomuotoisten aineistojen, esimerkiksi datakehikoiden käsittelyyn liittyvä kohtuullisen usein vastaantuleva knoppi (R Inferno, 8.1.44) on, että kun datakehikosta poimitaan yksi rivi tai sarake, ei tuloksena saatava objekti enää olekaan datakehikko, vaan vektori. Esimerkiksi:
xdf2 <- data.frame(a=1:4, b=42:45, row.names=LETTERS[1:4]) xdf2[,1] [1] 1 2 3 4 class(xdf2[,1]) [1] "integer"
Ongelma voidaan korjata lisäämällä poiminnan hoitavan alaindeksin sisään argumentti drop=FALSE:
xdf2 <- data.frame(a=1:4, b=42:45, row.names=LETTERS[1:4]) xdf2[,1, drop=FALSE] a A 1 B 2 C 3 D 4 class(xdf2[,1, drop=FALSE]) [1] "data.frame"
Lisähyötynä argumentin drop=FALSE käytöstä on, että sarakkeiden ja rivien nimet saadaan mukaan muodostuvaan objektiin, sillä ne katoavat, jos tuloksena on vektori (vrt. yllä), ja ne pitää mahdollisesti lisätä jälkikäteen erikseen.
Labyrinttien generointi
Etsin R-ohjelmointikurssille ideoita, ja löysin Buckblogs-blogin, joka luettelee suuren määrän erilaisia labyrinttien generointiin soveltuvia algoritmeja. Labyrintillä tarkoitan siis sellaista sokkelotehtävää, jossa pyritään löytämään reitti alkupisteestä johonkin kohteeseen, esimerkiksi labyrintin uloskäynnille.
Robert McGehee näyttää jo ohjelmoineen R:llä pari hyvin toimivaa funktiota, joilla voidaan tuottaa halutun kokoisia, neliön muotoisia labyrinttejä. Hänen ideoistaan laajentamalla tai keskeisiä algoritmin osia korvaamalla voisi hyvinkin saada toimivan harjoitustehtävän, tai ainakin ajanvietettä joulun pyhiksi.
Tilastolliset tietosuojamenetelmät R:ssä
Tilastollisten tietosuojamenetelmien tarkoituksena on estää esimerkiksi henkilö-, yritys- tai muun yksikkötasoisen tiedon paljastuminen. Paljastuminen voi koskea esimerkiksi yksikön identiteettiä (suora tunnistaminen) tai ominaisuuksia (epäsuora tunnistaminen). Epäsuora tunnistus voi lisäinformaation kanssa yhdistettynä johtaa myös yksikön suoraan tunnistamiseen. Tunnistus voi tapahtua jopa hyvin vähäisiä tietomääriä käyttäen (esimerkiksi syntymävuosi ja sukunimi), ja voi perustua myös esimerkiksi taulukoituihin tietoihin. Aineisto voidaan suojata kahdella tavalla, joko naamioimalla alkuperäistä aineistoa sopivasti tai luomalla sen perusteella täysin synteettinen, simuloitu aineisto. Tilanne voi tulla vastaan esimerkiksi viranomaistoiminnassa tai epidemiologisessa tutkimuksessa, jos joku projektin ulkopuolinen tutkija haluaisi saada aineiston analysoitavakseen vaikkapa artikkelin tulosten toistamiseksi.
Tietosuojan tilastolliseen hallintaan on kehitetty useita tietokoneohjelmia, esimerkiksi micro-Argus ja tau-Argus, mutta R:stäkin löytyy työkaluja, joista tärkeimpiä ovat laajennuspaketit sdcMicro ja sdcTable. Paketti sdcMicro on tarkoitettu yksikötasoisen tiedon suojaamiseen ja sdcTable puolestaan taulukkomuotoisille aineistoille. Seuraavassa havainnollistetaan vain paketin sdcMicro toimintaan, sillä useimmiten tutkijoiden käytössä on jonkinlaista yksikkötason tietoa, joka voi olla tarpeen suojata. Paketin mukana tulee micro-Arguksestakin löytyvä aineisto free1, jota esimerkissä käytetään.
# Ladataan data data(free1) # Perusta, yhden tai kahden yksikön ryppäiden lukumäärät # Komennossa pitää määrittää ne muuttujat, jotka katsotaan keskeisiksi # Tässä ne ovat sarakkeissa 1-3 (alue, sukupuoli ja ikä) f <- freqCalc(free1, keyVars=1:3, w=30) f -------------------------- 2664 observation with fk=1 788 observation with fk=2 -------------------------- # Lasketaan yksiköiden tunnistustodennäköisyyksiä ind <- indivRisk(f) ind method=approx, qual=1 --------------------------- 534 obs. with much higher risk than the main part # Suojataan aineisto siten, ettei yksittäisiä yksiköitä # voi enää tunnistaa. Laskenta voi kestää pitkään! # Parametrin k arvolla määritetään, kuinka monta yksikköä # kussakin suojatun aineiston ryhmässä saa olla. free2<-localSupp2(free1, keyVars=1:3, w=30, k=3, method="fastSupp") [1] "estimated time for method minimizeSupp: 2131 (in seconds)" [1] "estimated time for method fastSupp: 1385 (in seconds)" [1] "3-anonymity is not achieved. Re-apply localSupp2 on this result or use function localSupp2Wrapper " # Täydellistä anonymiteettiä ei ilmoituksen mukaansaavutetty, # mutta vilkaistaan tulos silti f2<-freqCalc(free2$xAnon, keyVars=1:3, w=30) f2 -------------------------- 0 observation with fk=1 0 observation with fk=2 -------------------------- # Myös kunkin ryhmän havaintojen määrä tukee tulosta min(f2$fk) [1] 21
Pienehköillä aineistoilla paikallinen suppressointi on vielä mahdollista, mutta suuria, kenties satoja tuhansia havaintoja sisältävien aineistojen käsittely tulee nopeasti mahdottomaksi, etenkin jos anonymisoinnin kannalta keskeisiä muuttujia on paljon. Laajennuspaketista löytyy paikallisen suppression lisäksi myös muita suojaamiseen soveltuvia menetelmiä, kuten mikroaggregointi, havaintoarvojen vaihto havaintojen kesken, kohinan lisääminen ja muuttujien arvojen uudelleenkoodaus, ja osa niistä toimii huomattavasti nopeammin kuin tässä esitelty paikallinen suppressointimenetelmä.
Aineiston suojaaminen on usein pitkällinen prosessi, jos sen haluaa tehdä hyvin, ja vaatii tyypillisesti useampien menetelmien käyttämistä yhtäaikaisesti ja kenties muuttujakohtaisesti sopivimman menetelmän valitsemista.