R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

SAS ja R: esimerkki SAS-koodin toteuttamisesta R:llä

Pidin maaliskuun lopussa turkulaisille tilastotieteilijöille ja biostatistikoille esityksen R:stä ja sen mahdollisuuksista datan käsittelyssä. Turussa on vahvaa SAS-osaamista etenkin lääketieteellisen tutkimuksen puolella; R-osaamistakin toki löytyy, mutta paino on vahvasti SAS:n puolella.

Valmistelimme esityksen yhdessä turkulaisen biostatistikon Lauri Sillanmäen kanssa. Johtoajatuksena oli se, että SAS:n hyvin tuntevat tilastotieteilijät näkisivät miten SAS:lla tehdyn asian voi tehdä R:llä. Usein on myös niin, että oppikirjoissa ja oppaissa data on usein valmiissa muodossa analyointia varten, mutta oikeassa elämässä dataa täytyy monesti pyöritellä paljon ennen kuin data on analysoitavassa muodossa.

Esitystä varten Lauri valmisteli SAS:lla esimerkkiaineiston ja pyörittelyt ja minä toteutin vastaavat pyörittelyt käyttäen R:ää. Tavoitteena oli pysytellä lähinnä datan käsittelyssä eikä niinkään vertailla tilastollisia menetelmiä tai mahdollisuuksia. Tarkoitus ei myöskään ollut tarkastella välineiden paremmuutta vaan pyrkiä saamaan käsitys siitä, miten asia voidaan toteuttaa eri välineillä. Esimerkkidataan generoitiin myös puuttuvia arvoja, jotta näkisimme miten näiden käsittely sujuu eri välineillä. Lopputulokseen tietysti vaikutti se, että miten olemme tottuneet käyttämään ohjelmia: joku toinen olisi saanut SAS-/R-skriptin näyttämään aivan erilaiselta päätyen kuitenkin suurinpiirtein samoihin lopputuloksiin.

Sitten asiaan:

Muodostetaan aluksi SAS:lla aineisto, jota molemmat ohjelmat käsittelevät:

/*ESIMERKKIDATOJEN GENEROINTI*/
data perus1 (drop= i j);
     do id=1 to 50;
        *Height;
        array height(*) height1 height2 height3;
        do i=1 to 3;
           height(i)=ranuni(34523)*8.5+(5*i)+150;
           if ranuni(123)/i<0.1 then height(i)=.;
        end;
        *Weight;
        array weight(*) weight1 weight2 weight3;
        do j=1 to 3;
           weight(j)=ranuni(34523)*15+(10*j)+20;
           if ranuni(673)/j<0.1 then weight(j)=.;
        end;
        *Gender;
        sp=ranbin(123,1,0.4);
        if ranuni(456)<0.1 then sp=.;
 
        rannum=ranuni(999);
        output;
     end;
run;
 
proc sort data=perus1 out=perus1(drop=rannum); by rannum; run;
 
data diagnoosi0;
     do id=1 to 50;
        random=scan('1 2 3 4',rand('TABLE',0.55,0.3,0.2,0.05));
        output;
     end;
run;
data diagnoosi1 (drop=random nro);
     set diagnoosi0;
     do nro=1 to random;
           diag=catt(byte(floor(ranuni(2349)*3)+65),put(ranuni(489)*9,z2.));
           pp=ceil(ranuni(11)*30);
           kk=ceil(ranuni(22)*12);
           yy=ceil(ranuni(33)*20)-10;
           if yy<0 then yy=100+yy;
        output;
     end;
run;
 
/*Tallennetaan datat R-käyttöä varten*/
libname temp "c:\temp";
 
 
/*Pysyvät SAS-datat*/
data temp.perus1;
     set perus1;
run;
 
data temp.diagnoosi1;
     set diagnoosi1;
run;
 
/*CSV-datat*/
proc export data=perus1
            outfile="c:\temp\perus1.csv"
            dbms=DLM replace;
     delimiter=','; 
     putnames=yes;
run;
 
proc export data=diagnoosi1
            outfile="c:\temp\diagnoosi1.csv"
            dbms=DLM replace;
     delimiter=','; 
     putnames=yes;
run;

Tämän jälkeen käsitellään dataa:

/*DATOJEN VALMISTELU*/
proc print data=perus1 (obs=10); 
run;
 
proc print data=diagnoosi1 (obs=10); 
run;
 
proc format; 
     value spf 
        0="Male"
        1="Female";
     value bmif 
        low-18.5="Underweight"
       18.5-25.0="Normal"
       25.0-30.0="Overweight"
       30.0-high="Obese";
run;
 
/*Lasketaan painoindeksit, lisätään sukupuoli- ja BMI-formaatit ja muuttujalabelit*/
data perus2;
     set perus1;
     bmi1=weight1/(height1/100)**2;
     bmi2=weight2/(height2/100)**2;
     bmi3=weight3/(height3/100)**2;
     format sp spf. bmi1--bmi3 bmif.;
     label bmi1="Body Mass Index at time 1"
           bmi2="Body Mass Index at time 2"
           bmi3="Body Mass Index at time 3";
run;
 
/*Muodostetaan päivämäärämuuttuja päivä-, kuukausi- ja vuositiedon pohjalta ja
tehdään diagnoositiedosta erillinen muuttuja*/
data diagnoosi2 (drop=pp kk yy);
     set diagnoosi1;
   	 if diag ne "" then pvm=mdy(kk,pp,yy);
     format pvm EURDFDD10.;
     if diag ne "" then diagc=substr(diag,1,1);
run;
 
proc sort data=diagnoosi2; by id; run;
 
proc sort data=perus2; by id; run;
 
/*Yhdistetään sukupuolitieto diagnoosidataan*/
data diagnoosi3;
     merge diagnoosi2 perus2(keep=id sp);
     by id;
run;
 
proc print noobs; run;
 
/*DATOJEN TARKASTELUA*/
/*Kohta 1: Lasketaan jokaisella henkilölle A, B ja C-pääryhmän diagnoosien lkm*/
proc sort data=diagnoosi3; by diagc; run;
 
proc freq data=diagnoosi3 noprint;
     tables id*diagc /out=diagnoosi3f(drop=percent) sparse;
run;
 
proc sort data=diagnoosi3f; by diagc; run;
 
/*Tulostetaan frekvenssit*/
proc freq data=diagnoosi3f; 
     tables count;
     by diagc;
run;
 
/*Tehdään frekvensseistä pylväskuvio*/
proc gchart data=diagnoosi3f;
     hbar count /type=freq subgroup=diagc discrete ;
run;

/*Lisätään frekvenssitieto perusdataan*/
proc sort data=diagnoosi3f; by id; run;
 
proc transpose data=diagnoosi3f
               out=diagnoosi4f(drop=_name_ _label_)
               prefix=count_; 
     var count; 
     by id; 
     id diagc;
run;
 
data perus3;
     merge perus2 diagnoosi4f;
     by id;
run;
 
/*Kohta 2: Tehdään sukupuolittaiset pituuksien keskiarvokuvat ajan suhteen*/
proc sort data=perus2; by id sp; run;
 
proc transpose data=perus2
               out=perus2_trans(rename=(col1=height));
     by id sp;
     var height1 height2 height3;
run;
 
/*Valmistellaan dataa, jotta miesten ja naisten käyrät saadaan limitettyä horisontaalisesti*/
data perus2_trans;
     set perus2_trans;
     time=abs(substr(_name_,7,1));
     if sp=0 then time2=time-0.05;
     if sp=1 then time2=time+0.05;
run;
proc print; run;
 
/*Luodaan kuvaa varten erillinen formaatti*/
proc format;
     value timef (fuzz=0.1)
        1="1"
        2="2"
        3="3";
run;
 
/*Tuotetaan keskiarvoprofiilikuva, joissa janat kuvaavat keskihajontaa. X-akselilla
on vielä tarpeettomat arvot "1" ja "4", mutta niitä ei saa helposti pois*/
 
symbol1 v=none interpol=std1jt ci="blue" line=1  w=2 v=none cv="blue" h=0.1 mode=include;
symbol2 v=none interpol=std1jt ci="red"  line=33 w=2 v=none cv="red"  h=0.1 mode=include;
axis1 minor=none;
 
proc gplot data=perus2_trans;
     where sp in(0,1);
     plot height*time2=sp /haxis=axis1; 
     format time2 timef.;
run; quit;

/*Kohta 3: Lasketaan perustunnareita aikapisteittäin ja sukupuolittain*/
proc means data=perus2_trans n min max mean median stddev;
     var height; 
     class time; 
run;
 
proc means data=perus2_trans n min max mean median stddev;
     where sp in(0,1);
     var height; 
     class sp time; 
run;
 
/*TALLENNETAAN TÄRKEIMMÄT KÄYTTÖDATAT MYÖHEMPÄÄ KÄYTTÖÄ VARTEN*/
libname data "c:\data";
 
data data.diagnoosi3;
     set diagnoosi3;
run;
 
/*Ylimääräinen diagnoosidata, jossa vain miesten lääkkeet 21.vuosituhannelta*/
data diagnoosi3_A_males;
     set diagnoosi3;
     where sp=0 and pvm>='01Jan2000'd;
run;
 
data data.perus2;
     set perus2;
run;
 
data data.perus2_trans;
     set perus2_trans;
run;

Koitetaan tehdä vastaavat asiat R:llä. Olen pyrkinyt osoittamaan vastaavat aloituskohdat samoilla teksteillä kun SAS-ajojonossakin. Tällöin kohta, jossa sama asia tehdään, löytyy toivottavasti helpommin.

R:llä homma voisi siis mennä vaikka näin:

#----------------------------- ESIMERKKIDATOJEN TUONTI ----------------------------------
#Asetetaan työskentelyalue
setwd("C:/Temp/")
 
#Ladataan paketit
library(ggplot2)
library(gregmisc)
library(doBy)
 
#Muutetaan .csv-muotoon sassissa ja tuodaan sisään
diagnoosi1 <- read.csv("diagnoosi1.csv", na.strings=".")
perus1 <- read.csv("perus1.csv", na.strings=".")
 
#----------------------------- DATOJEN VALMISTELU ------------------------------
#Lasketaan painoindeksit, lisätään sukupuoli- ja BMI-formaatit ja muuttujalabelit
perus2 <- transform(perus1, bmi1 = weight1/(height1/100)**2,
                            bmi2 = weight2/(height2/100)**2,
                            bmi3 = weight3/(height3/100)**2,
                            spf = ifelse(sp == 0, "Male", "Female"))
head(perus2)
 
temp <- as.data.frame(lapply(perus2[, c("bmi1", "bmi2", "bmi3")], cut, breaks = c(-Inf, 18.5, 25, 30, 50), labels =
c('Underweight', 'Normal', 'Overweight', 'Obese')))
perus2 <- data.frame(cbind(perus2, temp))
head(perus2)
 
#Muodostetaan päivämäärämuuttuja päivä-, kuukausi- ja vuositiedon pohjalta ja tehdään diagnoositiedosta erillinen muuttuja
g <- with(diagnoosi1, paste(pp, ".", kk, ".", yy, sep=""))
 
diagnoosi2 <- transform(diagnoosi1, pvm = as.Date(strptime(g, "%d.%m.%y"), "%Y-%m-%d"),
                            diagc = substr(diag, 1, 1))
 
#Yhdistetään sukupuolitieto diagnoosidataan
diagnoosi2$spf <- perus2$spf[match(diagnoosi2$id, perus2$id)]
diagnoosi3 <- diagnoosi2[order(diagnoosi2$id),]
head(diagnoosi3)
 
#----------------------------- DATOJEN TARKASTELU ------------------------------
#Kohta 1: Lasketaan jokaisella henkilölle A, B ja C-pääryhmän diagnoosien lkm
#Jätetään laskematta prosentit, kumulatiivinen frekvenssi ja kumulatiivinen prosentti
g <- with(diagnoosi3, table(id, diagc))
g1 <- do.call(smartbind, apply(g, 2, table))
g1
 
#Tehdään frekvensseistä pylväskuvio
g1[is.na(g1)] <- 0
g1$diagc <- rownames(g1)
m <- melt(g1, id.vars = "diagc")
m$variable <- factor(m$variable, rev(unique(m$variable)))
ggplot(m, aes(x = variable, y = value, fill = diagc)) + geom_bar() + coord_flip() + scale_x_discrete("Count") + scale_y_continuous("Frequency")

#Lisätään frekvenssitieto perusdataan
g <- reshape(as.data.frame.table(g), v.names="Freq", idvar="id", timevar="diagc", direction="wide")
perus3 <- merge(perus2, g)
head(perus3)
 
#Kohta 2: Tehdään sukupuolittaiset pituuksien keskiarvokuvat ajan suhteen
perus2_trans <- reshape(perus2, varying = list(c("height1", "height2", "height3"),
                               c("weight1", "weight2", "weight3"),
                               c("bmi1", "bmi2", "bmi3")),
                v.names=c("heigth","weigth", "bmi"),
                timevar = "time",
                times = factor(c(1, 2, 3)),
                idvar = c("id", "sp"),
                direction = "long",
                new.row.names = 1:(nrow(perus2)*3))
perus2_trans <- perus2_trans[order(perus2_trans$id, perus2_trans$time),]
 
#Tehdään kuva (käytetään tässä luottamusvälejä)
erb <- summaryBy(heigth ~ spf + time, data = subset(perus2_trans, !is.na(spf)), FUN = ci, na.rm=T)
erb$spf <- factor(erb$spf, rev(unique(erb$spf)))
limits <- aes(ymax = erb[, "heigth.CI upper"], ymin = erb[, "heigth.CI lower"])
p <- ggplot(erb, aes(colour=spf, y=heigth.Estimate, x=time))
p + geom_line(aes(group=spf), position = position_dodge(0.4)) + geom_errorbar(limits, width=0.2, position = position_dodge(0.4)) + scale_colour_manual(values = c("blue","red")) + scale_y_continuous(limits=c(155, 174), breaks=155:174)

#Kohta 3: Lasketaan perustunnareita aikapisteittäin ja sukupuolittain.;
statfun <- function(x) c(n = nobs(x),
                         min = min(x, na.rm=T),
                         max = max(x, na.rm=T),
                         mean = mean(x, na.rm=T),
                         median = median(x, na.rm=T),
                         stdev = sd(x, na.rm=T))
 
summaryBy(heigth ~ time, data = subset(perus2_trans, !is.na(spf)), FUN = statfun)
summaryBy(heigth ~ spf + time, data = subset(perus2_trans, !is.na(spf)), FUN = statfun)
 
#---------TALLENNETAAN TÄRKEIMMÄT KÄYTTÖDATAT MYÖHEMPÄÄ KÄYTTÖÄ VARTEN----------
kirjoita <- function(x, nimi, ...) {
          write.csv2(x, file=paste(getwd(), "/", Sys.Date(), "_", nimi, ".csv", sep=""), ...)
          }
 
kirjoita(diagnoosi3, "diagnoosi3", row.names=FALSE)
 
#Ylimääräinen diagnoosidata, jossa vain miesten lääkkeet 21.vuosituhannelta
diagnoosi3_A_males <- subset(diagnoosi3, (spf == "Male" & pvm >= "2000-01-01"))
kirjoita(diagnoosi3_A_males, "diagnoosi3_A_males", row.names=FALSE)
 
kirjoita(perus2, "perus2", row.names=FALSE)
 
kirjoita(perus2_trans, "perus2_trans", row.names=FALSE)

Ohessa SAS:n ja R:n output:
SAS-output
R-output

Ohessa lähtödatat, joita näissä ajoissa käytettiin:
perus1
diagnoosi1