1 R amb RStudio

Autors: Muñoz R. & Gascón S.

En aquesta secció d’introducció a R amb RStudio veureu com carregar dades, comprovar que s’ha fet correctament i exploració incial de les dades. A més veurem com adequar i comprovar que les dades tenen el format necessari per a realitzar una anàlisi d’ANOVA¹.

A l’exemple que anem a treballar farem una ANOVA per saber si hi ha diferències d’alçada (height; numèric continu) en funció del sexe (sex; categòric i binari).

Figure 1.1: A) histograma de les alçades per gènere (binari). B) boxplots de les alçades per gènere (binari).

Abans de començar a treballar cal fer algunes passos:

Establir el que s’anomena directori de treball o alternativament també podeu utlitzar Projecte d’RStudio: L’espai on l’R efectuarà tota la feina. És on heu de guardar les dades, perquè R les trobi i les pugui carregar
Carregar les dades a l’R

1.1 Establir directori de treball

Per tal que R trobe l’arxiu a cada sessió cal definir quina és la nostra carpeta/directori de treball o bé generar un projecte d’Rstudio.

Podeu saber, en quin directori de treball està funcionant l’R si piqueu getwd(). Llavors, R us retorna la ruta del directori de treball actual.

La carpeta/directori de treball es defineix mitjançant la funció set working directory setwd() o bé pel menú superior d’RStudio (allà on posa, File,Edit,Code) on posa Session> Set Working directory>Choose Directory

Ull !!!! al establir la ruta (el camí de carpetes fins arribar a la carpeta que voleu utilitzar com a directori de treball), heu de vigilar de no tenir carpetes amb símbols estranys (per R els accents són símbols estranys) i **tampoc espais en els noms de les carpetes. Si us passa això, heu de canviar el nom d ela carpeta, ja que sinó teniu moltes probabilitats de que R no pugui trobar el camí* fins la vostra carpeta que serà el directori de treball.

Per solucionar problemes amb els path o rutes, l’RStudio té la opció de definir Projectes. Això permet que tot el que poseu dins la carpeta del projecte directament la podreu llegir sense haver d’anar canviant el directori de treball.

1.1.1 Per Windows

setwd("C:/Documents/AAD")

En el vostre cas haureu de substituir: “C:/Documents/AAD” per la vostra ruta. Alerta! Els símbols / han d’estar orientats cap a la dreta.

1.1.2 Per Mac

setwd("/AAD")

En el vostre cas haureu de substituir: “/AAD” per la vostra ruta. Alerta! Els símbols / han d’estar orientats cap a la dreta.

Una vegada indicat el directori de treball podem carregar les dades a R.

1.2 Crear i treballar amb Projecte d’RStudio

Quan es crea un Projecte tots els arxius queden vinculats al projecte, i per tant no cal establir cap directori de treball.

Com crear un nou Projecte:

Heu de seleccionar al menú File > New Project
S’obrirà una finestra demanant si es vol:

New Directory (creeu el projecte en una carpeta nova). Llavors seleccioneu el tipus de Projecte, en el nostre cas Empty Projecti assignem un nom al directori (carpeta) que es crearà i que serà el nom del Projecte. Per acabar, heu de fer clic a Create Project. Això, crea una carpeta i a dins un fitxer amb el mateix nom que la carpeta però amb extensió .Rproj
Existing Directory (creeu el projecte en una carpeta que ja existeix). Sel·leccioneu la carpeta on voleu posar el vostre projecte utilitzant el botó Browse. Un cop escollida feu clic a Create Project.
Version Control (no la escolliu)

Per obrir un projecte feu doble clic al arxiu amb extensió .Rproj, o el podeu obrir des del menú d’RStudio seleccionant la opció de File>Open Project

AVANTATGE de treballar amb Projectes, qualsevol arxiu o fitxer que creem queda vinculat automàticament a la carpeta del Projecte i quan guardem es guardarà allà.

1.3 Carregar dades per a l’anàlisi

1.3.1 format dades d’Excel .csv

Per obtenir dades d’Excel en format .csv, cal que guardeu el vostre arxiu d’excel com a text delimitat per comes, trobareu la opció si feu Guardar com

Un cop guardades en aquest format ja podeu utilitzar les següents instruccions per carregar les dades.

Les dades per aquest exercici les trobareu a l’Excel “dades_base.csv”. Cal comprovar/preparar les dades com en vist anteriorment.

Cada columna serà una variable.
Eviteu barrejar nombres i text dintre d’una columna.
Quan una dada és desconeguda, es deixa la cel·la en blanc o bé s’escriu NA (Not Available).
Eviteu apòstrofs (’), espais, #, etc. als encapçalaments.
Cada columna ha de tindre únicament una cel·la com a encapçalaments/nom.
No deixeu files en blanc entre fila de encapçalament i els valors.

pes_alumnes.data <- read.csv("dades_base.csv", header=T, sep=";")

Pel que fa als arguments:

pes_alumnes.data és el nom (objecte) que li donem a la taula de dades o data frame (Nom amb que és designen aquest tipus de taules a R).
<- és el símbol d’assignació. Assignem el què apareix a la dreta a un objecte (nom) que apareix a l’esquerre. Podeu comprovar que l’objecte ha estat assignat mirant si apareix a la finestra superior dreta d’RStudio anomenada Global Environment.
read.csv() és el nom de la funció que usem per a llegir l’arxiu de dades “dades_base.csv”.
header=T indica que la nostra taula de dades té una fila a dalt indicant el nom de cada columna/variable.
sep=“;” indica quin és el separador entre columnes de dades. Depenent de com tingueu configurat l’Excel pot ser “;” o “,”.

Podeu accedir a l’ajuda d’una funció (e.g. read.csv()), per veure quins arguments té, picant a la consola:

help(read.csv)

?read.csv

També podeu picar sobre la funció a l’script i després picar F1.

1.3.2 format dades d’R .RData

Si teniu dades en format R, les reconeixereu per l’extensió .RData, també les podeu carregar molt fàcilment.

Passos a seguir:

Guardeu l’arxiu al vostre directori de treball o en una carpeta coneguda²
Escriviu load i entre () el nom de l’arxiu. Recordeu posar el nom entre cometes

load("pes_alumnes.RData")

I ja ho teniu carregat. Si s’ha carregat us apareixerà a la Pestanya d’Environment de l’RStudio (part superior dreta), sota el títol Data

1.4 Comprovar que s’ha fet correctament

Una vegada han estat carregades i assignades les dades cal comprovar que no ens hem equivocat. Existeixen diverses formes de saber si les dades han estat carregades correctament. A banda de mirar les característiques de l’objecte a la finestra superior dreta d’RStudio ,com hem comentat abans, des de la consola R tenim diverses opcions.

Podem usar la força bruta i picar (o bé copiar i enganxar) el nom de l’objecte (en aquest cas un data frame) pes_alumnes.data a la consola:

pes_alumnes.data

##    weight height shoe.size exercise diet sex
## 1    42.0    156      35.0      2.0    3   D
## 2    44.0    155      35.0      5.0    1   D
## 3    45.0    154      35.0      5.0    2   D
## 4    45.0    157      37.0      4.0    2   D
## 5    45.0    156      37.0      2.0    0   D
## 6    48.0    160      37.0      7.0    2   D
## 7    48.0    165      39.0      3.0    1   D
## 8    48.0    166      39.0      3.0    1   D
## 9    49.0    164      39.0      2.0    1   D
## 10   50.0    170      39.0      2.0    2   D
## 11   50.0    162      36.0      0.0    3   D
## 12   51.0    160      38.0      8.0    1   D
## 13   51.0    167      37.0      1.0    6   D
## 14   51.5    164      38.0      2.0    0   D
## 15   52.0    168      37.0      5.0    1   D
## 16   52.0    165      39.0      4.0    1   D
## 17   52.0    163      39.0      3.0    2   D
## 18   53.0    162      36.0      0.0    5   D
## 19   54.0    170      40.0      0.0    0   H
## 20   54.0    160      37.0      4.0    3   D
## 21   55.0    169      39.0      1.0    2   D
## 22   55.0    156      37.0      2.0    0   D
## 23   55.0    168      39.0      7.0    0   D
## 24   56.0    166      37.0      5.0    3   D
## 25   56.0    170      40.0      0.0    0   D
## 26   57.0    160      39.0      5.0    0   D
## 27   58.0    175      40.0      4.0    1   D
## 28   58.0    161      41.0      3.0    5   H
## 29   58.0    172      38.0      5.0    0   D
## 30   58.5    175      43.0      5.3    2   H
## 31   60.0    170      39.0      6.0    2   D
## 32   60.0    178      42.0      2.0    4   H
## 33   60.0    161      38.0      0.0    0   D
## 34   60.0    173      43.0      3.0    0   H
## 35   60.0    178      42.0      4.0    2   D
## 36   61.0    168      39.0      2.0    0   D
## 37   61.5    179      41.5      6.0    2   D
## 38   62.0    165      38.0      5.0    0   D
## 39   62.0    170      39.0      5.0    0   D
## 40   63.0    175      40.5      9.0    3   H
## 41   63.0    167      37.0      0.0    0   D
## 42   63.0    180      39.0      5.0    0   D
## 43   63.0    169      38.0      4.5    0   D
## 44   64.0    174      40.0      5.0    1   D
## 45   64.0    174      41.0      5.0    3   D
## 46   65.0    170      41.0      5.0    0   D
## 47   65.0    170      40.0      3.0    0   D
## 48   65.0    173      42.0      8.0    0   H
## 49   65.0    169      40.0      0.0    3   H
## 50   66.0    174      42.0      9.0    0   H
## 51   66.0    179      45.0      5.0    4   H
## 52   66.0    175      42.0      4.0    5   H
## 53   67.0    190      43.0      6.0    2   H
## 54   67.0    175      42.0      0.0    1   H
## 55   67.0    173      40.0      3.5    4   D
## 56   67.0    182      44.0      6.0    2   H
## 57   67.0    175      43.0      3.0    3   H
## 58   68.0    160      37.0      1.0    3   D
## 59   69.0    173      41.0      4.0    0   H
## 60   70.0    181      43.0      6.0    2   H
## 61   70.0    172      41.0      0.0    0   H
## 62   70.0    170      43.0      8.0    3   H
## 63   70.0    160      39.0      2.0    0   D
## 64   70.0    155      36.0      2.0    1   D
## 65   70.0    176      43.0     15.0    0   H
## 66   70.0    158      38.0      4.0    1   D
## 67   70.0    171      40.0      5.0    0   D
## 68   70.0    172      41.0      0.0    2   H
## 69   70.0    174      42.0      5.0    0   H
## 70   71.0    167      42.0      0.0    0   H
## 71   72.0    182      43.0      3.0    0   H
## 72   72.0    165      38.0      3.0    0   D
## 73   72.0    179      43.0      0.0    0   H
## 74   73.0    187      43.0      3.0    0   H
## 75   74.0    176      41.0      0.0    0   H
## 76   74.0    178      45.0      7.0    2   H
## 77   75.0    167      39.0      0.0    0   D
## 78   75.0    178      42.0      0.0    0   H
## 79   75.0    184      44.0      4.0    2   H
## 80   75.0    170      40.0      2.0    0   D
## 81   76.0    177      44.0     15.0    7   H
## 82   76.0    170      41.0      0.0    0   H
## 83   77.0    182      44.0      3.0    1   H
## 84   78.0    185      44.0      4.0    1   H
## 85   78.0    182      44.0      0.0    0   H
## 86   78.0    190      47.0      4.0    0   H
## 87   79.0    165      40.0      0.0    4   D
## 88   80.0    174      43.0      1.0    0   H
## 89   80.0    194      45.0      2.0    0   H
## 90   80.0    183      45.0      2.0    2   H
## 91   80.0    179      42.0      5.0    0   H
## 92   82.0    168      39.0      3.0    0   D
## 93   82.0    183      46.0      8.0    2   H
## 94   82.0    174      42.0      2.0    6   H
## 95   85.0    173      41.0      1.0    2   D
## 96   85.0    176      43.0      0.0    0   H
## 97   85.0    183      46.0      3.0    0   H
## 98   90.0    180      45.0      0.0    1   H
## 99  101.0    173      43.0      7.0    0   H

També podem seleccionar l’objecte a l’script i picar CTR + Intro. Aquesta opció mostra l’objecte (i.e. taula de dades) complet però no resulta massa pràctic quan les taules de dades són grans.

La taula de dades (data.frame) pes_alumnes.data conté 6 columnes anomenades:

pes
height
shoe.size
exercici
dieta
sex

i 99 dades.

Per veure solament al nom de les columnes/variables de la taula podem usar la funció names():

names(pes_alumnes.data)

## [1] "weight"    "height"    "shoe.size" "exercise"  "diet"      "sex"

No obstant però, aquesta opció sols ens indica el nom de les columnes/variables i no sabem si les dimensions de la taula són les correctes. Tampoc veiem el contingut/natura de la informació de cada columna. Per a conèixer les dimensions podem usar la funció dim():

dim(pes_alumnes.data)

## [1] 99  6

El primer valor correspon al nombre de files (99) i el segon de columnes (6). Es recomanable comprovar a l’Excel (dades_base.csv) si aquestes són les dimensions que ha de tindre la taula de dades (recordeu que l’encapçalament no conta com a línia una vegada carregada la taula a R). El nom de les columnes i les dimensions de la taula poden ser correctes però per a saber si la natura de la informació que hi conté és la correcta cal usar altres funcions. Entre les diferents opcions disponibles summary() i str() són les més habituals.

str(pes_alumnes.data)

## 'data.frame':	99 obs. of  6 variables:
##  $ weight   : num  42 44 45 45 45 48 48 48 49 50 ...
##  $ height   : int  156 155 154 157 156 160 165 166 164 170 ...
##  $ shoe.size: num  35 35 35 37 37 37 39 39 39 39 ...
##  $ exercise : num  2 5 5 4 2 7 3 3 2 2 ...
##  $ diet     : int  3 1 2 2 0 2 1 1 1 2 ...
##  $ sex      : Factor w/ 2 levels "D","H": 1 1 1 1 1 1 1 1 1 1 ...

summary(pes_alumnes.data)

##      weight           height        shoe.size        exercise     
##  Min.   : 42.00   Min.   :154.0   Min.   :35.00   Min.   : 0.000  
##  1st Qu.: 56.50   1st Qu.:165.0   1st Qu.:39.00   1st Qu.: 1.500  
##  Median : 66.00   Median :171.0   Median :40.00   Median : 3.000  
##  Mean   : 65.35   Mean   :171.2   Mean   :40.52   Mean   : 3.508  
##  3rd Qu.: 73.50   3rd Qu.:176.5   3rd Qu.:43.00   3rd Qu.: 5.000  
##  Max.   :101.00   Max.   :194.0   Max.   :47.00   Max.   :15.000  
##       diet       sex   
##  Min.   :0.000   D:53  
##  1st Qu.:0.000   H:46  
##  Median :1.000         
##  Mean   :1.343         
##  3rd Qu.:2.000         
##  Max.   :7.000

En general str() dóna una informació més acurada de l’estructura de la informació i no tant del contingut, però per a nouvinguts pot ser difícil d’interpretar. És per això que, en el context d’AAD, la combinació de dim() i summary() és més recomanable.

La funció summary dóna una descripció de la informació continguda a cada columna que permet deduir la seua natura o classe. Quan per a una columna dona estadístics de posició com el mínim (Min.), el primer quantil (1st Qu.) etc. podem inferir que es tracta d’una variable contínua. Mentre que quan dóna un compte del nombre de vegades que apareix una o més categories podem assumir que és una variable categòrica o factor.

Com veurem a continuació, és molt importat que la variable que defineix els grups que volem comparar siga un factor o l’anàlisi serà incorrecta!!!

1.5 Exploració de les dades

1.5.1 Selecció de dades/files/columnes

Com podeu veure a continuació existeixen diverses formes d’accedir a una única columna d’un data.frame.

La primer usa el nom de l’objecte (data.frame) i el nom de la columna. Entremig cal posar el símbol del dòlar ($).
La segona també usa el nom de l’objecte (data.frame) i el nom de la columna però es basa en l’estructura de les dades. Els símbols [_ i _] amb una coma entremig (,) fan referència a l’estructura de files i columnes de la taula de dades. L’espai abans de la coma ens permet accedir files determinades mentre que el segon ho permet per columnes. Com que les columnes tenen noms, pes_alumnes.data[ ,“sex”] ens permet accedir a la columna sex.
Finalment la tercera forma no usa el nom de la columna i es basa en l’ordre de les columnes. sex és la sisena columna de la taula de dades. Llavors, pes_alumnes.data[ , 6] permet accedir a aquesta columna.

pes_alumnes.data$sex

##  [1] D D D D D D D D D D D D D D D D D D H D D D D D D D D H D H D H D H D
## [36] D D D D H D D D D D D D H H H H H H H D H H D H H H H D D H D D H H H
## [71] H D H H H H D H H D H H H H H H D H H H H D H H D H H H H
## Levels: D H

pes_alumnes.data[ ,"sex"]

##  [1] D D D D D D D D D D D D D D D D D D H D D D D D D D D H D H D H D H D
## [36] D D D D H D D D D D D D H H H H H H H D H H D H H H H D D H D D H H H
## [71] H D H H H H D H H D H H H H H H D H H H H D H H D H H H H
## Levels: D H

pes_alumnes.data[ , 6]

##  [1] D D D D D D D D D D D D D D D D D D H D D D D D D D D H D H D H D H D
## [36] D D D D H D D D D D D D H H H H H H H D H H D H H H H D D H D D H H H
## [71] H D H H H H D H H D H H H H H H D H H H H D H H D H H H H
## Levels: D H

Com a exercici vos recomanem que piqueu:

pes_alumnes.data[1, ]

pes_alumnes.data[1, 6]

pes_alumnes.data[1:3, ]

pes_alumnes.data[ ,5:6]

pes_alumnes.data[1:3, 5:6]

Quines dades es mostren?

També es poden seleccionar files i columnes en base als valors d’una altra fila o columna:

pes_alumnes.data[pes_alumnes.data$sex=="H", c("height", "sex")]

##    height sex
## 19    170   H
## 28    161   H
## 30    175   H
## 32    178   H
## 34    173   H
## 40    175   H
## 48    173   H
## 49    169   H
## 50    174   H
## 51    179   H
## 52    175   H
## 53    190   H
## 54    175   H
## 56    182   H
## 57    175   H
## 59    173   H
## 60    181   H
## 61    172   H
## 62    170   H
## 65    176   H
## 68    172   H
## 69    174   H
## 70    167   H
## 71    182   H
## 73    179   H
## 74    187   H
## 75    176   H
## 76    178   H
## 78    178   H
## 79    184   H
## 81    177   H
## 82    170   H
## 83    182   H
## 84    185   H
## 85    182   H
## 86    190   H
## 88    174   H
## 89    194   H
## 90    183   H
## 91    179   H
## 93    183   H
## 94    174   H
## 96    176   H
## 97    183   H
## 98    180   H
## 99    173   H

pes_alumnes.data[pes_alumnes.data$sex=="H", c(1, 6)]

##    weight sex
## 19   54.0   H
## 28   58.0   H
## 30   58.5   H
## 32   60.0   H
## 34   60.0   H
## 40   63.0   H
## 48   65.0   H
## 49   65.0   H
## 50   66.0   H
## 51   66.0   H
## 52   66.0   H
## 53   67.0   H
## 54   67.0   H
## 56   67.0   H
## 57   67.0   H
## 59   69.0   H
## 60   70.0   H
## 61   70.0   H
## 62   70.0   H
## 65   70.0   H
## 68   70.0   H
## 69   70.0   H
## 70   71.0   H
## 71   72.0   H
## 73   72.0   H
## 74   73.0   H
## 75   74.0   H
## 76   74.0   H
## 78   75.0   H
## 79   75.0   H
## 81   76.0   H
## 82   76.0   H
## 83   77.0   H
## 84   78.0   H
## 85   78.0   H
## 86   78.0   H
## 88   80.0   H
## 89   80.0   H
## 90   80.0   H
## 91   80.0   H
## 93   82.0   H
## 94   82.0   H
## 96   85.0   H
## 97   85.0   H
## 98   90.0   H
## 99  101.0   H

En aquest cas s’han seleccionat totes les dades de les columnes “height” i “sex” que tenen com a sex la categoria “H” (Home).

Alerta! Fixeu-vos que “height” i “sex” i 2 i 3 estan continguts dins d’un vector de dades c(,). De no ser així la selecció no funcionarà perquè poden haver taules de dades amb més de dues dimensions. Llavors pes_alumnes.data[, 2, 6] indica tots els elements en la dimensió 1, la columna 1 de la dimensió 2 i la 6 en la dimensió 3. També és molt important que la categoria vaja entre " " o R tractarà de trobar un objecte al Global Environment de nom H.

1.5.2 Exploració gràfica

La funció summary() ens dóna informació sobre el contingut i natura de cada variable. No obstant però, és habitual realitzar una exploració gràfica de cada variable per separat i de la relació entre variables. L’opció bàsica per a la darrera anàlisi consisteix a realitzar un gràfic, mitjançant la funció plot(), amb dues variables d’interès per observar la presència d’alguna mena de relació entre ambdues:

plot(x=pes_alumnes.data$height,
     y=pes_alumnes.data$shoe.size,
     col="black",
     xlab="height (cm)", ylab="shoe.size (#)",
     bty="n")

Pel que fa als arguments:

x és la variable usada a l’eix x.
y és la variable usada a l’eix y.
col determina color dels punts, línies, etc.
xlab determina l’etiqueta (label) usat per a l’eix x.
ylab determina l’etiqueta (label) usat per a l’eix y.
bty elimina la caixa que envolta el plot. És un argument merament estètic.

Pel que fa a la relació entre les dues variables, el gràfic mostra una correlació clara entre l’height dels i les alumnes i la mida de shoe.size que tenen.

El gràfic anterior mostra les dues categories emprades juntes. És possible diferenciar-les usant diferents colors i/o tipus de punts per a cadascuna.

plot(x=pes_alumnes.data$height,
     y=pes_alumnes.data$shoe.size,
     col=ifelse(pes_alumnes.data$sex=="D", "orangered", "dodgerblue2"),
     pch=ifelse(pes_alumnes.data$sex=="D", 16, 1),
     xlab="height (cm)", ylab="shoe.size (#)", bty="n")

legend("topleft",
       legend=c("Women (D)", "Men (H)"),
       col=c("orangered", "dodgerblue2"),
       pch=c(16, 1), 
       bty="n")

Figure 1.2: Scatter plot entre el height i la mida de shoe.size diferenciat per sex (binari)

Pel que fa als arguments:

col determina el color dels punts. En aquest cas s’ha usat la funció ifelse() per a modificar els valors segons els valors de la variable sex. ifelse(pes_alumnes.data$sex==“D”, “orangered”, “dodgerblue2”) pot llegir-se com: si (if) pes_alumnes.data$sex és (==) “D” (dona) usa el color “orangered” (carabassa) o bé usa “dodgerblue2” (blau).
pch determina el tipus de punt utilitzat. ifelse(pes_alumnes.data$sex==“D”, 16, 1) pot llegir-se de forma anàloga com: si (if) pes_alumnes.data$sex és (==) “D” (dona) usa 16 (punts sòlids) o bé usa 1 (buits).

A continuació es mostra un exemple de l’ús de la funció legend(). Serveix per a afegir la llegenda a la cantonada superior esquerre (“topleft”).

Pel que fa als arguments:

“topleft” indica la posició.
legend el text/etiquetes.
col el color dels punts.
pch determina el tipus de punt utilitzat.
bty elimina la caixa que envolta el plot. És un argument merament estètic.

Com que sols hi ha dues etiquetes col té dos valors (“orangered” i “dodgerblue2”) i pch també. Fixeu-vos que tots dos es troben indicats mitjançant un vector (c( , )).

La resta d’arguments ha estat explicada al plot anterior.

Quan la variable x és categòrica (pes_alumnes.data$sex), R automàticament genera un boxplot per a cada categoria.

plot(x=pes_alumnes.data$sex,
     y=pes_alumnes.data$height,
     col=c("orangered", "dodgerblue2"),
     xlab="sex", ylab="height (cm)", frame=F)

legend("topleft",
       legend=c("Women (D)", "Men (H)"),
       fill=c("orangered", "dodgerblue2"), 
       col=c("black"),
       bty="n")

La funció boxplot() permet generar el mateix tipus de plot. Mireu l’ajuda (help(boxplot)) i tracteu de reproduir la mateixa figura usant la funció boxplot().

Els arguments han estat explicats als plots anteriors.

1.5.3 Exploració numèrica

De vegades és necessari conèixer el tipus de distribució que tenen les dades. En moltes ocasions les dades han de ser normals per tal de complir amb els supòsits de les anàlisis estadístiques. Per tal d’investigar si una sèrie de dades s’ajusta a una distribució normal es pot aplicat el test de normalitat de Shapiro-Wilk usant la funció shapiro.test().

shapiro.test(pes_alumnes.data$height)

## 
## 	Shapiro-Wilk normality test
## 
## data:  pes_alumnes.data$height
## W = 0.98816, p-value = 0.5278

hist(pes_alumnes.data$height,
     main="Histograma del height",
     xlab="height (cm)",
     col="dodgerblue3",
     bty="n")

Figure 1.3: Histograma del height dels alumnes

Com que el p-value > 0.05 es pot concloure que la distribució de les alçades és normal. L’histograma de les dades ja ens indicava que amb tota probabilitat el test anava a donar aquest resultat.

En altres ocasions podem estar interessats en la relació entre dues variables continues (e.g. correlació). Per tal de cheightular la correlació entre dues variables la funció més emprada és cor().

cor(x = pes_alumnes.data$height,
    y = pes_alumnes.data$shoe.size,
    method = "pearson")

## [1] 0.8680414

Com era d’esperar la correlació entre l’height dels i les alumnes i la seua mida de shoe.size és molt alta i positiva (i.e. 0.8680414). Es a dir, quan més alt és l’individu, major és la mida de shoe.size. Quan les variables no són normals és millor usar l’índex de correlació Spearman (i.e. method = “spearman”).

Pel que fa als arguments:

x és la 1a variable.
y és la 2a variable.
method és el mètode emprat per cheightular la correlació: “pearson”, “kendall” o “spearman”.

La funció cor() en dóna el valor de la correlació però no indica la significació estadística (i.e. credibilitat) d’aquest valor. Per tal d’obtindre-la podem usar la funció cor.test().

cor.test(x = pes_alumnes.data$height,
         y = pes_alumnes.data$shoe.size,
         method = "pearson",
         alternative = "two.sided",
         conf.level = 0.95)

## 
## 	Pearson's product-moment correlation
## 
## data:  pes_alumnes.data$height and pes_alumnes.data$shoe.size
## t = 17.219, df = 97, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8093152 0.9095852
## sample estimates:
##       cor 
## 0.8680414

Pel que fa als arguments:

x és la 1a variable.
y és la 2a variable.
method és el mètode emprat per cheightular la correlació: “pearson”, “kendall” o “spearman”.
alternative indica la hipòtesi alternativa.
conf.level determina el nivell de confiança.

1.6 Adequació de dades

Com veurem a continuació, en alguns casos cal que les variables tinguen una classe específica per a realitzar certes anàlisis. La funció summary() ja ens ha indicat que la variable sex és categòrica. No obstant, podem veure quina classe té usant la funció class(). També podem preguntar si la columna/variable conté factors usant la funció is.factor().

class(pes_alumnes.data$sex)

## [1] "factor"

is.factor(pes_alumnes.data$sex)

## [1] TRUE

La resposta és TRUE (cert) per la qual cosa sabem que la variable és categòrica (factor).

Si calgués convertir una variable (e.g. continua) a categòrica (factor) s’usa la funció factor().

factor(c(1:5))

## [1] 1 2 3 4 5
## Levels: 1 2 3 4 5

1.7 Principis bàsics de la regressió lineal (RLS)

En ocasions a més del valor de la correlació estem interessats a matematitzar aquesta relació. El model resultant permet explorar numèricament la relació entre ambdues variables obtenint estadístics diferents. També permet fer prediccions sobre valors no observats. Un dels models més simples i habituals és la regressió lineal entre dues variables on s’obtenen els paràmetres corresponents a una fórmula de tipus:

\[y = A \cdot x + B\]

on A és el pendent de la recta i B la intercepció o ordenada a l’origen. El signe d’A pot esser interpretat d’una forma anàloga al signe de la correlació (i.e. relació positiva o negativa) entre les dues variables.

Continuant amb l’exemple anterior, la relació entre l’height (x/variable independent) i la mida de shoe.size dels (y/variable dependent) i les alumnes és positiva i correspon a la línia roja que es pot observar a la figura següent. Com es veurà a la pràctica corresponent, els models lineals han de complir amb una sèrie de supòsits per tal d’ésser considerats vàlids. Una de molt important es bassa en la Homoscedasticitat, Normalitat i manca de no-linealitat dels residus. Els residus corresponen a la diferència entre el valor modelitzat/predit (línia roja) i l’observat (fletxes blaves a la figura següent). Per simplicitat la figura sols inclou 5 fletxes però hi ha tants residus com dades (en aquest cas 99).

Figure 1.4: A) regressió entre height i la mida de shoe.sizes dels alumnes; B) residus de la RLS

A banda dels tests estadístics que s’apliquen per comprovar l’Homoscedasticitat i la Normalitat és habitual fer un gràfic amb els residus. Quan observem alguna mena de patró als residus, ens indica que el model és incorrecte i que cal fer alguna cosa per millorar-lo (e.g. transformació de les dades). A la figura següent es mostren diferents patrons idealitzats. Únicament la figura A mostra un patró en els residus correcte. Les figures B, E i F mostren biaix, les figures C i F mostren no-linealitat i les figures D, E i F mostren Heteroscedasticitat (no Homoscedasticitat): la variància dels residus és major per als valors superiors del rang de la variable en comparació amb els menors. Fixeu-vos que poden haver combinacions (e.g. la figura E mostra no linealitat i Heteroscedasticitat en els residus).

Figure 1.5: Patrons desitjats (verds) i incorrectes (roig) als residus

Tot i que les causes d’eixos patrons poden ser diverses. Per exemple:

B i C poden ser causats per la manca d’algun predictor. Existeixen models multivariants del tipus: \[z = A \cdot x + B \cdot y + C\]
B també pot ocórrer per causa d’una mala optimització dels paràmetres del model (i.e. A).
a banda C pot estar per una mala parametrització del model ja que la relació s’ajusta més a una funció quadràtica.

1.8 Com podem fer una ANOVA

L’anàlisi de la variància (ANOVA) s’utilitza per analitzar les diferències entre les mitjanes d’un seguit de grups (i.e. 2 o 2+ grups). A l’exemple que anem a treballar farem una ANOVA per saber si hi ha diferències en l’height en funció del sex (binari). Com que no tenim el mateix # de casos entre les dues categories de la variable categòrica o factor (i.e. sexe (binari) dels alumnes), per tal de realitzar l’anàlisi d’ANOVA utilitzem la funció lm() (linear model):

anova_height <- lm(height~sex, data=pes_alumnes.data)

Pel que fa als arguments:

anova_height és l’objecte que contindrà l’anàlisi.
height és la variable depenent (y).
~ és el signe que indica la relació o =.
sex és la variable independent (x).
data=pes_alumnes.data indica d’on s’han de extraure les variables (i.e. columnes de dades).

Amb certes particularitats com ara que la variable independent és categòrica (i.e. factor), cal acomplir amb els mateixos supòsits que han estat descrits per a la regressió lineal. Des d’un punt de vista gràfic les distribucions de valors semblen diferents (A - ANOVA). Així els valors modelitzats (i.e. mitjanes) per a cada categoria (i.e. D i H) amb el model de l’ANOVA també semblen diferents (quadrats carabassa i blau) (A - ANOVA). La variància dels residus de totes dues categories (fletxes blaves) són semblants (B - RESIDUS), hi ha Homoscedasticitat, i la distribució dels residus és normal perquè al C - Normal Q-Q Plot els residus s’alineen al llarg de la línia teòrica dels quantils que tindria una distribució normal amb la mitjana i desviació estàndard dels residus de l’ANOVA.

Figure 1.6: A) ANOVA (regressió) entre el gènere (binari) i l’height dels i les alumnes. B) residus de l’ANOVA (regressió). C - Normal Q-Q Plot.

Des d’un punt de vista numèric devem inspeccionar els resultats mitjançant la funció summary(). Ens hem de fixar en els coeficients (Coefficients:). El valor estimat (Estimate) per a la intercepció ((Intercept)) correspon a la mitjana dels valors que alhora correspon al valor predit pel model. El valor calculat és de 165.8491 aquest correspon al punt carbassa de la figura anterior. A la columna Pr(>|t|) podem trobar la significació d’aquest paràmetre. El valor (< 2e-16) i el nombre d’estreles (_***_) ens indiquen que el valor es significatiu i distint de zero.

El valor estimat (Estimate) per a la categoria H (sexH) + el valor calculat per a la categoria D ((Intercept)) correspon a la mitjana dels valors d’H que alhora correspon al valor predit pel model (punt blau de la figura anterior = 165.8491 + 11.4988).

Des del punt de vista de l’objectiu de l’anàlisi, a la columna Pr(>|t|) podem trobar la significació d’aquest paràmetre. El valor (2.93e-14) i el nombre d’estreles (_***_) ens indiquen que el valor es significatiu i distint de zero. El supòsits de l’anàlisi han estat comprovats a la figura anterior per la qual cosa podem concloure que hi ha diferències significatives entre l’alçada (height) de les persones classificades com a Dones (D) i Hòmens (H) i que aquesta té una mitjana d’11.4988 cm.

summary(anova_height)

## 
## Call:
## lm(formula = height ~ sex, data = pes_alumnes.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.3478  -4.3478  -0.3478   4.1509  16.6522 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 165.8491     0.8791 188.654  < 2e-16 ***
## sexH         11.4988     1.2897   8.916 2.93e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.4 on 97 degrees of freedom
## Multiple R-squared:  0.4504,	Adjusted R-squared:  0.4447 
## F-statistic: 79.49 on 1 and 97 DF,  p-value: 2.935e-14

Com heu observat al llarg de la pràctica programar figures en R no és gaire amigable. R a més d’esser una eina d’anàlisi és un llenguatge de programació. Llavors, quan hem de repetir una anàlisi moltes vegades és recomanable crear una funció que agafant com a argument el nostre model gènere els plots de diagnosi de forma automatitzada.

En aquest exemple anem a generar una funció anomenada mcheck() que tindrà com a argument el nostre model (el.meu.model) i que generarà tots dos plots de diagnosi.

 mcheck <- function(el.meu.model){                                                                                                                                              
   residus <- el.meu.model$resid # crea un objecte anomenat 'residus' amb els residus del meu model 
   fitted.values <- el.meu.model$fitted # crea un objecte anomenat 'fitted.values' amb els valors modelitzats pel meu model  
   
   par(mfrow=c(1, 2), bty="n") # paràmetres del plot. mfrow=c(1,2) indica dos plots en una única fila. bty="n" elimina la caixa que envolta el plot. És un argument merament estètic.
   
   plot(fitted.values, residus,
        xlab="Fitted values", ylab="Residuals") # plot fitted.values vs. residus
   abline(h=0, lty="dashed") # genera un línia horitzontal de valor constant 0 (h=0). lty="dashed" indica que la línia ha de ser de punts.
   
   qqnorm(residus, 
          xlab="Normal scores", ylab="Ordered residuals") # funció que genera un 'Q-Q plot normal'.
   qqline(residus, lty="dashed") # funció que genera la línia que tindria una distribució un 'Q-Q plot normal'.
  
   }

Ara podem executar la funció i veure els resultats.

mcheck(el.meu.model=anova_height)

Figure 1.7: Residus de l’ANOVA (esquerra). C - Normal Q-Q Plot (dreta).

A la propera pràctica contrastos no-paramètrics veurem com realitzar una comprovació numèrica dels supòsits de l’ANOVA mitjançant testos estadístics. També veurem quins proves cal utilitzar quan no es compleixen els supòsits de per a realitzar una anàlisi de l’ANOVA.

L’anàlisi de la variància (ANOVA) s’utilitza per analitzar les diferències entre les mitjanes d’un seguit de grups (i.e. 2 o 2+ grups). En la seva forma més simple, ANOVA proporciona una prova estadística indicant si dues o més mitjanes de població són iguals, i per tant generalitza la prova t més enllà de dues mitjanes.↩︎
Si guardeu l’arxiu en una carpeta coneguda, després per RStudio el podeu anar a buscar. Mirant a la part inferior dreta, a la pestanya Files podeu navegar per els vostres arxius. Un cop localitzades les dades que voleu carregar, si cliqueu 2 vegades amb el ratolí, es carreguen automàticament les dades. Ull! si ho feu així, no queda registrat a l’Script d’R, i després podeu tenir problemes per saber quines dades vau utilitzar per fer aquelles anàlisis.↩︎