Introduzione a R

Autore/Autrice

Affiliazione

Ottavia M. Epifania ottavia.epifania@unitn.it

Università di Trento

Introduzione

Prima di iniziare: Scaricate la cartella compressa data dentro cui si trovano tutti i dataset che sono usati in questo file.

R è un software open source per il calcolo statistico, la grafica e molto altro

Installazione

È necessario installare R (disponibile https://cran.r-project.org/bin/windows/base/)

RStudio è l’IDE perfetto per R $\rightarrow$consente un utilizzo di R migliore e più semplice

RStudio si può installare dopo aver installato R (disponibile https://posit.co/download/rstudio-desktop/)

R funziona su Windows, MacOs, Unix

RStudio IDE

La schermata di RStudio è suddivisa in 4 pannelli:

Alto sx: è l’area dedicata agli script di codice
Alto dx: l’area dove è possibile vedere quello che è presente all’interno dell’ambiente di R (i.e., tutti gli oggetti, le variabili, i dataframe che sono stati caricati/creati per le analisi)
Basso sx: la console di R, dove vengono mostrati gli output delle operazioni e i codici eseguiti
basso dx: l’area dei grafici

Console

I comandi nella console vengono eseguiti ma il codice non viene salvato e scompare sequenzialmente.

Possono essere recuperati utilizzando la freccia verso l’alto.

Per eseguire qualsiasi comando nella console $\rightarrow$ Invio (o Enter)

L’output appare immediatamente nella console.

Script

Si possono salvare gli script con tutte e righe di codice che servono per le analisi

Per far “correre” qualsiasi comando nello script $\rightarrow$ Ctrl + Invio (or ctrl + Enter or cmd + Enter) oppure utilizza semplicemente il pulsante “Run” in cima allo script.

L’output appare nella console.

Si possono inserire commenti (i.e., righe di codice che non vengono eseguite) inserendo # $\rightarrow$ tutto il testo dopo # non viene eseguito

Per fare lo switch alla console $\rightarrow$ ctrl + 2

Per fare lo switch allo script $\rightarrow$ ctr + 1

Working directories

Bisogna conoscere le working directories (da qui in avanti indicate come wd), ovvero le cartelle dove si lavor, dove si salvano i dati ecc.

getwd() # per sapere in quale wd si sta lavorando

dir() # lista degli oggetti all'interno della wd

Per cambiare la wd:

setwd("C:/Users/huawei/OneDrive/Documenti/GitHub/RcouRse")

Basics

Calcolat`R`ice

3 + 2   # più
3 - 2   # meno
3 * 2   # per
3 / 2   # diviso
sqrt(4) # radice quadrata
log(3)  # logaritmo naturale
exp(3)  # esponenziale

Parentesi e simili:

() [] {} "" : ; ,

L’uso delle parantesi funziona come se si stesse scrivendo un’equazione:

(3 * 2)/ sqrt(25 + 4) # Look at me!

R ignora tutto quello dopo # (è un commento)

Assign

I risultati delle operazioni possono essere “memorizzati” in oggetti (variabili, in senso informatico) con nomi specifici definiti dagli utenti

Per assegnare un valore a un oggetto, si possono utilizzare due operatori:

x = exp(2^2)
X <- log(2^2)

Gli elementi a destra del simbolo vengono assegnati all’oggetto a sinistra

Attenzione! R e case sensitive: x e X sono due oggetti diversi!!!

Nomi delle variabili

Nomi di variabili validi sono lettere, numeri, punti, underscore (e.g., variable_name)

I nomi delle variabili non possono iniziare con numeri

L’idea per i nomi delle variabili è di dare indicazione di cosa contengono senza stare a scrivere divine commedie, ad esempio

Buoni nomi per le variabili:

mean_rt
sum_age
prop_item

Pessimi nomi per variabili:

x
y
ciccio

Chiedere aiuto

R community is the best feature of R

Basta copiare & incollare qualsiasi messaggio di errore o avviso in Google o chiedere a Google “come fare [qualcosa] in R”.

Chiedete a R di aiutarvi! Basta scrivere ? nella console seguito dal nome della funzione che si vuole conoscere:

?mean()

Mostra la help page della funzione mean()

Ordine, ordine e ancora ordine

R projects for the win

Gestire i directory di lavoro è un vero e proprio incubo. Per evitare di doverle gestire, si possono usare i porgetti di R

Un progetto R crea automaticamente una sorta di isola sul vostro pc con le sue sottoisole (che sono le sottocartelle che create voi)

Sono davvero comodi almeno per due motivi:

Permettono di avere diverse istanze di R aperte allo stesso tempo $\rightarrow$ Potete lavorare a progetti diversi contemporaneamente
Tutti i file sono sempre a disposizione, senza preoccuparsi delle wd

Creare un nuovo `R` project

File $\rightarrow$ New project e scegliete quello che fa al caso vostro (a meno che non abbiate già inizializzato un directory per il tuo progetto, selezionate una nuova directory):

R project “basic”
R package
Shiny project

e tanto altro

Creare un nuovo R project

File $\rightarrow$ New project:

Selezionare il tipo di progetto

Assegnargli una directory e dargli un nome

Buttare via la spazzatura

L’ambiente R dovrebbe essere sempre ordinato

Per pulirlo:

ls() # lista tutti gli oggetti nell'ambiente
rm(A) # rimuove l'oggetto A dalla directory
rm(list=ls()) # rimuove tutto dall'ambiente

Salvare l’ambiente

Per salvare tutti i risultati delle analisi in formato ambiente di R:

save.image("my-computations.RData")

Per ricaricare l’ambiente in R:

load("my-computations.RData")

Attenzione!

I codici di prima, dove viene salvato e ricaricato l’ambiente, presuppone che l’ambiente di R venga salvato nella stessa wd dove è presente lo script

Quando salvare l’ambiente?

I calcoli sono lenti e hai bisogno che siano sempre e facilmente accessibili

La pratica migliore è salvare lo script e documentarlo in un file RMarkdown $\rightarrow$ Riproducibilità!

Your turn!

Crea un progetto R per questo corso nella tua cartella “documenti” (scegli un nome carino :)
Crea un nuovo script R (shift + ctrl + n)
Calcolatrice utilizzando lo script:
- $\sqrt{15} \times 14 - \frac{22}{4}$ [48.72177]
- $\frac{\sqrt{7-\pi}}{3 \times (45-34)}$ [0.05952372]
Salva lo script
Assegna i risultati del primo calcolo a una variabile chiamata my_results

See solution

# Calculator
sqrt(15) *14 - 22/4

[1] 48.72177

# Calculator 
sqrt(7 - pi)/(3*(45 -34))

[1] 0.05952372

# assign
my_results = sqrt(15) *14 - 22/4

Strutture in `R`

Funzioni e loro argomenti

Quasi tutto si fa attraverso le funzioni, composte da:

Un nome: mean
un paio di parentesi: ()
degli argomenti: na.rm = TRUE

mean(1:5, trim = 0, na.rm = TRUE)

[1] 3

Gli argomenti possono essere impostati con valori predefiniti; la funzione di ogni argomento è spiegata esaustivamente nell’help della funzione stessa, e.g., ?mean()

Gli argomenti possono essere specificati:

senza un nome (seguendo il loro ordine di default)
con il loro nome (in qualsiasi ordine) $\rightarrow$ keyword matching

mean(x, trim = 0.3, na.rm = TRUE)

Alcune funzioni funzionano anche senza argomenti:

ls(), dir(), getwd()

Per vedere il codice di una funzione, basta scrivere il suo nome senza parentesi nella console e premere invio

chisq.test

`concatenate`, `c()`

Una delle funzioni più usate di R.

Concatena diverse variabili insieme per creare una nuova unica variabile

x = c(1, 2, 3) # crea un vettore concatenando 1,2,3
x

[1] 1 2 3

X = 1:3 # crea lo stesso identico vettore
X

[1] 1 2 3

x == X

[1] TRUE TRUE TRUE

Vettori

Tipi di variabili diverse $\rightarrow$ Tipi di vettori diversi:

int: Numeri interi (1,2, 0, -5)
num: Numeri (-0.56, 3, 4.5)
logi: Logico (TRUE vs FALSE)
chr: characters ("a", "b", "c")
factor: factor con diversi livelli

`int` and `num`

int: numeri interi: -3, -2, -1, 0, 1, 2, 3

months = c(5, 6, 8, 10, 12, 16)

[1]  5  6  8 10 12 16

num: tutti i numeri da $-\infty$ a $\infty$: 1.1335688, -0.6872158, -0.9923375, 0.5400471, 1.211567, 0.4292457

weight = seq(3, 11, by = 1.5)

[1]  3.0  4.5  6.0  7.5  9.0 10.5

`logi`

Valori loigici, TRUE (T) o FALSE (F)

v_logi = c(TRUE, TRUE, FALSE, FALSE, TRUE)

[1]  TRUE  TRUE FALSE FALSE  TRUE

I vettori logici si possono ottenere attraverso il confronto dei valori in vettori di altro tipo, ad esempio:

months > 12

[1] FALSE FALSE FALSE FALSE FALSE  TRUE

restituisce un vettore logico dove FALSE indica tutti i valori all’interno di months che sono inferiori o uguali a 12, TRUE indica tutti i valori maggiori di 12

`chr` and `factor`

chr: characters: a, b, c, D, E, F

v_chr = c(letters[1:3], LETTERS[4:6])

[1] "a" "b" "c" "D" "E" "F"

factor: usa u numeri o i characters per identificare i livelli della variabile

ses = factor(rep(c("low", "medium", "high"), each = 2))

[1] low    low    medium medium high   high  
Levels: high low medium

Si può cambiare l’ordine dei livelli:

ses1 = factor(ses, levels = c("medium", "high", "low"))

[1] low    low    medium medium high   high  
Levels: medium high low

Creare i vettori

I vettori si possono creare (anche) concatenando insieme variabili diverse

Si possono creare con la funzione c()

Tutti gli elementi all’interno di c() devono essere separati da una virgola

Concantenare gli oggetti con la funzione c(): vec = c(1, 2, 3, 4, 5)

Utilizzare le sequenze:

-5:5 # vettore di numeri da -5 a 5

 [1] -5 -4 -3 -2 -1  0  1  2  3  4  5

seq(-2.5, 2.5, by = 0.5) # sequenza di numeri da -2.5 a 2.5 in step da 0.5

 [1] -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5

Ripetere gli elementi:

# ripete gli elementi da 1a 3 per 4 volte
rep(1:3, 4)

 [1] 1 2 3 1 2 3 1 2 3 1 2 3

rep(c("condA", "condB"), each = 3)

[1] "condA" "condA" "condA" "condB" "condB" "condB"

rep(c("on", "off"), c(3, 2))

[1] "on"  "on"  "on"  "off" "off"

paste0("item", 1:4)

[1] "item1" "item2" "item3" "item4"

Non mischiare tra loro i vettori

int + num $\rightarrow$ num

int/num + logi $\rightarrow$ int/num

int/num + factor $\rightarrow$ int/num

int/num + chr $\rightarrow$ chr

chr + logi $\rightarrow$ chr

Vettori e operazioni

I vettori possono essere sommati/sottratti/divisi e moltiplicati tra loro.

a = c(1:8)
a

[1] 1 2 3 4 5 6 7 8

b = c(4:1)
b

[1] 4 3 2 1

a - b

[1] -3 -1  1  3  1  3  5  7

Se non hanno la stessa lunghezza e un vettore non ha una lunghezza che è un multiplo del’altro vettore si ottiene un messaggio di avvertimento

Se si applica una funzione, questa viene applicata a tutti gli elementi all’interno del vettore

sqrt(a)

[1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427

Allo stesso modo, la stessa operazione può essere applicata a tutti gli oggetti all’interno del vettore

(a - mean(a))^2 # deviazione dalla media al quadrato

[1] 12.25  6.25  2.25  0.25  0.25  2.25  6.25 12.25

Your turn!

Crea un nuovo script e salvalo all’interno del progetto di R
Assegna i seguenti valori a una variabile chiamata my_var (i valori devono essere concatenati tra loro)

23, 24, 25, 27, 28, 29, 30

Calcola la media di my_var usando la funzione mean()
Calcola la media di my_varusando le funzioni sum() e length()
Trova il valore minimo (min()) e massimo (max()) di my_var

See solution

# crea il vettore e assegnalo a my_var
my_var = c(23, 24, 25, 27, 28, 29, 30)
# calcola la media
mean(my_var)

[1] 26.57143

# calcola la media usando altre funzioni
sum(my_var)/length(my_var)

[1] 26.57143

# trova il minimo e il masssimo
min(my_var); max(my_var)

[1] 23

[1] 30

Indicizzare gli oggetti all’interno di un vettore

In questo ambito, indicizzare è vuol dire aver accesso a un particolare elemento all’interno del vettore, con una sua particolare posizione

names = c("Pasquale", "Egidio", "Giulia", "Livio", "Andrea")

Pasquale	Egidio	Giulia	Livio	Andrea
1	2	3	4	5

I vettori si indicizzano usando il numero corrispondente alla posizione dell’oggetto che si vuole estrarre all’interno di parentesi quadre: vector_name[index]

names[1] $\rightarrow$ Pasquale

names[3] $\rightarrow$ Giulia

names[seq(2, 5, by = 2)] $\rightarrow$ Egidio, Livio

Esempi di indicizzazione di vettori

weight

[1]  3.0  4.5  6.0  7.5  9.0 10.5

weight[2] # secondo elemento di weight

[1] 4.5

(weight[6] = 15.2) # sostituisce il sesto elemento di weight

[1] 15.2

weight[seq(1, 6, by = 2)] # elementi 1, 3, 5

[1] 3 6 9

weight[2:6] # elementi dal 2 al 6 (incluso)

[1]  4.5  6.0  7.5  9.0 15.2

weight[-2] # rimuove il secondo elemento

[1]  3.0  6.0  7.5  9.0 15.2

Indicizzare usando la logica

weight

[1]  3.0  4.5  6.0  7.5  9.0 15.2

Quali sono i valori $> 7$ all’interno di weight?

weight > 7

[1] FALSE FALSE FALSE  TRUE  TRUE  TRUE

“Filtra” il vettore con la logica

weight[weight > 7] # solo i valori > di 7

[1]  7.5  9.0 15.2

weight[weight >= 4.5 & weight < 8] # valori compresi tra 4.5 e 8

[1] 4.5 6.0 7.5

Your turn!

Considerando my_var
- Estrarre il terzo elemento
- Estrai tutti gli elementi dispari e assegnali a una nuova variabile my_vector1
- Estrai tutti gli element $> 25$ da my_vector1

See solution

# Terzo elemento di my_var
my_var[3]

[1] 25

# Estrai tutti gli elementi dispari da my_var e assegnali a my_vector1
my_vector1 = my_var[seq(1, length(my_var), by = 2)]
# Estrai tutti gli elementi maggiori di 25 da my_vector1
my_vector1[my_vector1 > 25]

[1] 28 30

Matrici ed arrays

Per creare una matrice:

matrix(data, nrow, ncol, byrow = TRUE)

Crea una matrice $3 \times 4$:

A = matrix(1:12, nrow=3, ncol = 4, byrow = TRUE)

Etichette:

rownames(A) = paste("a", 
                    1:nrow(A), sep ="_") 
colnames(A) = paste("b", 
                    1:ncol(A), sep ="_")
A

    b_1 b_2 b_3 b_4
a_1   1   2   3   4
a_2   5   6   7   8
a_3   9  10  11  12

Trasposizione:

t(A)

    a_1 a_2 a_3
b_1   1   5   9
b_2   2   6  10
b_3   3   7  11
b_4   4   8  12

Le matrici si possono creare attraverso la concatenazione di colonne e righe:

cbind(a1 = 1:4, a2 = 5:8, a3 = 9:12) # concatenazione per colonna

     a1 a2 a3
[1,]  1  5  9
[2,]  2  6 10
[3,]  3  7 11
[4,]  4  8 12

rbind(a1 = 1:4, a2 = 5.8, a3 = 9:12) # concatenazione per riga

   [,1] [,2] [,3] [,4]
a1  1.0  2.0  3.0  4.0
a2  5.8  5.8  5.8  5.8
a3  9.0 10.0 11.0 12.0

Arrays

Per cerare un array:

array(data, c(nrow, ncol, ntab))

my_array = array(1:30, c(2, 5, 3)) # 2 x 5 x 3 array

, , 1

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    5    7    9
[2,]    2    4    6    8   10

, , 2

     [,1] [,2] [,3] [,4] [,5]
[1,]   11   13   15   17   19
[2,]   12   14   16   18   20

, , 3

....

Indicizzare gli elementi in una matrice:

	[,1]	[,2]	[,3]
[1,]	$[1, 1]$	$[1, 2]$	$[1, 3]$
[2,]	$[2, 1]$	$[2, 2]$	$[2, 3]$
[3,]	$[3, 1]$	$[3, 2]$	$[3, 3]$

Stessa logica usata per i vettori ma aggiungendo la seconda dimensione:

matrix_name[row, column]

A[2, 3] # cella della seconda riga nella terza colonna

[1] 7

A[2, ] # seconda riga

b_1 b_2 b_3 b_4 
  5   6   7   8

A[, 3] # terza colonna

a_1 a_2 a_3 
  3   7  11

Indicizzare gli elementi in un array

Oltre alle righe e alle colonne, va considerata anche la terza dinmensione, il tab

array_name[row, col, tab]

Codice

my_array[2, 1, 3] # cella nella seconda riga della prima colonna della terza tab

[1] 22

Codice

my_array[, , 3] # terza tab

     [,1] [,2] [,3] [,4] [,5]
[1,]   21   23   25   27   29
[2,]   22   24   26   28   30

Codice

my_array[1, ,2] # prima riga della seconda tab

[1] 11 13 15 17 19

Your turn!

Crea una matrice $3 \times 3$ con la tabella del 3 fino a 24
Assegna la matrice alla variabile my_mat
Dai un nome alle righe come “row” e alle colonne come “column”
Trasponi my_mat e assegnala alla variabile my_t
Indice da my_t:
- La prima riga
- La seconda colonna
- La cella $[1,2]$

See solution

# Crea una matrice 3x3 con la tabella del 3 fino a 24
my_mat = matrix(seq(0, 24, by = 3), nrow = 3, ncol = 3)

# Assegna la matrice alla variabile `my_mat`
my_mat

     [,1] [,2] [,3]
[1,]    0    9   18
[2,]    3   12   21
[3,]    6   15   24

# Dai un nome alle righe come "row" e alle colonne come "column"
rownames(my_mat) = paste("row", 1:nrow(my_mat), sep = "_")
colnames(my_mat) = paste("column", 1:ncol(my_mat), sep = "_")

# Trasponi `my_mat` e assegna la trasposta alla variabile `my_t`
my_t = t(my_mat)

# Indice da `my_t`: 
  # - La prima riga
my_t[1, ]

row_1 row_2 row_3 
    0     3     6

  # - La seconda colonna
my_t[,2]

column_1 column_2 column_3 
       3       12       21

  # - La cella [1,2]
my_t[1,2]

[1] 3

Liste

Possono contenere diversi tipi di oggetti (e.g., vettori, data frames, altre liste):

my_list = list(w = weight, m = months, s = ses1, a = A)

Gli elementi nelle liste possono essere indicizzati con $ oppure [[]] e il nome o la posizione dell’elemento che si vuole estrarre:

Indicizzare months:

my_list[["m"]] # my_list$m

[1]  5  6  8 10 12 16

Indicizzare weight:

my_list[[1]] # my_list$weight oppure my_list[["w"]]

[1]  3.0  4.5  6.0  7.5  9.0 15.2

Your turn!

Creare una lista con i seguenti elementi
- my_mat
- my_t
- I valori $> 25$ di my_var
- Assegnare la lista all’oggetto my_list

See solution

# - Creare una lista con i seguenti elementi
#   - `my_mat`
#   - `my_t`
#   - I valori $> 25$ di `my_var`
#   - Assegnare la lista all'oggetto `my_list`

my_list = list(my_mat, my_t, my_var[my_var > 25])

Data frames

I data frames sono liste che contengono vettori anche di tipo diverso ma di stessa lunghezza:

id = paste0("sbj", 1:6)
babies = data.frame(id, months, weight)

babies

    id months weight
1 sbj1      5    3.0
2 sbj2      6    4.5
3 sbj3      8    6.0
4 sbj4     10    7.5
5 sbj5     12    9.0
6 sbj6     16   15.2

Working with data frames

Gli elementi nei data frame si possono indicizzare in modo diverso (e simile a quanto visto per le matrici):

babies$months # colonna "months"

[1]  5  6  8 10 12 16

babies$months[2] # secondo elemento della colonna months

[1] 6

babies[, "id"] # colunna id

[1] "sbj1" "sbj2" "sbj3" "sbj4" "sbj5" "sbj6"

babies[2, ] # seconda riga

    id months weight
2 sbj2      6    4.5

Si può usare la logica per indicizzare gli elementi ed estrarre (od escludere) le osservazioni che soddisfano un determinato criterio:

Codice

babies[babies$weight > 7, ] # tutte le osservazioni che mostrano weight > 7

    id months weight
4 sbj4     10    7.5
5 sbj5     12    9.0
6 sbj6     16   15.2

Codice

babies[babies$id %in% c("sbj1", "sbj6"), ] # ossrevazioni dei sbj1 e sbj6

    id months weight
1 sbj1      5    3.0
6 sbj6     16   15.2

Alcune funzioni utili da applicare ai dataframe

dim(babies) # mostra le dimensioni (numero righe e colonne) del data frame

[1] 6 3

names(babies) # nomi delle variabili

[1] "id"     "months" "weight"

str(babies) # mostra le caratteristiche delle variabili nel data frame

'data.frame':   6 obs. of  3 variables:
 $ id    : chr  "sbj1" "sbj2" "sbj3" "sbj4" ...
 $ months: num  5 6 8 10 12 16
 $ weight: num  3 4.5 6 7.5 9 15.2

View(babies) # apre il data viewer

plot(babies) # pariwise plot

summary(babies) # calcola le statistiche descrittive

      id                months         weight      
 Length:6           Min.   : 5.0   Min.   : 3.000  
 Class :character   1st Qu.: 6.5   1st Qu.: 4.875  
 Mode  :character   Median : 9.0   Median : 6.750  
                    Mean   : 9.5   Mean   : 7.533  
                    3rd Qu.:11.5   3rd Qu.: 8.625  
                    Max.   :16.0   Max.   :15.200

Riordinare i valori (e il data frame)

order():

babies[order(babies$weight), ] # riordina le righe del dataframe in ordine

    id months weight
1 sbj1      5    3.0
2 sbj2      6    4.5
3 sbj3      8    6.0
4 sbj4     10    7.5
5 sbj5     12    9.0
6 sbj6     16   15.2

                # crescente di weight

babies[order(babies$weight,    # riordina le righe in ordine decrescente
             decreasing = T), ]

Aggregate

Aggrega e applica una funzione su una variabile di confronto a seconda di una o più variabili di raggruppamento (per esempio è utile per calcolare le medie dei tempi di risposta in diverse condizioni sperimentali):

# Una variabile di confronto e una variabile di raggruppamento
aggregate(y ~ x, data = data, FUN, ...)

# Più variabili di confronto e di raggruppamento
aggregate(cbind(y1, y2) ~ x1 + x2, data = data, FUN, ...)

Un esempio

L’esempio è basato su un dataset interno a R, ToothGrowth, dove si investiga la crescita dei denti delle cavie in base alla somministrazione di vitamina C

ToothGrowth # Crescita dei denti nelle cavie con la somministrazione di Vitamina C

    len supp dose
1   4.2   VC  0.5
2  11.5   VC  0.5
3   7.3   VC  0.5
4   5.8   VC  0.5
5   6.4   VC  0.5
6  10.0   VC  0.5
7  11.2   VC  0.5
8  11.2   VC  0.5
9   5.2   VC  0.5
10  7.0   VC  0.5
11 16.5   VC  1.0
12 16.5   VC  1.0
13 15.2   VC  1.0
14 17.3   VC  1.0
15 22.5   VC  1.0
16 17.3   VC  1.0
17 13.6   VC  1.0
18 14.5   VC  1.0
19 18.8   VC  1.0
20 15.5   VC  1.0
21 23.6   VC  2.0
22 18.5   VC  2.0
23 33.9   VC  2.0
24 25.5   VC  2.0
25 26.4   VC  2.0
26 32.5   VC  2.0
27 26.7   VC  2.0
28 21.5   VC  2.0
29 23.3   VC  2.0
30 29.5   VC  2.0
31 15.2   OJ  0.5
32 21.5   OJ  0.5
33 17.6   OJ  0.5
34  9.7   OJ  0.5
35 14.5   OJ  0.5
36 10.0   OJ  0.5
37  8.2   OJ  0.5
38  9.4   OJ  0.5
39 16.5   OJ  0.5
40  9.7   OJ  0.5
41 19.7   OJ  1.0
42 23.3   OJ  1.0
43 23.6   OJ  1.0
44 26.4   OJ  1.0
45 20.0   OJ  1.0
46 25.2   OJ  1.0
47 25.8   OJ  1.0
48 21.2   OJ  1.0
49 14.5   OJ  1.0
50 27.3   OJ  1.0
51 25.5   OJ  2.0
52 26.4   OJ  2.0
53 22.4   OJ  2.0
54 24.5   OJ  2.0
55 24.8   OJ  2.0
56 30.9   OJ  2.0
57 26.4   OJ  2.0
58 27.3   OJ  2.0
59 29.4   OJ  2.0
60 23.0   OJ  2.0

Per vedere la lunghezza media len (variabile di confronto) a seconda di due variabili di raggruppamento (supp + dose)

aggregate(len ~ supp + dose, data = ToothGrowth, mean)

  supp dose   len
1   OJ  0.5 13.23
2   VC  0.5  7.98
3   OJ  1.0 22.70
4   VC  1.0 16.77
5   OJ  2.0 26.06
6   VC  2.0 26.14

Reshaping: Long to wide

Solitamente i dati sono organizzati o in wide format (i.e., una riga corrisponde a una singola unità statistica) o in long format (i.e., una riga corrisponde a una singola osservazione).

head(Indometh) # Long format

  Subject time conc
1       1 0.25 1.50
2       1 0.50 0.94
3       1 0.75 0.78
4       1 1.00 0.48
5       1 1.25 0.37
6       1 2.00 0.19

Per passare da long format a wide forma (e viceversa) si utilizza la funzione reshape()

Long to wide

# From long to wide
df.w <- reshape(Indometh, v.names="conc", 
                timevar="time",
                idvar="Subject", direction="wide")
df.w

   Subject conc.0.25 conc.0.5 conc.0.75 conc.1 conc.1.25 conc.2 conc.3 conc.4
1        1      1.50     0.94      0.78   0.48      0.37   0.19   0.12   0.11
12       2      2.03     1.63      0.71   0.70      0.64   0.36   0.32   0.20
....

Reshaping: Wide to long

# From wide to long
df.l <- reshape(df.w, varying=list(2:12), 
                v.names="conc", 
                idvar="Subject", direction="long",
times=c(.25, .5, .75, 1, 1.25, 2, 3, 4, 5, 6, 8))

       Subject time conc
1.0.25       1 0.25 1.50
2.0.25       2 0.25 2.03
3.0.25       3 0.25 2.72
....

df.l[order(df.l$Subject), ] # rirodina il data farme per soggetti

       Subject time conc
1.0.25       1 0.25 1.50
1.0.5        1 0.50 0.94
1.0.75       1 0.75 0.78
....

Your turn!

Crea un data frame con 10 osservazioni e le seguenti colonne:
- id: character, ID dei soggetti
- ses: factor, livello socio economico con tre livelli, low, medium, high (3 low, 5 medium, 2 high)
- income: numeric (e.g., runif(10, min = 1110, max = 5430))
Applicare i seguenti filtri:
- Soggetti con ses == high
- Soggetti con income > 2000

See solution

set.seed(999) # imposta un seed in modo che tutti si abbia gli stessi risultati
ids = paste("sbj", 1:10, sep = "_")
ses = factor(c(rep("low", 3), rep("medium", 5), rep("high", 2)), 
             levels = c("low", "medium", "high"))
income = runif(10, min = 1110, max = 5430)
income = income[order(income)]
dati = data.frame(ids, ses, income)

dati[dati$ses == "high", ]

      ids  ses   income
9   sbj_9 high 4508.746
10 sbj_10 high 4793.367

dati[dati$income > 2000, ]

      ids    ses   income
4   sbj_4 medium 2790.788
5   sbj_5 medium 2798.137
6   sbj_6 medium 3628.822
7   sbj_7 medium 3729.851
8   sbj_8 medium 3786.121
9   sbj_9   high 4508.746
10 sbj_10   high 4793.367

Data input and output

Leggere file esterni

File di testo ASCII o spread sheet possono essere imporati con la funzione read.table()

data = read.table("C:/RcouRse/file.txt", header = TRUE)

Argomenti utili:

header: Se la rpima riga è occupata dai nomi delle variabili, va messo TRUE
sep: separatore delle colonne (e.g., virgola, \t)
dec: Se i deciminali sono nel formato 1.2 oppure1,2?

data = read.csv("C:/RcouRse/file.csv", 
                header = TRUE, sep = ";", 
                dec = ",")

In alternativa alla funzione read.csv(), si può utilizzare la funzione read.csv2(), dove il default per il separatore di colonna è “;”. Per vedere se i dati sono stati importati correttamente:

data = read.csv("data/benessere.csv")
head(data) # mostra le prime 6 righe del dataset

  benessere stipendio genere
1         5 1461.0983      m
2         7 1132.3637      f
3         7 1675.9004      m
4         2  328.9587      f
5         6 1370.0146      m
6         5  954.3540      f

I dati sono importati correttamente! Si hanno le variabili nelle rispettive colonne.

Se i dati non vengono importati correttamente, si ha una situazione del genere:

data = read.csv2("data/benessere.csv")
head(data) # mostra le prime 6 righe del dataset

  benessere.stipendio.genere
1       5,1461.09828023079,m
2       7,1132.36368361099,f
3       7,1675.90040479853,m
4       2,328.958701913838,f
5       6,1370.01460952768,m
6       5,954.354030915469,f

Dove tutte le variabili vengono “schiacciate” in un’unica colonna

Da SPSS (file con estensione .sav):

install.packages("foreign")
library(foreign)
data = read.spss("data.sav", to.data.frame = TRUE)

Combinare data frame diversi

Se hanno lo stesso numero di colonne o lo stesso numero di righe:

all_data = rbind(data, data1, data2) # stesso numero di colonne
all_data = cbind(data, data1, data2) # stesso numero di righe

Se hanno numeri di righe/colonne diverse ma condividono almeno una carattertistica contenuta in una colonna (e.g., ID):

all_data = merge(data1, data2, 
                 by = "ID")

Se ci sono IDs differenti nei due dataset, la funzione aggiunge automaticamente delle righe, mettendo NA nelle colonne del dataframe in cui l’osservazione non era presente

all_data contiene tutte le colonne che erano presenti in data1 e data2.

Esportare i dati

Per esportare i dati in formato txt o csv:

write.table(data, # il dataset che si vuole salvare 
            file = "data/mydata.txt", # directory/nome_del_file.estensione
            header = TRUE, # La prima riga contiene i nomi delle variabili?
            sep = "\t",  # separatore di colonna
            ....) # altri argomenti (leggere la documentazione)

Si può esportare direttamente tutto l’ambiente:

save(dat, file = "exp1_data.rda") # Salva solo un oggetto specifico
save(file = "the_earth.rda")      # salva tutto l'ambiente
load("the_earth.rda")             # importa l'mabiente salvato

Basics of Programming

Bisogna essere pronti a sbagliare, TANTO

Coding is ~~hard~~ art

L’ideale è suddividere il problema in sotto-problemi

Da tenere a mente: You’re not alone! $\rightarrow$ stackoverflow (o Google o CHATGPT) sono degli ottimi aiuti

`ifelse()`

Esecuzione condizionata:

ifelse(test, if true, if false)

ifelse(weight > 7, "big boy", "small boy")

[1] "small boy" "small boy" "small boy" "big boy"   "big boy"   "big boy"

Pros:

è davvero facile da usare
Si possono concatenare diversi ifelse()

Cons:

Funziona bene finché la condizione da testare è semplice

`if () {} else {}`

Nel caso si abbiano condizioni più complesse

if (test_1) {
  command_1
} else {
  command_2
}

`if () {} else {}`

Si possono mettere più condizioni all’interno dello stesso ciclo

if (test_1) {
  command_1
} else if (test_2) {
  command_2
} else {
  command_3
}

test_1 (e test_2) devono risultare in TRUE o FALSE, altrimenti non funziona. Se ci sono valori particolari tipo i missing, codifictai come NA, va specificato

if(!is.na(x)) y <- x^2 else stop("x is missing")

Loops

for() e while()

Ripete un comando per un numero specifico di iterazioni

# Don't do this at home
x <- rnorm(10)      
y <- numeric(10)   # crea un contenitore vuoto per mettere i risultati
for(i in seq_along(x)) { # itera tante volte quanti sono gli elementi inx
y[i] <- x[i] - mean(x)
}

In realtà quel codice poteva essere scritto in maniera molto più semplice:

y = x - mean(x)

Sempre meglio evitare i loop

I loop possono essere evitati attraverso l’useo della funzione apply()!

apply(X, margin, fun)

X: l’elemento su cui vi vuole ripetere la funzione

margin: può essere 1 (la funzione viene ripetuta attraverso le colonne, si ottengono i valori marginali di riga) oppure 2 (la funzione viene ripetuta attraverso le righe, si ottengono i valori marginali di colonna)

X <- matrix(rnorm(20),
            nrow = 5, ncol = 4)
apply(X, 2, max) # si ottiene ik valore massimo per ogni colonna (vettore di lunghezza 4)

[1] -0.5660237  1.3254637  0.9013448  0.6430443

apply(X, 1, max) # si ottiene ik valore massimo per ogni riga (vettore di lunghezza 5)

[1] 1.3254637 0.6430443 0.9387494 0.9013448 0.9576504

Alternativa (sconsigliata) utilizzando un ciclo for()

y = NULL
for (i in 1:ncol(X)) {
  y[i] = max(X[, i])
}
y

[1] -0.5660237  1.3254637  0.9013448  0.6430443

Importare i dati

`Import csv`

Attraverso la funzione read.csv()

data = read.csv("data/benessere.csv", # directory/nome_file.csv
                header = TRUE,  # prima riga occupata dai nomi delle variabili
                sep =",", # separatore di colonne  
                dec = ".") # separatore di decimali
head(data)

  benessere stipendio genere
1         5 1461.0983      m
2         7 1132.3637      f
3         7 1675.9004      m
4         2  328.9587      f
5         6 1370.0146      m
6         5  954.3540      f

Calcolo dei punteggi somma

L’assunto è che il data set sia in wide format, ovvero che ogni riga corrisponda a una unità statistica

Somma attraverso le colonne $\rightarrow$ Punteggi somma marginali di riga

rowSums() (rowMeans() se si vuole calcolare la media)

Somma attraverso le righe $\rightarrow$ Punteggi somma marginali di colonna

colSums() (colMeans() se si vuole calcolare la media)

Dati del benessere

# va installato il pacchetto readxl attraverso il codice install.packages("readxl")
# da far correre solo una volta in console
library(readxl)  
benessere = read_xlsx("data/datiBenessere.xlsx")

head(benessere)

# A tibble: 6 × 19
     ID   età genere  frat item1 item2 item3 item4 item5   au1   au2   au3   au4
  <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1     1    16      1     0     1     2     4     3     4     5     4     5     2
2     2    21      1     1     2     2     3     4     3     5     4     5     2
3     3    28      2     4     2     3     5     1     2     4     4     4     4
4     4    15      2     2     3     4     2     2     2     4     5     5     3
5     5    23      1     3     4     5     1     5     2     3     2     3     2
6     6    31      1     2     2     3     2     1     2     5     1     5     1
# ℹ 6 more variables: au5 <dbl>, au6 <dbl>, au7 <dbl>, au8 <dbl>, au9 <dbl>,
#   au10 <dbl>

Se si volessero calcolare i punteggi marginali di riga, si sarebbe tentati di fare:

rowSums(benessere)

  [1]  65  79  88  74  75  83  78  79  85  87  80  88  83  84  92 103  88 101
 [19]  79  83 112 105  98  86  94  94 120 108  98 100 107 111 113 111 112 104
 [37] 100 108 107 126 126 116 116 110 117 117 106 111 122 123 129 146 134 118
 [55] 121 120 121 148 141 129 130 128 132 148 145 152 136 137 138 145 135 136
 [73] 140 140 143 151 148 158 156 143 149 157 165 162 159 170 163 167 163 144
 [91] 149 168 173 176 162 171 174 169 176 172 175 175 172 182 191 186 175 187
[109] 178 190 189 186 184 174 187 189 190 187 199 192 191 194 199 188 206 198
[127] 207 214 208 213 214 205 205 207 227 212 213 208 209 209 216 211 211 216
[145] 208 216 215 219 222 210

Il problema è che ha poco senso fare una cosa del genere! Nei punteggi somma sono inclusi anche le età, gli id, il genere eccetera

Condizioni basate sulle etichette delle variabili

Deve essere presente una regular expression, una sorta di radice comune delle etichette del dataframe:

colnames(benessere)

[1] "ID" "età" "genere" "frat" "item1" "item2"
"item3" "item4"
[9] "item5" "au1" "au2" "au3" "au4" "au5" "au6"
"au7"
[17] "au8" "au9" "au10"

I nomi delle variabili che contengono la stringa item $\rightarrow$ sono gli item del benessere

I nomi delle variabili che contengono la stringa au $\rightarrow$ sono gli item che valutano l’autostima

grep()e grepl(): sono le funzioni che permettono di filtrare i dati a seconda delle regular expression (regex)

grep("regex", vector)

"regex": regular expression che si vuole trovare nel vettore

vector: vettore entro cui si vuole trovare la regex desiderata

La funzione grepl("regex", vector) funziona in modo analogo

Dall’esempio:

# vettore con i nomi delle variabili del dataset
(my_vector = colnames(benessere))

 [1] "ID"     "età"    "genere" "frat"   "item1"  "item2"  "item3"  "item4" 
 [9] "item5"  "au1"    "au2"    "au3"    "au4"    "au5"    "au6"    "au7"   
[17] "au8"    "au9"    "au10"

grep()

# ottiene le posizioni delle colonne che hanno la stringa au
grep("au", my_vector)

 [1] 10 11 12 13 14 15 16 17 18 19

grepl()

# testa quali colonne contengono la stringa au
grepl("au", my_vector)

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

Per calcolare i punteggi somma per i soggetti (riportati sulle righe), va usata la funzione rowSums() condizionata su specifiche colonne:

rowSums(benessere[, grep("item", 
                             colnames(benessere))])

  [1] 14 14 13 13 17 10 16 16  9 13 14 16 17 12  9 20 14 15 13 10 17 16 18 11 14
 [26] 17 20 13 13 18 12 11 14 17 14 17 11 14 11 14 18  9 11 11 16 18 12 12 18 14
 [51] 19 19 13  9 13 10 11 18 17 17 18 12 12 16 14 16 15 17 18 12 11 11 10 11 16
 [76] 19 10  9 12 12 14 14 20 13 12 20 14 20 16 11 10 17 19 26 15 12 20 16 16 14
[101] 12 11 11 16 18 15 10 16 11 18 17  9 11 13 16 19 13 11 15 14 21 22 14  9 16
[126] 13 16 19 13 14 16 12 11 18 18 17 12 12 14 14 18 14 17 15 10 12 15 19 18  9

Per assegnare i punteggi somma a una nuova variabile all’interno del dataframe

benessere$score_ben = rowSums(benessere[, grep("item", colnames(benessere))])

`summary()`

Permette di ottenere delle prime statistiche descrittive di una variabile o di una dataframe

babies <- read.table("data/babies.tab")
summary(babies$peso)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.411   7.809   9.429   9.970  11.268  17.343

data set

summary(babies)

      id               genere               peso           altezza     
 Length:10          Length:10          Min.   : 5.411   Min.   :46.80  
 Class :character   Class :character   1st Qu.: 7.809   1st Qu.:62.90  
 Mode  :character   Mode  :character   Median : 9.429   Median :78.73  
                                       Mean   : 9.970   Mean   :74.88  
                                       3rd Qu.:11.268   3rd Qu.:84.98  
                                       Max.   :17.343   Max.   :99.22

`table()`

Permette di ottenere delle tavole di frequenze per le variabili categoriali

table(babies$genere)


f m 
6 4

Tabelle di frequenza a doppia entrata

benessere = read.csv("data/benessere.csv", header = T)
# viene creata una nuova variabile dicotomica di benessere basso (se sotto la media)
# o alto (se sopra la media)
benessere$new_benessere = with(benessere, 
                               ifelse(benessere > mean(benessere), 
                                      "low", 
                                      "high"))
with(benessere, table(genere, new_benessere))

      new_benessere
genere high low
     f   32  22
     m   28  18

`table()` e percentuali

Nel caso di una singola variabile, basta dividere le frequenze di ogni categoria per il numero totale di osservazioni e moltiplicare per 100

(table(babies$genere)/nrow(babies))*100


 f  m 
60 40

Nel caso di tabelle a doppia entrata, vanno calcolati prima i valori marginali di riga

my_perc = table(benessere$genere, benessere$new_benessere)
# calcola i valori marinali di riga per il genere
(my_perc = cbind(my_perc, rowSums(my_perc)))

  high low   
f   32  22 54
m   28  18 46

(my_perc[,1]/my_perc[,3])*100 # percentuale di maschi e femmine con benessere alto

       f        m 
59.25926 60.86957

(my_perc[,2]/my_perc[,3])*100 # percentuale di maschi e femmine con benessere basso

       f        m 
40.74074 39.13043

`aggregate()`

Aggrega una variabile (o più variabili) di confronto a seconda di una (o più) variabile di raggruppamento, permettendo di applicare delle funzioni sui singoli gruppi. Ad esempio, è utile per ottenere la media dei tempi di risposta in diverse condizioni sperimentali:

# una variabile di confronto y e una variabile di raggruppamento x
aggregate(y ~ x, data = data, FUN, ...)

# più variabili di confronto y1 e y2 e più variabile di confronto x1 e x2
aggregate(cbind(y1, y2) ~ x1 + x2, data = data, FUN, ...)

Un esempio sui punteggi del benessere

benessere = read.csv("data/benessereScores.csv", 
                     header = T, sep =",")
head(benessere)

  ID età genere frat item1 item2 item3 item4 item5 au1 au2 au3 au4 au5 au6 au7
1  1  16      1    0     1     2     4     3     4   5   4   5   2   3   3   1
2  2  21      1    1     2     2     3     4     3   5   4   5   2   4   4   5
3  3  28      2    4     2     3     5     1     2   4   4   4   4   4   4   4
4  4  15      2    2     3     4     2     2     2   4   5   5   3   5   4   2
5  5  23      1    3     4     5     1     5     2   3   2   3   2   4   3   1
6  6  31      1    2     2     3     2     1     2   5   1   5   1   4   5   1
  au8 au9 au10 score_ben score_au
1   4   2    4        14       33
2   5   1    5        14       40
3   4   2    4        13       38
4   4   2    4        13       38
5   4   1    3        17       26
6   5   1    5        10       33

Si può calcolare la media del punteggio di autostima (variabile di confronto) a seconda del genere (variabile di raggruppamento a due livelli):

aggregate(score_au ~ genere, data = benessere, mean)

  genere score_au
1      1 33.08750
2      2 32.62857

Oppure si può calcolare la media sia del punteggio di benessere che di autostima contemporaneamente a seconda del genere

aggregate(cbind(score_ben, score_au) ~ genere, 
          data = benessere, mean)

  genere score_ben score_au
1      1  14.46250 33.08750
2      2  14.41429 32.62857

Your turn!

Utilizzano il data set “benessereScores.csv”

Ricodificare la variabile frat in modo che diventi una variabile dicotomica a due livelli:
- livello siblings se frat == 0
- livello no_siblings se frat != 0 Suggerimento: funzione ifelse()
calcolare la media di score_ben a seconda della variabile siblings
calcolare la media di score_ben e score_au a seconda della variabile sibilings e gender (assegnare il risultato alla variabile mean_dep)
calcolare la deviazione standard di score_ben e score_au a seconda della variabile sibilings e gender (assegnare il risultato alla variabile sd_dep)
Unire (con la funzione merge()) mean_dep e sd_dep e assegnare il risultato alla variabile descr

WARNING!

Quando si utilizza la funzione merge() i nomi delle variabili devono essere diversi

See solution

# ricodifica la variabile frat in modo che dica "no" quando la persona ha 0 fratelli
# e assegnala a una nuova colonna nel dataframe chiamata "siblings"
benessere$siblings = ifelse(benessere$frat == 0, 
                            "no", "yes")
# calcola la media per score_ben e score_au in base a siblings e genere
mean_dep = aggregate(cbind(score_ben, score_au) ~ siblings + genere, 
                     data = benessere, mean)
# rinomina le colonne in modo che abbiano l'etichetta "mean"
colnames(mean_dep)[3:4] = paste("mean", colnames(mean_dep)[3:4], sep = "_")
# calcola la deviazione standard per score_ben e score_au in base a siblings e genere
sd_dep = aggregate(cbind(score_ben, score_au) ~ siblings + genere, 
                     data = benessere, sd)
# rinomina le colonne in modo che abbiano l'etichetta "sd"
colnames(sd_dep)[3:4] = paste("sd", colnames(sd_dep)[3:4], sep =  "_")
# unisci i due dataframe
descr = merge(mean_dep, sd_dep)
descr

  siblings genere mean_score_ben mean_score_au sd_score_ben sd_score_au
1       no      1       13.90909      33.27273     3.250042    4.452734
2       no      2       14.51852      32.11111     3.309315    3.004270
3      yes      1       14.67241      33.01724     3.347625    3.743961
4      yes      2       14.34884      32.95349     3.228472    4.396717

`tidyverse()`

tidyverse() è un pacchetto specifico pensato per il data science che permette un utilizzo più agevole e forse inutitivo delle funzioni di base di R.

install.packages("tidyverse") # installa il pacchetto ( va fatto solo una volta)
library(tidyverse) # carica la libreria e la rende disponibile

Tutto il suo funzionamento si base sull’operatore %>% (Pipe)

Lo si ottiene con la combinazione di tasti shift + ctrl + M

La logica del suo funzionamento è molto semplice:

object %>%   # oggetto su cui si vogliono svolgere le operazioni 
  grouping %>% # eventuali variabili di ragguppamento
  function # la o le funzioni che si vogliono applicare

Statistice descrittive con `tidyverse`

Per ottenere le statistiche descrittive (media e ds) del punteggio di benessere e autostima considerando le variabili di raggruppamento siblings e genere

benessere %>%  # oggetto (dataframe) su cui si vuole lavorare
  group_by(siblings, genere) %>%  # variabili di raggruppamento
  summarise(m_benessere = mean(score_ben),  # funzioni che si vogliono applicare e variabili di confronto
            sd_benessere = sd(score_ben), 
            m_au = mean(score_au), 
            sd_au = sd(score_au))

# A tibble: 4 × 6
# Groups:   siblings [2]
  siblings genere m_benessere sd_benessere  m_au sd_au
  <chr>     <int>       <dbl>        <dbl> <dbl> <dbl>
1 no            1        13.9         3.25  33.3  4.45
2 no            2        14.5         3.31  32.1  3.00
3 yes           1        14.7         3.35  33.0  3.74
4 yes           2        14.3         3.23  33.0  4.40

Your turn!

Calcolare minimo, massimo e mediana di score_au e score_ben con tidyverse
Importare il dataset babies e calcolare media e deviazione standard delle variabili peso e altezza con tidyverse

See solution

benessere %>% 
    summarise(min_ben = min(score_ben),  
            max_ben = max(score_ben), 
            median_ben = median(score_ben),
            min_au = min(score_au), 
            max_au = max(score_au), 
            median_au = median(score_au))

  min_ben max_ben median_ben min_au max_au median_au
1       9      26         14     25     48        32

# importa il dataset babies
babies = read.table("data/babies.tab")
# calcola le statistiche descrittive
babies %>%  
  summarise(mean_weight = mean(peso), 
            sd_weight = sd(peso), 
            mean_height = mean(altezza), 
            sd_height = sd(altezza))

  mean_weight sd_weight mean_height sd_height
1    9.969853  3.274042    74.87698  16.31357

Grafici

Si hanno a disposizione i grafici “tradizionali” per cui non è necessario scaricare ulteriori pacchetti e i grafici prodotti con ggplot, per cui è necessario scaricare il pacchetto ggplot2

Per entrambe le tipologie:

Funzioni di alto livello $\rightarrow$ Quelle che effettivamente producono il grafico
funzioni di basso livello $\rightarrow$ Quelle che lo rendono più carino

Grafici tradizionali

Funzioni di alto livello:

plot()      # scatter plot, specialized plot methods
boxplot()
hist()      # histogram
qqnorm()    # quantile-quantile plot
barplot()
pie()       # pie chart
pairs()     # scatter plot matrix
persp()     # 3d plot
contour()   # contour plot
coplot()    # conditional plot
interaction.plot()

demo(graphics) produce un tour guidato di tutte le funzioanalità grafiche!

Low level functions

points()       # add points
lines()        # add lines
rect()
polygon()
abline()       # add line with intercept a, slope b
arrows()
text()         # add text in plotting region
mtext()        # add text in margins region
axis()         # customize axes
box()          # box around plot
legend()

Il layout dei plot è molto importante ed è costituito da due regioni principali:

La regione del plot (quella che effettivamente contiene il plot)
La regione dei margini (dove sono contenuti i margini e le etichette degli assi)

Ad esempio, per ottenere uno scatterplot:

x <- runif(50, 0, 2) # genera 50 valori casuali da una distribuzione uniforme
y <- runif(50, 0, 2) # genera 50 valori casuali da una distribuzione uniforme
# produce il plot
plot(x, y, # variabili sull'asse x e y, rispettivamente 
     main="Title", # titolo del grafico
     sub="Subtitle", # sottotitolo 
     xlab="x-label", # etichetta asse x
     ylab="y-label") # etichetta asse y
text(0.6, 0.6, "Text at (0.6, 0.6)") # aggiunge del testo al grafico
# aggiunge linee orizzontal e verticali
abline(h=.6, v=.6, lty=2)

`plot()`

plot(x) # una variabile
plot(x, y) # scatter plot
plot(y ~ x) # scatter plot (codice alternativo)

Si uò cambiare il colore degli elementi del grafico anche condizionando su specifiche variabili di raggruppamento:

with(benessere, 
     plot(score_au, 
          col = ifelse(genere == 1, "blue", "pink"), # condiziona il colore tale per cui se è maschio è blu
          pch = ifelse(genere == 1, 3, 16))) # condiziona la forma del puntino 
# inserisce la legenda
legend(x = 100, y = 48, # coordinate specifiche per la legenda 
       c("Male", "Female"), # etichette 
       pch = c(3, 16),  # mappe dei puntini
       col =c("blue", "pink"), # mappe del colore 
       cex = 2) # grandezza del testo

Esempio di scatterplot: `plot(x, y)`

with(benessere, 
     plot(score_au, score_ben))

Esempio di scatterplot: `plot(y ~ x)`

with(benessere, 
     plot(score_ben ~ score_au))

Agli scatterplot si può aggiungere anche la linea che indica il modello di regressione

with(benessere, 
     plot(score_ben ~ score_au))
# aggiunge la linea
abline(lm(score_ben ~ score_au, data = benessere),# specifica il modello di regressione per le coordinate della linea
          col = "blue", 
       lty = 3, # tipo di linea 
       lwd = 3) # spessore della linea

Esempio di scatterplot quando x è una variabile categorica: `plot(y ~ x)`

benessere$genere <- factor(ifelse(benessere$genere == 1, 
                           "maschio", "femmina"))
plot(score_au ~ genere, data = benessere)

equivale ad usare la funzione boxplot(y ~ x)

boxplot(score_au ~ genere, data = benessere)

`hist()`

Serve per produrre l’istogramma delle frequenze per variabili continue

hist(benessere$score_au, # la variabile che si vuole plottare
     breaks = 20) # il numero di bins

Si può plottare la funzione di densità della variabile plottat in istogramma (solo per variabili conitnue)

Codice

hist(benessere$score_au, 
     density=20, # va specificato che si vuole aggiungere l'istogramma 
     breaks=20, 
     prob=TRUE, # va specificato che sull'asse delle y dovranno esserci valori di probabilità 
     col = "darkblue")
curve(dnorm(x, mean=mean(benessere$score_au),  # si utilizza la funzione dnorm() per avere la densità
            sd=sd(benessere$score_au)), 
      col="darkblue", lwd=4, 
      add=TRUE, # fa in modo che la funzione di densità venga aggiunta al grafico
      yaxt="n")

Modificare il layout

Si possono plottare più grafici in un’unica finesra modifcando il layout

par(mfrow=c(nrows, ncolumns)) #  i pannelli vengono riempiti per riga
par(mfcol=c(nrows, ncolumns)) # i pannelli vengono riempiti per colonna

Multi plot: Righe

Codice

par(mfrow=c(1, 2)) # specifica che si vuole una riga e due colonne

hist(benessere$score_au,density=50, breaks=20, prob=TRUE, 
     main = "Score au")
curve(dnorm(x, mean=mean(benessere$score_au), 
            sd=sd(benessere$score_au)), 
      col="springgreen4", lwd=2, add=TRUE, yaxt="n")

hist(benessere$score_ben,density=50, breaks=20, prob=TRUE, 
     main = "Score benessere")
curve(dnorm(x, mean=mean(benessere$score_ben), 
            sd=sd(benessere$score_ben)), 
      col="royalblue", lwd=2, add=TRUE, yaxt="n")

Multi plot: colonne

Codice

par(mfcol= c(2,2))
# primo pannello in alto a sinistra
hist(benessere$score_au,density=50, breaks=20, prob=TRUE, 
     main = "Score au")
curve(dnorm(x, mean=mean(benessere$score_au), 
            sd=sd(benessere$score_au)), 
      col="springgreen4", lwd=2, add=TRUE, yaxt="n")
# secondo pannello in basso a sinistra
hist(benessere$score_ben,density=50, breaks=20, prob=TRUE, 
     main = "Score benessere")
curve(dnorm(x, mean=mean(benessere$score_ben), 
            sd=sd(benessere$score_ben)), 
      col="royalblue", lwd=2, add=TRUE, yaxt="n")
# terzo pannello in alto a destra
with(benessere, 
     plot(score_au, score_ben, frame = FALSE))
# quarto pannello in basso a adestra
with(benessere, 
     plot(score_au, score_ben, frame = FALSE))
abline(lm(score_ben ~ score_au, data = benessere), col = "blue", lwd = 2)

`barplot()`

Serve per ottenere il plot delle frequeunze per le categorie di una variabile categoriale

freq_item1 = table(benessere$item1)
barplot(freq_item1, 
        main = "Item 1 - Frequencies")

`barplot()`: Frequenze relative

è necessario calcolare prima le proporzioni sul totale per ogni categoria di risposta

Codice

perc_item1 = freq_item1/sum(freq_item1)
barplot(perc_item1, 
        ylim = c(0, 1), # setta l'asse delle y ad essere compreso tra 0 e 1 
        main = "Item 1 - Realtive frequencies")

Nel caso si voglia vedere la frequenza relativa in vase a una variabile di raggruppamento, vanno prima calcolate le frequenze relative per ognuno dei gruppi considerati

Codice

# tabella a doppia entrata delle frequeunze di risposta all'item 1 in base al genere
perc_item1_gender = table(benessere$genere, benessere$item1)
# calcola le frequeunze relative sul totale delle risposte per il livello 1
perc_item1_gender[1,] = perc_item1_gender[1,]/table(benessere$item1)
# calcola le frequeunze relative sul totale delle risposte per il livello 2
perc_item1_gender[2,] = perc_item1_gender[2,]/table(benessere$item1)
# produce il plot, automaticamente suddiviso nei due gruppi
barplot(perc_item1_gender, ylim=c(0,1), 
        legend = rownames(perc_item1_gender))
abline(h = .5, lty = 2, col = "red")

Multi plot: Esempio su tutti gli item

Si può ottenere un grafico come quello illustrato in precedenza per tutti gli item, ripetendo il procedimento per oguno degli item attraverso l’applicazione di un ciclo for

Codice

# isola le etichette degli item di interesse
item_ben = benessere[, grep("item", colnames(benessere))]
# predispone il layout per ospitare tutti i grafici
par(mfrow = c(2, round(ncol(item_ben)/2  + 0.2)))
temp = NULL # crea una variabile temporanea dove mettere le tabelle delle frequeunze necessarie per il plot
for (i in 1:ncol(item_ben)) { # il ciclo itera per ognuno degli item
  temp = table(benessere$genere, item_ben[,i]) # calcola le frequeunze per ogni item a sceonda del genere
  for (j in 1:nrow(temp)) {
    temp[j,] = temp[j,]/table(item_ben[,i]) # calcola le frequeunze percentuali
  }
barplot(temp, ylim=c(0,1), legend = rownames(temp), # produce il plot
        main = colnames(item_ben)[i])
abline(h = .5, lty = 2, col = "red")
}

`interaction.plot()`

è possibile anche ottenere un grafico che mostra la relazione tra la variabile $x$ e $y$ dato l’effetto della variabile $z$

interaction.plot(x, z, y)

Riprendendo l’esempio dei bambini, è possibile che la relazione tra peso e altezza vari in base al genere?

Codice

babies = read.table("data/babies.tab")
# plot base relazione peso e altezza
plot(altezza ~ peso, data = babies, 
     col = ifelse(genere == "f", "pink", "blue"), 
     pch = 16, main = "Weight - height", 
     cex = 3)

Codice

babies$cat_weight = with(babies, 
                         ifelse(
                           peso <= quantile(babies$peso)[2], 
                           "light", ifelse(peso > quantile(babies$peso)[2] & 
                                    peso > quantile(babies$peso)[4], "medium", 
                                    "heavy")))
babies$cat_weight = factor(babies$cat_weight,
                           levels = c("light", 
                                      "medium", 
                                      "heavy"))


babies

       id genere      peso  altezza cat_weight
1   baby1      f  7.424646 62.07722      light
2   baby2      m  7.442727 58.18877      light
3   baby3      f  9.512598 84.52737      heavy
4   baby4      f 11.306349 85.13573     medium
....

with(babies, 
     interaction.plot(cat_weight, genere, altezza))

`ggplot2`

ggplot2 (Grammar of Graphics plot, Wickman, 2016) è uno dei pacchetti migliori per ottenere dei grafici d’effetto con (relativamente) poco sforzo

install.packages("ggplot2") # installa il pacchetto (da far correre solo la priuma volta)
library(ggplot2) # carica il pacchetto

La funzione ggplot() è las base del funzionamento del pacchetto. Il suo funzionamento prevede fondamentalmente 3 elementi principali:

ggplot(data, # il dataframe da cui si vogliono prendere le variabili da plottare
       aes(x = x, y = y)) + # aes sta per aestetichs, definisce le variabili da plottare sulla x, y, i colori, le forme eccetera
  geom_graph.type() # definisce il tipo di grafico specifico che si vuole ottenere

(alcuni) Argomenti di aes():

x = variabile asse x
y = variabile asse y
col = variabile per distinguere i colori delle linee/contorni (factor o charcater)
fill = variabile per distinguere il riempimento delle figure (factor o charcater)
shape = variabile per distinguere la forma delle figure (factor o charcater)
size = variabile per distinguere la grandezza delle figure (numeric o integer)
….

Alcuni tipi di grafici:

geom_point() = plotta punti, se vengono date le variabili x e y è uno scatterplot
geom_line() = plotta le linee
geom_boxplot() = per ottenere i boxplot
geom_violin() = per ottenere i violinplot
geom_hist(): istogramma
geom_bar(): barplot
….

Un esempio di scatterplot basato sul dataset rock (è un dataset interno a R)

ggplot(rock, 
       aes(y=peri,
           x=shape, 
           color =shape, # il colore viene definito sulla base della forma delle rocce
           size = peri)) + # la grandezza dei punti è definita dalla lunghezza del perimentro
  geom_point() + # produce lo scatterplot
  theme_bw() + # cambia lo sfondo del grafico
  theme(legend.position = "none") # toglie la legenda

Un altro esempio, dove si crea un data frame dove si hanno due a disposizione i tempi medi di risposta di due gruppi (variabile gruppo, livelli a1 e a2) in tre wave temporali (variabile tempo, tre livelli, b1, b2,b3`):

dat <- read.table(header=TRUE, text="
gruppo tempo rt
a1 b1 825
a1 b2 792
a1 b3 840
a2 b1 997
a2 b2 902
a2 b3 786
")
# le variabili gruppo e tempo vengono forzate ad essere factor
dat[,1:2] = lapply(dat[,1:2], as.factor)

A livello grafico, si può investigare se la relazione tra i tempi di risposta e la wave temporale cambia in base all’appartenenza di gruppo (interaction plot):

Codice

ggplot(dat, 
       aes(x = tempo, 
           y = rt, 
           group = gruppo, 
           color= gruppo, 
           shape = gruppo, 
           linetype = gruppo)) +  
  geom_point(size = 5) + 
  geom_line(size=1)  + theme_classic() +
  ylab("RT") +  scale_linetype_manual("Task", values =c(3,4),
                                labels = c("A1", "A2")) +
  scale_x_discrete(labels = c("B1", "B2", "B3"))

Scatter plot con `ggplot2`

Codice

states = data.frame(state.x77)
ggplot(states,  # raw data
               aes(x = Illiteracy, y = Murder)) + 
  geom_point(size =3, pch=3) + theme_classic() + theme(legend.position = "none")

Per aggiungere la linea del modello di regressione con tanto di banda che indica l’errore standard della stima, si utilizza la funzione geom_smooth() (si possono specificare vari modelli oltre al modello lineare).

Codice

ggplot(states,  # raw data
               aes(x = Illiteracy, y = Murder)) + 
  geom_point(size =3, pch=10) + theme_classic() + 
  geom_smooth(method="lm", color="red", aes(fill="red"))

Multi plot

Per ottenere i multiplot combinando plot di analisi differenti (e.g, grafico delle statistiche descrittive e grafico dei risultati), è necessario installare il pacchetto gridExtra

install.packages("gridExtra") # installa il pacchetto, da usare solo una volta
library(gridExtra) # carica il pacchetto e lo rende usabile

I multi plot si ottengono con la funzione grid.arrange()

Codice

murder_raw = ggplot(states,  # scatterplot
               aes(x = Illiteracy, y = Murder)) + geom_point(size =3, pch=3) + theme_classic() + theme(legend.position = "none")

murder_lm = ggplot(states,  # scatterplot con  modello di regerssione
               aes(x = Illiteracy, y = Murder)) + geom_point(size =3, pch=3) + theme_classic() + 
  geom_smooth(method="lm", color="red", aes(fill="red")) + theme(legend.position = "none")
grid.arrange(murder_raw, # primo plot
             murder_lm, # secondo plot
             nrow=1) # numero di righe

Multi plot II

states = data.frame(state.x77, state.name = state.name,
                    state.region = state.region)

Nel caso in cui si voglia vedere la relazione tra due variabili nei vari livelli definiti da una terza variabile presente sempre nello stesso dataframe, si può usare la funzione facet_wrap() inclusa nel pacchetto ggplot2.

Ad esempio, prendendo come riferimento il dataset disponibile all’interno di R state.x77 si vuole vedere se la relazione tra illiteracy e numero di omicidi varia al variare della macro area geografica considerata (Northeast, South, North Central, West)

Codice

# si crea il dataframe corretto a partire da state.77
states = data.frame(state.x77, state.name = state.name,
                    state.region = state.region)

ggplot(states, 
       aes(x = Illiteracy, y = Murder )) + geom_point() +
  facet_wrap(~state.region) # imposta automatiacamente di suddividere i grafici a seconda della state region

`boxplot()` e `violinplot()`

Per ottenere i boxplot e i violin plot, i dati devono essere riorganizzati in formato long, come nell’esempio qui sotto:

    id condition mean_time
1 sbj1         A  3.860907
2 sbj1         B  2.131674
3 sbj2         A  3.423382
4 sbj2         B  3.919431
5 sbj3         A  4.241790
6 sbj3         B  4.002556

Prendiamo ad esempio i dati del benessere, facendo esclusivamente riferimento ai punteggi di scala dell’autostima e del benessere:

# seleziona solo le variabili di interesse
small = benessere[, c("ID", "score_au", "score_ben")]
# trasforma il dataframe in fromato long
score_long  = reshape(small, 
        idvar = "ID", 
        times =names(small)[-1], 
        timevar = "score", v.names = "value",
        varying = list(names(small)[-1]), 
        direction = "long")
head(score_long)

           ID    score value
1.score_au  1 score_au    33
2.score_au  2 score_au    40
3.score_au  3 score_au    38
4.score_au  4 score_au    38
5.score_au  5 score_au    26
6.score_au  6 score_au    33

Boxplot:

Codice

ggplot(score_long, 
       aes(x = score, y = value,
           fill = score)) + 
  geom_boxplot()

Violinplot:

Codice

ggplot(score_long, 
       aes(x = score, y = value,
           fill = score)) + 
  geom_violin()

score_long = merge(score_long, benessere[, c("ID", "genere")])

ggplot(score_long, 
       aes(x = score, y = value,
           fill = genere)) + geom_violin(trim = FALSE)

Grafici a barre ed istogrammi

geom_hist(): Istograma per variabili continue

geom_bar(): Grafico a barre per variabili categoriali

Argomenti

geom_bar(stat = "count"): conta automaticamente le frequeunze per ogni categoria, per cui non c’è bisogno di specificare una variabile y in aes()
geom_bar(stat = "identity"): associa un valore presente nel dataframe alla variabile y che deve essere specificata in aes()`

`geom_bar(stat = "count")`:

Ad esempio, per ottenere le frequenze di risposta di un item:

ggplot(benessere, 
       aes(x = item5)) + geom_bar(stat = "count") # conta automaticamente le frequeunze

`geom_bar(stat = "identity")`

Per lo stesso item, si vuole plottare la proporzione di risposte epr ogni categoria $\rightarrow$ si vuole specificare il valore delll’asse y

item_5 = data.frame(table(benessere$item5)/nrow(benessere))
ggplot(item_5, 
       aes(x = Var1, 
           y = Freq)) + # viene specificata la variabile da plottare sulla y
  geom_bar(stat = "identity") # va specificato che si ha una variabile y

Un esempio completo

Si possono mettere insieme le funzioni da tidyverse per ottenere più facilmente e velocemente grafici informativi sulle proprie variabili.

Ad esempio mettiamo di volere vedere la proporzione di risposte in ogni categoria di ogni item per la misurazione del benessere:

Codice

# seleziona solo le colonne con gli id dei soggetti e gli item di benessere
new_item = benessere[, c(grep("ID", colnames(benessere)), 
                        grep("item", colnames(benessere)))]
new_item = pivot_longer(new_item, # funzione di tidyverse per girare i dati 
             cols = !ID, 
             names_to = "item", 
             values_to = "value")

Vanno calcolate le proporzioni per ogni categoria per ogni item

Codice

proporzione = new_item %>% 
  group_by(item, value) %>% 
  summarise(prop = n()/nrow(benessere))
proporzione

# A tibble: 26 × 3
# Groups:   item [5]
   item  value  prop
   <chr> <int> <dbl>
 1 item1     1 0.193
 2 item1     2 0.293
 3 item1     3 0.227
 4 item1     4 0.133
 5 item1     5 0.153
 6 item2     1 0.187
 7 item2     2 0.313
 8 item2     3 0.187
 9 item2     4 0.18 
10 item2     5 0.133
# ℹ 16 more rows

e infine si può otttenere il grafico:

ggplot(proporzione, 
       aes(x = value, y = prop)) + 
  geom_bar(stat = "identity") + 
  facet_wrap(~item)

Esportare i grafici

Automaticamente:

pdf("nome_grafico.pdf")
png("nome_grafico.png")          
tiff("nome_grafico.tiff")
jpeg("nome_grafico.jpeg")
bmp("nome_grafico.bmp")

Ma poi bisogna ricordarsi di chiudere il device grafico con:

dev.off()

Ad esempio, se si vuole salvare il violin plot ottenuto dai punteggi di benessere e autostima:

pdf("il_mio_violin.pdf")
ggplot(score_long, 
       aes(x = score, y = value,
           fill = score)) + geom_violin(trim = FALSE) 
dev.off()

Manualmente:

	[,1]	[,2]	[,3]
[1,]	\([1, 1]\)	\([1, 2]\)	\([1, 3]\)
[2,]	\([2, 1]\)	\([2, 2]\)	\([2, 3]\)
[3,]	\([3, 1]\)	\([3, 2]\)	\([3, 3]\)

Introduzione

Installazione

RStudio IDE

Console

Script

Working directories

Basics

CalcolatRice

Assign

Nomi delle variabili

Chiedere aiuto

Ordine, ordine e ancora ordine

R projects for the win

Creare un nuovo R project

Creare un nuovo R project

Selezionare il tipo di progetto

Assegnargli una directory e dargli un nome

Buttare via la spazzatura

Salvare l’ambiente

Quando salvare l’ambiente?

Your turn!

Strutture in R

Funzioni e loro argomenti

concatenate, c()

Vettori

int and num

logi

chr and factor

Creare i vettori

Non mischiare tra loro i vettori

Vettori e operazioni

Your turn!

Indicizzare gli oggetti all’interno di un vettore

Esempi di indicizzazione di vettori

Indicizzare usando la logica

Your turn!

Matrici ed arrays

Arrays

Indicizzare gli elementi in una matrice:

Indicizzare gli elementi in un array

Your turn!

Liste

Your turn!

Data frames

Working with data frames

Riordinare i valori (e il data frame)

Aggregate

Un esempio

Reshaping: Long to wide

Long to wide

Reshaping: Wide to long

Your turn!

Data input and output

Leggere file esterni

Combinare data frame diversi

Esportare i dati

Basics of Programming

ifelse()

if () {} else {}

if () {} else {}

Loops

Sempre meglio evitare i loop

Importare i dati

Import csv

Calcolo dei punteggi somma

Dati del benessere

Condizioni basate sulle etichette delle variabili

summary()

data set

table()

Tabelle di frequenza a doppia entrata

table() e percentuali

aggregate()

Your turn!

tidyverse()

Statistice descrittive con tidyverse

Your turn!

Grafici

Grafici tradizionali

Funzioni di alto livello:

Calcolat`R`ice

Creare un nuovo `R` project

Strutture in `R`

`concatenate`, `c()`

`int` and `num`

`logi`

`chr` and `factor`

`ifelse()`

`if () {} else {}`

`if () {} else {}`

`Import csv`

`summary()`

`table()`

`table()` e percentuali

`aggregate()`

`tidyverse()`

Statistice descrittive con `tidyverse`

`plot()`

Esempio di scatterplot: `plot(x, y)`

Esempio di scatterplot: `plot(y ~ x)`

Esempio di scatterplot quando x è una variabile categorica: `plot(y ~ x)`

`hist()`

`barplot()`

`barplot()`: Frequenze relative

`interaction.plot()`

`ggplot2`

Scatter plot con `ggplot2`

`boxplot()` e `violinplot()`

`geom_bar(stat = "count")`:

`geom_bar(stat = "identity")`