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.
Datan lukeminen R:ään: Scan
Blogin päivittämisessä on ollut kohtalaisen pitkä tauko muiden kiireiden johdosta. Ajattelin, että seuraava koodinpätkä voisi olla jakamisen arvoinen, se osoittaa mielestäni hyvin R:n ilmaisuvoiman. Esimerkin lähtökohtana on data, jossa koko data on tekstitiedostossa yhdellä rivillä. Datan rakenne on seuraava:
Henkilö1 Hetu Henkilö2 Hetu…
Seuraava koodi muuttaa datan data.frame-muotoon:
x <- scan(file = textConnection("Henkilo1 XXXXXX-XXXX Henkilo2 XXXXXX-XXXX Henkilo3 XXXXXX-XXXX Henkilo4 XXXXXX-XXXX"), what = list(Henkilö="character", Hetu="character")) as.data.frame(x)
Ajettuna
> x <- scan(file = textConnection("Henkilo1 XXXXXX-XXXX Henkilo2 XXXXXX-XXXX Henkilo3 XXXXXX-XXXX Henkilo4 XXXXXX-XXXX"), what = list(Henkilö="character", Hetu="character")) Read 4 records > as.data.frame(x) Henkilö Hetu 1 Henkilo1 XXXXXX-XXXX 2 Henkilo2 XXXXXX-XXXX 3 Henkilo3 XXXXXX-XXXX 4 Henkilo4 XXXXXX-XXXX >
Komento order()
Komento order() lajittelee havainnot nousevaan tai laskevaan järjestykseen. Esimerkiksi seuravassa muodostetettavan taulukon d rivit saadaan järjestetty muuttujan x suhteen numeeriseen järjestykseen yhdistämällä havaintojen järjestäminen ja rivien poiminta alaindeksillä:
d<-data.frame(x = c(1,1,3:1,1:4,3), y = c(9,9:1), z = c(2,1:9) ) d[order(d$x),] x y z 1 1 9 2 2 1 9 1 5 1 6 4 6 1 5 5 4 2 7 3 7 2 4 6 3 3 8 2 8 3 3 7 10 3 1 9 9 4 2 8
Argumentilla decreasing voi määrittää, haluaako rivit nousevaan vai laskevaan järjestykseen:
> d[order(d$x, decreasing=TRUE),] x y z 9 4 2 8 3 3 8 2 8 3 3 7 10 3 1 9 4 2 7 3 7 2 4 6 1 1 9 2 2 1 9 1 5 1 6 4 6 1 5 5
Lajittelusuunta voidaan määritellä myös yksittäisille sarakkeille sijoittamalla niiden eteen miinus-merkki tai jättämällä se pois esimerkiksi seuraavasti:
d[order(d$x, d$y),] x y z 6 1 5 5 5 1 6 4 1 1 9 2 2 1 9 1 7 2 4 6 4 2 7 3 10 3 1 9 8 3 3 7 3 3 8 2 9 4 2 8 d[order(d$x, -d$y),] x y z 1 1 9 2 2 1 9 1 5 1 6 4 6 1 5 5 4 2 7 3 7 2 4 6 3 3 8 2 8 3 3 7 10 3 1 9 9 4 2 8
Jostakin syystä tuo miinus-merkin käyttö oli aiemmin jäänyt minulta huomaamatta, ja olekin aina lajitellut useampien sarakkeiden perusteella Excel-ohjelmaa käyttäen, jos sellaiseen on ollut tarvetta. Aina sitä oppii uutta!
Logistisen regression selitysaste
Tein erääseen julkaisuun logistisia regressiomalleja, ja käsikirjoituksen vertaisarvioija halusi nähdä mallien max rescaled R squared [sic] -arvon. Onneksi on Google, joka tiesi kertoa, että se on SAS:n logistiselle regressiolle ilmoittama pseudo-selitysaste, mutta sen laskeminen R:ssä olikin toinen juttu. Laajennuspaketista pscl löytyi nopeasti komento pR2(), joka laskee erilaisia pseudo-selitysasteita:
n <- 100 x <- c(rnorm(n), 1+rnorm(n)) y <- c(rep(0,n), rep(1,n)) fit <- glm(y~x, family=binomial) library(pscl) pR2(fit) # llh llhNull G2 McFadden r2ML r2CU # -121.2974283 -138.6294361 34.6640157 0.1250240 0.1591316 0.2121754
Mutta mikä pR2():n raportoimista arvoista on sama kuin SAS:n Max-rescaled Rsquared? Puoli päivää asiaa tutkittuani sain viimein selvitettyä, että yllä olevassa tulosteessa r2ML on Cox-Snell:in selitysaste, joka vastaa SAS:ssa R Squared -arvoa. Arvo r2CU on Cragg-Uhlerin selitysaste, jota myös Nagelkerken selitysasteeksi kutsutaan, ja joka vastaa tuota etsimääni SAS:n max-rescaled selitysastetta. Mainittakoon täydellisyyden vuoksi, että tuo arvo McFadden on McFaddenin selitysaste, jota toisinaan kutsutaan myös nimellä deviance pseudo-R squared.
Olisin väittänyt referee:lle vastaan koko selitysasteen käytöstä yleistettyjen lineaaristen mallien tapauksessa, mutta palautteessa oli niin monta muutakin taistelua vaativaa seikkaa, että jossain piti antaa periksikin.
lmer() ja p-arvot
R:ään on saatavilla useampiakin laajennuspaketteja, joilla voidaan sovittaa aineistoon lineaarisia sekamalleja. Eräs ensimmäisiä paketteja oli nlme, josta on sittemmin jatkokehitetty paketti lme4. Monia käyttäjiä harmittava piirre on, ettei lme4-paketin funktio lmer() anna mallin parametrien estimaateille p-arvoja. Syynä on, ettei paketin kehittäjä Douglas Bates luota p-arvojen laskemiseksi tarvittavien vapausasteiden nykyisiin määrittämismenetelmiin, eikä siten ole halunnut tuoda ko. toiminnallisuutta pakettiinsa. Perusteltua, mutta siltikin käyttäjän kannalta hankalaa.
SAS-ohjelmisto käyttää sekamallien F-testiin perustuvien p-arvojen määrittämiseen Kenward-Rogersin approksimaatiota. Samainen toiminnallisuus löytyy nykyisin myös R:stä laajennuspaketin doBy funktiosta KRmodcomp(). Kannattaa tutustua, jos p-arvot ovat oleellisia tuloksia raportoitaessa!
Louhos-blogi
Ohessa linkki uuteen mielenkiintoiseen suomenkieliseen blogiin, jossa käsitellään asioita mm. R-kieltä käyttäen:
Kannattaa pitää blogia silmällä!
Karttapallo
Paul Murrell on julkaissut alunperin vuonna 2005 painetusta kirjastaan R Graphics uuden päivitetyn ja laajennetun version.
Aiempikin versio oli kultakaivos, mutta tästä opin jälleen uusia temppuja, nimittäin interaktiivista grafiikkaa käyttäen toteutetun karttapallon. Siinä pallon pinnan tekstuurina käytetään erillistä png-muotoista kuvatiedostoa. Toteutus menee pääpiirteissään seuraavasti.
library(rworldmap) library(rgl) # Luodaan tekstuurina käytettävä maailmankartta png("worldsmall.png", width=2000, height=1000) par(mar=c(0,0,0,0)) plot(getMap(resolution="medium"), lwd=1, col="#009900", bg="#00CCFF44") dev.off() # Muodostetaan pallo, ja päällystetään se tekstuurilla lat <- matrix(seq(90,-90, len=200)*pi/180, 200, 200, byrow=TRUE) long <- matrix(seq(-180, 180, len=200)*pi/180, 200, 200) r <- 6378.1 # radius of Earth in km x <- r*cos(lat)*cos(long) y <- r*cos(lat)*sin(long) z <- r*sin(lat) # Plotataan pallo rgl-paketin funktioita käyttäen open3d(windowRect=c(100, 100, 300, 300)) clear3d("all") light3d() persp3d(x, y, z, col="white", texture="worldsmall.png", specular="black", axes=FALSE, box=FALSE, xlab="", ylab="", zlab="", normal_x=x, normal_y=y, normal_z=z) par3d(userMatrix=rotationMatrix(-pi/2, 1, 0, 0)%*% rotationMatrix(-30/180*pi, 0, 0, 1)%*% rotationMatrix(45/180*pi, 1, 0, 0), zoom=2/3)
Tuloksena on siis interaktiivinen grafiikka, jota voi pyöritellä ja lähentää ja loitontaa. Tulos näyttää ruutukaappauksena seuraavalta:
