Manual de R (Estadística)/Introducció. Primers passos

El R és un llenguatge orientat al processament i anàlisi estadística de dades. És un derivat d'un llenguatge anterior anomenat S. Des de fa uns vint anys el desenvolupa un grup de voluntaris de diferents països. És un programari lliure i gratuït, publicat sota llicència GNU GPL.

El R pot arribar a semblar un infern. De fet hi ha un llibre sobre R amb aquest títol (The R Inferno) , escrit, com no podia ser d'altre manera, pel Sr. Burns (Patrick Burns) amb l'ajuda del famós escriptor italià el Sr. D. Alighieri . El llibre comença: "Si utilitzeu el R i creieu que esteu a l'infern, [...]".

Aquest llibre pretén que, mitjançant unes explicacions (presumptament) clares i amb molts exemples (en un llibre electrònic les limitacions d'espai no són tant importants), el R no arribi a ser un infern.

Particularitats del R

R en certa manera és molt diferent d'altres programaris estadístics com el SPSS o el SAS. Algunes particularitats que convé tenir en ment al hora d'utilitzar-lo (pàg. 6):

• Els altres programaris estadístic solen proporcionar un gran nombre de resultats i hom ha de cercar entre ells els que li interessen. R per omissió mostra molts pocs resultats i és l'usuari qui te que demanar més detalls. Com es veurà més endavant, algunes de les ordres per obtenir més detalls són (pàg. 46):
  • print: mostra un resum simple dels resultats
  • summary: mostra resum més detallat
  • df.residual: mostra el nombre de graus de llibertat residuals
  • coef: mostre els coeficients estimats
  • residuals: mostra els residuals
  • deviance: mostra la llunyania
  • fitted: mostra els valors ajustats pel model
  • logLik: calcula el logaritme de la versemblança i el nombre de paràmetres del model
  • AIC: calcula el criteri d’informació d’Akaike o AIC (depèn de logLik())
• El R és el fruit de la col·laboració descentralitzada de moltes persones. En principi, tothom pot fer un paquet per ell i el nombre d'autors que ha fet almenys un paquet és considerable. Les avantatges són immenses ja que les contribucions intel·lectuals i de coneixement segurament no es troben en cap altre programa estadístic. Te el problema de que, malgrat existir un nucli central de desenvolupadors, els paquets poden ser molt heterogenis quant al seu comportament.
• El resultat de la gran majoria d'ordres es poden utilitzar com l'entrada d'altres ordres.
• Els resultats no es guarden en cap finestres de resultats com ho fa, per exemple, el SAS. Si no es dóna l'ordre de salvar-los en algun calaix (que els informàtics anomenen objecte), no es pot accedir a ells més endavant. Això es fa assignant el resultat a un calaix amb un nom determinat. Per exemple, si es volgués salvar els resultats d'una prova t de la comparació de l'edat entre homes i dones (d'unes dades que estan al fitxer DADES) a un calaix amb el nom DADES.PROVA_T s'executaria l'ordre:

DADES.PROVA_T = t.test(DADES2$EDAT ~ DADES2$SEXE)

En qualsevol moment es pot accedir als resultats:

print(DADES.PROVA_T)

• En comparació amb SAS (que té el senzill mètode de ods pdf) o el SPSS, la capacitat del R per traslladar els resultats en un informe de qualitat és molt limitada o molt complicada, encara que aquest aspecte s'està millorant amb el paquet R2HTML i amb altres desenvolupaments bastants complicats (p. ex., Sweave o knitr).

Primers passos

A continuació es comenta una petita mostra de com s'utilitza el R. Es mostra com realitzar tasques molt senzilles. En primer lloc cal tenir-lo instal·lat i executar-lo. Si s'utilitza Windows o Mac OS, la solució més fàcil és utilitzar la interfície gràfica d'usuari R (fer clic sobre la seva icona). Si es fa servir Linux, s'obre un terminal i es tecleja R en l'indicatiu.

En general, quan s'obre R en una terminal, veureu un missatge semblant al següent:

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>

El R funciona a base d'instruccions. Les instruccions s'ha d'escriure després de l'indicador ">". En els exemples presentats aquí, el codi que s'ha d'escriure i executar és el que hi ha després del símbol ">".

Si hom només necessita l'estadística convencional (p. ex., descriptiva, khi quadrat, t-Student, Wilcoxon o regressió lineal), i no ha de fer manipulacions complicades de les dades, probablement amb el R Commander o el RKWard en tindrà més que suficient. Funciona amb menús i no representa cap complicació.

Tant en Windows com en Linux, pot utilitzar el R Commander que en realitat és un conjunt de diversos paquets i que facilita utilitzar aquelles instruccions més complicades. Primer s'ha d'instal·lar el paquet Rcmdr (veure Instal·lar un paquet) i un cop instal·lat es pot executar amb la ordre library(Rcmdr):

> library(Rcmdr)

Al executar-se aquesta ordre, s'obre la finestra del R Commander i ja es pot començar a treballar mitjançant menús.

Si treballa amb linux, també pot utilitzar el RKWard, una interfície gràfica del R sota KDE que funciona amb menús que potser és més amigable que el R Commander. Te grans avantatges:

• Emmagatzema els resultats en format html (al fitxer ~/.rkward/rk_out.html) i per tant, són fàcils de compartir ja que te l'opció de copiar els resultats en un altre fitxer i en el directori que se li indiqui. Aquest fitxer es pot compartir amb altres persones.
• La gestió de l'espai de treball és molt senzilla.
• Integra un editor molt adequat (Kate).

En Debian (i suposo que en Ubuntu), rkward és un paquet que s'instal·la com qualsevol altre paquet de la distribució.

Si no voleu treballar amb menús, sinó que voleu escriure els programes com es sol fer, p. ex., en SAS, ho podeu fer amb el R normal escrivint el codi amb un editor i fent copia-enganxar o bé amb l'editor del RKWard. Aquesta darrera manera de treballar és molt còmode. Per poder fer-ho, cal conèixer el llenguatge del R en el que, a grans trets, es poden diferenciar dos grups:

• Els elements del llenguatge de programació amb, p. ex., if, for. Si al indicador del R es tecleja if i després Enter, apareix un signe "+" indicant que està esperant més ordres.
• Les ordres (funcions) que corresponent a les funcionalitats del R. Una funció és, p. ex., la plot que realitza, p. ex., un diagrama de dispersió. Si es tecleja al indicador plot, es mostre el codi de la funció.

Com tot programa té la seva sintaxi i la seva manera d'emmagatzemar les dades i és el que explicarem.

En el indicador de la terminal del R (el símbol ">") s'escriuen expressions, que el R les interpreta, avalua i, finalment, les executa. Les expressions més utilitzades són:

• Els noms que corresponent al noms dels calaixos (o objectes) on s'han anat guardant coses. En el R quan es prem la tecla Enter després d'haver escrit un nom, mostrar el contingut del calaix.
• Les crides a ordres (que en R s'anomenen funcions) que realitzen els càlculs especificats.

Una taula de dades

Suposem que es guarda en el calaix DADES, el sexe d'un conjunt d'individus. Si després es tecleja el nom del calaix, el R mostre el seu contingut:

# Es crea una taula de dades:
> DADES <- data.frame(SEXE=c("Dona", "Home", "Home", "Dona"))

> DADES
  SEXE
1 Dona
2 Home
3 Home
4 Dona

Una ordre es crida habitualment teclejant el seu nom seguit d'una llista d'arguments entre parèntesis. Per exemple, l'ordre table mostra la freqüència d'homes i dones:

> table(DADES$SEXE)

Dona Home 
   2    2

Del codi anterior es dedueix que en R per fer referència a la variable SEXE de la taula de dades DADES, s'escriu DADES$SEXE.

Hem de saber que R si fa diferència entre majúscules i minúscules i per tant per a ell les següents variables són diferents: "AA", "Aa", "aA" i "aa".

Quan s'assigna un valor a una variable, el valor que aquesta pren no surt per pantalla i per tant s'ha de cridar una altra vegada la variable per visualitzar el valor que ha pres.

> EDAT = 45                       # Al assignar el valor, aquest no es visualitzar

> EDAT                            # Cal cridar la variable per visualitzar el valor de l'edat
[1] 45

Si es deixa de banda (de moment) el llenguatge de programació, treballar amb el R consistirà simplement en identificar l'ordre (funció) que realitza la tasca que es vol fer (p. ex., una prova de khi quadrat) i determinar com s'ha d'executar a fi de que la faci. Es veurà amb un exemple.

En la taula de dades Aids2 del paquet MASS hi ha el sexe i l'estat (viu o mort) d'un conjunt de malalts amb sida. Si es vol visualitzar les dades, primer s'han de carregar en memòria (ordre data()) i després es pot visualitzar l'estructura de la taula de dades (ordre str()) i els sis primer registres (ordre head()):

> data(Aids2, package="MASS")                 # És carrega el paquet on esta les dades de mostra
> str(Aids2)                                  # Mostra l'estructura del fitxer 
'data.frame':	2843 obs. of  7 variables:
 $ state  : Factor w/ 4 levels "NSW","Other",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ sex    : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
 $ diag   : int  10905 11029 9551 9577 10015 9971 10746 10042 10464 10439 ...
 $ death  : int  11081 11096 9983 9654 10290 10344 11135 11069 10956 10873 ...
 $ status : Factor w/ 2 levels "A","D": 2 2 2 2 2 2 2 2 2 2 ...
 $ T.categ: Factor w/ 8 levels "hs","hsid","id",..: 1 1 1 5 1 1 8 1 1 2 ...
 $ age    : int  35 53 42 44 39 36 36 31 26 27 ...
> head(Aids2)                                 # Llista els 6 primer registres del fitxer
  state sex  diag death status T.categ age
1   NSW   M 10905 11081      D      hs  35
2   NSW   M 11029 11096      D      hs  53
3   NSW   M  9551  9983      D      hs  42
4   NSW   M  9577  9654      D    haem  44
5   NSW   M 10015 10290      D      hs  39
6   NSW   M  9971 10344      D      hs  36

Una prova de khi quadrat

Si es vol fer una prova de khi quadrat per estudiar la relació entre el sexe i l'estat, primer s'ha de saber amb quina ordre el R calcula aquesta prova. Cal anar a l'ajuda del R. Amb l'ordre:

help.start()

s'obre una pàgina web amb enllaços als documents d'ajuda. També es pot demanar ajuda amb l'ordre:

?<ordre>

Per exemple, si es vol que es mostre l'ajuda sobre l'ANOVA:

?anova

Pel cas que ens ocupa. si es vol buscar quelcom similar a chis en el R existeixen dues formes:

1. apropos("chis"): cerca entre els noms de les ordres (anomenades en R funcions), la paraula chis.
2. help.search("chis"): cerca les ordres que en que a la seva pàgina d'ajuda apareix el text chis.

Hi han cinc funcions (o ordres) amb el text chis en el seu nom:

> apropos("chis")
[1] "chisq.test" "dchisq"     "pchisq"     "qchisq"     "rchisq"

Hi han 12 ordres en les que el text chis apareix en la seva pàgina d'ajuda:

> help.search("chis")
Help pages:
coin::chisq_test  Independence in Three-Way Contingency Tables
Hmisc::biVar  Bivariate Summaries Computed Separately by a Series of Predictors
multcomp::summary.glht  Methods for General Linear Hypotheses
sfsmisc::p.dnorm  Plot Parametric Density Functions
survey::pchisqsum  Distribution of quadratic forms
survey::svreptable  Contingency tables for survey data
TeachingDemos::chisq.detail  Print details of a chi-squared test
timeDate::DaylightSavingTime  Daylight Saving Time Rules
VGAM::chisq  Chi-squared Distribution
stats::chisq.test  Pearson's Chi-squared Test for Count Data
stats::Chisquare  The (non-central) Chi-Squared Distribution
stats::stats-defunct  Defunct Functions in Package 'stats'

Clarament sembla que l'ordre que s'està cercant és chisq.test (Pearson's Chi-squared Test for Count Data).

Ara cal saber com s'ha d'executar aquesta ordre. Es pot cerca ajuda escrivint help(chisq.test) al indicador del R i prement la tecla Enter. Informa de:

Description: chisq.test performs chi-squared contingency table tests and goodness-of-fit tests.

Usage:

chisq.test(x, y = NULL, correct = TRUE,

p = rep(1/length(x), length(x)), rescale.p = FALSE,

simulate.p.value = FALSE, B = 2000)

No és molt clar el que vol dir. Al final de l'ajuda apareixent varis exemples del tipus:

## From Agresti(2007) p.39
M <- as.table(rbind(c(762, 327, 468), c(484,239,477)))
dimnames(M) <- list(gender=c("M","F"),
                    party=c("Democrat","Independent", "Republican"))
(Xsq <- chisq.test(M))  # Prints test summary
Xsq$observed   # observed counts (same as M) 
Xsq$expected   # expected counts under the null
Xsq$residuals  # Pearson residuals
Xsq$stdres     # standardized residuals

Continua, almenys per a mi, sense quedar gaire clar com s'ha de fer un simple khi quadrat amb R. Per aquest motiu s'escriu aquest llibre, per proporcionar explicacions quelcom més (esperem) entenedores. Aquestes explicacions per la prova de khi quadrat es poden trobar al apartat corresponent.

Però tornant a la prova de khi quadrat de les dades Aids2, hom es pot tornar a fixar en les primeres línies de l'ajuda, allà on posa:

chisq.test(x, y = NULL, correct = TRUE,

p = rep(1/length(x), length(x)), rescale.p = FALSE,

simulate.p.value = FALSE, B = 2000)

Es pot interpretar que x i y son dos variables i, deixant de banda la resta de codi, executar:

> chisq.test(Aids2$sex, Aids2$state)

	Pearson's Chi-squared test

data:  Aids2$sex and Aids2$state 
X-squared = 5.8235, df = 3, p-value = 0.1205

Fent-ho així, s'obté el que es volia: la prova de khi quadrat. Com ja s'ha comentat, una manera de fer referència a una variable que està en una taula de dades és <Nom taula de dades>$<nom de variable>.

Els valors desconeguts i altres caràcters especials

Sovint per alguna de les variables no es coneix el valor d'algun individu. Per exemple, és possible que en un estudi la pressió arterial d'algun malalt no s'hagi registrat. En R els valors desconeguts es representen per NA ("not available"). A diferència de, p. ex., el SAS, en R s'utilitza el mateix símbol per les variables quantitatives i qualitatives.

S'ha de tenir cura amb els valors desconeguts, ja que poden provocar que els resultats siguin estranys. P. ex., si es té una taula de dades amb el sexe i edat de quatre persones, però l'edat de la segona és desconeguda (valor NA):

> DADES <- data.frame(SEXE=c("Dona", "Home", "Home", "Dona"), EDAT=c(20,NA,30,35))  # Dades de 4 persones
> head(DADES)                                                                       # Visualitzar les dades
  SEXE EDAT
1 Dona   20
2 Home   NA
3 Home   30
4 Dona   35

Si es vol conèixer l'edat mitjana dels quatre individus, es pot utilitzar l'ordre mean(), però el resultat no és l'esperat:

> mean(DADES$EDAT)
[1] NA

El R no calcula la mitjana, ja que un dels valors és desconegut. Per evitar aquest tipus de problemes, la majoria de funcions numèriques tenen l'opció na.rm (NA remove) que si se li assigna el valor de TRUE elimina els valors desconeguts en l'anàlisi i només fa els càlculs amb els valors restants (, pàg. 80):

> mean(DADES$EDAT, na.rm=TRUE)
[1] 28.33333

Cura, doncs amb els NA, ja que a diferència de, p. ex., el SAS, hom pot tenir la necessitat d'excloure's de forma explícita al cridar una ordre.

Altres caràcters especials són:

• NaN (not a number): els valors impossibles com, p.ex., zero dividit per zero:

 > 0/0
 [1] NaN

• Inf: infinit. Per exemple, 1 dividit per zero dona $+\infty$:

 > 1/0
 [1] Inf

• -Inf: menys infinit.

Recodificar un valor a desconegut

Si en les dades, el valor desconegut de, p. ex., l'edat està codificat amb un "-9", per indicar-li al R que el "-9" és el valor pel desconegut i no una edat correcta, es pot fer:

> DADES <- data.frame(SEXE=c("Dona", "Home", "Home", "Dona"), EDAT=c(20,-9,30,35))  # Dades
> head(DADES)
  SEXE EDAT
1 Dona   20
2 Home   -9
3 Home   30
4 Dona   35

> DADES[DADES$EDAT==-9, "EDAT"] <- NA                  # Es recodifica el "-9" com a desconegut per la variable EDAT
> head(DADES)
  SEXE EDAT
1 Dona   20
2 Home   NA
3 Home   30
4 Dona   35

Mes detall sobre els valors desconeguts es comentaran en un apartat específic posterior.

Sortir de R

Ara ja pot sortir del R teclejant:

q()

Si s'especifica l'argument no, no es guarda la sessió de R:

q("no")

Tipus de variables estadístiques

Les variables, fent una abusiva, però de moment útil simplificació, poden ser de dos tipus segons el tipus de característica mesurada:

• Qualitatives: són les que fan referència a una qualitat dels objectes estudiats (p. ex., sexe, grup sanguini o estat civil) i prenen valors que no són numèrics.
• Quantitatives: són les que prenen valors numèrics (p. ex, edat, pes, alçada o nombre de leucòcits).

Segons el seu paper en l'anàlisi, una variable pot ser:

• Depenent: La variable que es pretenen estimar. És la d'interès de l'estudi. Per exemple, si es vol estimar quin és l'efecte del nivell socioeconòmic sobre la mortalitat, la variable dependent és la defunció.
• Independent: La que defineix els diferents grups en els que es vol mesurar la variable depenent. Per l'exemple anterior, la variable independent és el nivell socioeconòmic, ja que és vol determinar la mortalitat en els diferents grups definits pel nivell socioeconòmic.

Com utilitzar el llibre

El llibre està ordenat per tasques i les proves estadístiques més bàsiques, pel nombre de variables a tractar i el tipus de les mateixes.

Una primera part està dedicada a les dades. En primer lloc es presenta com llegir les dades externes i com tornar-les a desar en un fitxers extern. En segon lloc, com es poden modificar les dades que s'han d'analitzar.

Una segona part presenta l'anàlisi de les dades. Es presenten els gràfics, després els estadístics univariats i finalment las proves d'hipòtesi. Aquestes darreres estan ordenades segons el tipus de variable a analitzar i del seu nombre (primer el cas de dos variables i després els models amb més de dos variables).

Si el problema només te dos variables, existeixen fonamentalment tres possibilitats ( pàg. 270-1):

• Les dos són qualitatives. Anar a l'apartat de Relació entre dos variables qualitatives. Explica com realitzar la prova de la khi quadrat.
• Una és qualitativa però l'altre quantitativa. Anar a l'apartat de Relació entre una variable qualitativa i una quantitativa.
• Les dues són quantitatives. Anar a l'apartat de Relació entre dos variables quantitatives.

Si el problema té una variable dependent i diverses independents, s'anirà a l'apartat Relació entre una variable dependent i varies independents. Però si la variable dependent és del tipus temps a un esdeveniment, anar a l'apartat Anàlisi de la supervivència. Si les dades són més complexes, segurament caldrà anar a mètodes d'anàlisi més esotèrics, com els models lineals generalitzats, els models lineals mixtes o els models mixtes lineals generalitzats, segons la nomenclatura utilitzada en SAS (XXX desconec si és aplicable al R XXX).

En una tercera part, es presenta com fer informes de qualitat, els elements dels llenguatge del R i l'ús d'aquest per escriure funcions.