Il linguaggio S in 60 minuti
(forse anche meno!)
Federico M. Stefanini
Dipartimento di Statistica ‘G. Parenti’
Università degli Studi di Firenze
May 9, 2003
F.M.Stefanini - Facoltà di Agraria 2003
1
In questa introduzione è descritto l’uso della finiestra comandi di S, con
qualche riferimento specifico ad SPlus.
Nella finestra comandi è possibile richiamare le righe battute in precedenza usando la tastiera, in particolare il tasto con la freccia rivolta verso
l’alto (o verso il basso). Inoltre, per spostare il cursore da un punto all’altro
della stessa riga bisogna tenere premuto il tasto ”CTRL” e quindi premere
il tasto recante la freccia verso sinistra ( o destra).
1
Uso della finestra comandi
Operazioni elementari: + − ∗
2+3
Potenze, Radice quadrata, logaritmo naturale e base10, esponenziale
23
sqrt(9)
log(123)
log10(123)
exp(log(3.14))
Espressioni complesse
123*(2342/7 + (log(123)/2 -128))
Sommatoria e produttoria
sum(1, 2, 3, 4)
prod(1, 2, 3)
Per ottenere aiuto digitare, per esempio,
help(”log”)
1.1
Spazio di lavoro
Con il lancio del programma si crea uno spazio di lavoro in memoria che
contiene oggetti di vario genere:
• numeri, stringhe, vettori, matrici
• dataframes, liste
• funzioni
F.M.Stefanini - Facoltà di Agraria 2003
2
Ogni oggetto ha un nome, e digitando il nome di un oggetto compare
il suo contenuto (la sua definizione). Per creare un oggetto si assegna ad
un nome il risultato di un’operazione. Ad esempio uno scalare chiamato
primoscalare è creato:
primoscalare <− 234*34/7
Per visionare il contenuto
primoscalare
Si possono creare anche oggetti che contengono una stringa di caratteri
ad esempio
nomcamposper <− ”Landriano D’Adda”
Per visualizzare l’elenco degli oggetti nello spazio lavoro corrente si usa
una funzione. Si noti che le parentesi tonde denotano l’esecuzione di un
oggetto (che è eseguibile), a differenza del nome di una funzione senza paretesi
tonde.
ls()
ls
Per rimuovere da uno spazio di lavoro un oggetto, ad esempio un oggetto
di nome valorealtezza si procede con la funzione
rm(valorealtezza)
Per chiudere una sessione di lavoro si usa la funzione
q()
Per leggere ed eseguire il contenuto di un file ASCII contenente comandi
in linguaggio S
source(”c:/miacartella/giorno1.s”)
1.2
Elenchi ad una o più dimensioni
In un elenco i valori mancanti sono indicati da Na, oppure NaN. Per ottenere
informazioni sui valori mancanti
is.na(nomeoggetto)
1.2.1
Una dimensione
È possibile creare un oggetto di tipo elenco attraverso la funzione c(). Un
elenco di numeri
primoele <− c( 3.14,728, 6.28)
F.M.Stefanini - Facoltà di Agraria 2003
3
Un elenco di stringhe
secoele <− c(”pippo”, ”pluto”, ”paperino”)
Per trasformare un elenco di numeri in elenco di stringhe
risultato <− paste(primoele,sep=””,collapse=”,”)
risultato
Un elenco di elementi di tipo TRUE e FALSE, abbreviati T ed F.
vlogico <− c(T,T,F,T,F)
Un elenco di elementi vero-falso può essere generato da una espressione
test.
ele2 <− seq(-10,10,2)
ele3 <− (ele2 > 3)
ele3
Ci sono diversi modi per costruire un elenco numerico dotato di una
qualche regolarità:
Interi in un intervallo
5:25
-17:11
Numeri reali in un intervallo con incremento assegnato, con la funzione
seq(),
seq(-5,25,0.5)
Elenchi di uno stesso numero ripetuto, con la funzione rep(),
rep(3.14, 5)
Elenchi ripetuti di numeri
rep(c(1,0), c(3,5))
Molte funzioni operano sugli elenchi un elemento per volta
v3 <− c(4,9,16)
sqrt(v3)
log(v3)
Le operazioni elementari binarie sono effettuate elemento ad elemento,
perciò le dimensioni dei due operandi deve essere identica. Se la dimensione
non è la stessa, S cerca di ricondursi a tale condizione usando ciclicamente i
numeri disponibili. SI SCONSIGLIA L’IMPIEGO QUESTA CARATTERISTICA DI S, INIZIALMENTE.
F.M.Stefanini - Facoltà di Agraria 2003
4
v3 <− c(4,9,16)
v4 <− c(1,2,3)
v3 / v4
Si noti la differenza tra le due espressioni
2:(5-1)
2:5-1
Si possono selezionare sottoinsiemi di elementi di una lista. Si impiegano,
a tal fine, le parentesi quadre.
v7 <− -5:5
v7
v7[ c(1,3,7) ]
v7[2:5]
Per selezionare elementi di un elenco si può impiegare un vettore di elementi vero-falso
v8 <− c(1,2,3,4,5)
v8[ c(F,F,T,T,F) ]
Anolagamente è possibile specificare una condizione per la selezione.
v8 <− c(1,2,3,4,5)
v8[ v8 > 3 ]
Per escludere alcuni elementi si usa il segno meno.
v8 <− c(1,2,3,4,5)
v8[ - c(1,3,5) ]
SI NOTI che le due parentesi senza argomento comportano la restituzione
di tutto l’elenco
v7[]
Per conoscere la dimensione di un elenco, si usa length().
v10 <− (-7:436)
length(v10)
Per conoscere lo stato di un oggetto
is.vector(m3)
F.M.Stefanini - Facoltà di Agraria 2003
1.2.2
5
Due dimensioni
Per elenchi a due dimensioni si usa la funzione matrix(). Richiede gli argomenti: elenco numeri, numero righe, numerocolonne
primamatrice <− matrix(0, ncol= 2, nrow=3)
secondamatrice <− matrix(c(1,2,3,4,5,6), ncol= 2, nrow=3)
terzamatrice <− matrix(c(1,2,3), ncol= 2, nrow=3)
quartamatrice <− matrix(c(1,2,3,4,5,6), ncol= 2, nrow=3, byrow=TRUE)
Elenchi a due dimensioni possono derivare da ”postposizione” per riga o
per colonna.
v1 <− c(1,2,3)
v2 <− c(10,20,30)
m1 <− rbind(v1,v2,v1)
m2 <− cbind(v2,v2,v1,v1)
Per selezionare elementi di un elenco a due dimensioni si procede con le
parentesi quadre.
m3 <− matrix(1/(1:25), ncol= 5, nrow=5)
m3
Un elemento
m3[ 3, 2]
Un sottoelenco
m3[ c(2,3), c(1,2)]
Due righe
m3[c(1,3) ,]
Due colonne
m3[ , c(2,3)]
L’argomento entro parentesi quadre può essere generato da espressioni
logiche.
m3[ m3[,1] > 0.25, ]
Per conoscere la dimensione
dim(m3)
nrow(m3)
F.M.Stefanini - Facoltà di Agraria 2003
6
ncol(m3)
Righe e colonne possono avere un nome. Se si intende nominare una sola
dimensione, quella esclusa ha argomento NULL.
m3
dimnames(m3) <− list(c(”nome1”,”nome2”,”nome3”,”nome4”,”nome5”),NULL)
m3
Per interrogare lo stato di un oggetto
is.matrix(m3)
1.2.3
Più di due dimensioni
I principi visti nei paragrafi precedenti si ripetono con la sola differenza che
la funzione generatrice di un elenco multidimensionale è array()
ele222 <− array(1:8, c(2,2,2))
ele222
Per conoscere la dimensione
dim(ele222)
Per interrogare lo stato di un oggetto
is.array(ele222)
1.3
Operazioni su elenchi
Il linguaggio S consente di effettuare operazioni su elenchi ad una o più
dimensioni senza specificare gli indici su ogni elemento.
Le sommatorie e produttorie del vettore vettx si ottengono come
sum(vettx)
prod(vettx)
La somma e il prodotto cumulata con
cumsum(vettx)
cumprod(vettx)
Molte funzioni accettano come argomento un elenco ad una dimensione e
restituiscono il risultato della funzione iterata su ogni elemento dell’elenco.
sqrt(vettx)
F.M.Stefanini - Facoltà di Agraria 2003
7
Il prodotto interno (matriciale) di due elenchi a uno-due dimensioni, A e
B si ottiene con
A %*% B
È possibile applicare a due elenchi una funzione (di S o creata dall’utente)
elemento ad elemento. Ad esempio
outer(c(1,10,1100),c(2,4,8), ”*”)
In generale il nome di una funzione, ad esempio Fratio, non è racchiuso
tra virgolette
outer(c(1,10,1100),c(2,4,8), FRatio)
e la funzione deve richiedere due argomenti.
Il risultato ha dimensione pari alla somma delle dimensioni degli argomenti.
Per iterare l’applicazione di una funzione sugli elementi di un elenco a
due o più dimensioni si usa, ad esempio
apply( miamatrice,2,sum)
per sommare lungo le colonne (seconda delle due dimensioni) gli elementi.
1.4
Liste
Sono una generalizzazione degli elenchi in quanto una lista contiene oggetti
vari, numerici e non, anche di dimensione differente.
li1 <− list(eta=c(21,24,33),peso=c(74,76,80),nome=c(”Remo”,”Romolo”,”?”))
L’accesso ad una componente può avvenire in diversi modi.
li1[[2]]
li1$nome
Per aggiungere altri componenti alla lista
li1$occupazione <− c(1,1,0)
Per informazioni su di una lista
lenght(li1)
names(li1)
Per specificare che gli oggetti devono essere ricercati iniziando da una
lista li1, si procede con
attach(li1)
Per conoscere lo stato di un oggetto
F.M.Stefanini - Facoltà di Agraria 2003
8
is.data.frame(tmp)
Per compiere operazioni su di una lista applicando una funzione ad ogni
componente e ricevere una lista con i risultati
lapply(mialista, sqrt)
Funzioni simili sono sapply() e tapply().
1.5
Dataframes
Nella specifica di modelli statistici in generale sono richieste informazioni
aggiuntive per rendere i numeri da analizzare inequivocabili. Si consideri
l’elenco numerico c(1,2,3). Questo potrebbe riferirsi ad una variabile in scala
per rapporti, oppure ad una variabile su scala nominale. Di per sè, la codifica
numerica non implica la scelta di una specifica scala di misura. Inoltre, in
alcuni casi si desidera specificare se tali numeri sono stati generati da un
processo di campionamento. Si rimanda alla sezione sull’introduzione all’uso
dei modelli per approfondimenti.
Per creare un dataframe
miodf1 <− data.frame(eta=c(21,24,33),pae=as.factor(c(1,2,2)),
nomepaese=c(”USA”,”ITA”,”ITA”))
miodf1
La variabile pae è un fattore (scala nominale)
is.factor(miodf1$pae)
Se si desidera impiegare le etichette di paese
miodf1 <− data.frame(eta=c(21,24,33),pae=c(1,2,2),
nomepaese=factor(c(”USA”,”ITA”,”ITA”)))
I livelli di un fattore (elenco valori non ripetuti) si ottengono
levels(miodf1$nomepaese)
Se pae è random
is.random(miodf1$pae) <− T
Se fosse su scala ordinale con 2 precedente 1
ordered(miodf1$pae) <− c(”2”,”1”)
ordered(miodf1$pae)
Per interrogare lo stato di una componente
is.factor(miodf1$pae)
is.random(miodf1$pae)
F.M.Stefanini - Facoltà di Agraria 2003
9
is.ordered(miodf1$pae)
Per aggiungere componenti si procede come per le liste.
Per creare vettori con valori etichettati si procede con
fatuno <− cut(1:30, breaks=c(5,15,25),labels=c(”Basso”,”Alto”))
fatuno
1.6
Creare una funzione
Per evitare la ripetizione di sequenze di comandi ed implementare nuove
procedure è possibile definire nuove funzioni.
√
Ad esempio, una funzione con argomento a che restituisce a + 10 è
FPrimafun <− function(a) sqrt(a)+10
FPrimafun
Per modificare una funzione
fix(FPrimafun)
Se una funzione è composta da più linee, occorre impiegare le parentesi
graffe per raggruppare le linee eseguibili. Per calcolare la media, potremmo
scrivere
FSeconF <− function(v1)
{
tmp <− sum(v1)
lunghe <− lenght(v1)
tmp/lunghe
}
L’ultima istruzione è quella che compone il valore del risultato restituito
alla funzione chiamante.
Il linguaggio S offre funzioni di alto livello per un gran numero di quantità
usate in statistica. La funzione per calcolare la media aritmetica è mean().
Un modo estremamente produttivo di procedere nella scrittura di funzioni, prevede l’uso di nome <− edit(nome).
Per il debug di funzioni complesse si può ricorrere alla funzione debugger
inspect()
inspect(Miafunzione(23,3,”A”))
in cui Miafunzione(23,3,”A”) è la funzione da esaminare.
Questa funzione rende accessibile un ambiente in cui eseguire le istruzioni
contenute in una funzione passo a passo, e permette anche di esaminare il
F.M.Stefanini - Facoltà di Agraria 2003
10
contenuto delle variabili temporanee create durante la sua esecuzione. Tra i
comandi eseguibili dopo aver invocato inspect() vi sono:
• objects
• eval nomeviariabile
• step
• where
• mark
• enter
• quit
Le principali strutture di controllo sono
• if(condizione)espressione
• ifelse(condizione,expressione1,expressione2)
• for(aux in 1:10){ . . . }
• repeat{ ;if() break;}
• while(condizione) espressione
1.7
Input ed output
Le funzioni e gli oggetti prodotti una sessione di lavoro sono automaticamente
salvate in una directory di default. Per salvare su disco (dischetto) in formato
ASCII (testo) occore in primo luogo fare un elenco degli oggetti di interesse.
elencooggetti <− ls()
Si procede quindi con la rimozione di quelli non interessanti, attraverso
fix(elencooggetti)
Alternativamente si procede alla costruzione di un elenco
elencooggetti <− c(”v1”,”m2”, ”FSeconF”)
Quindi si scarica su file, fornendo il nome di percorso
dump(elencooggetti, fileout=”c:/miacartella/miofile.s”)
Oppure su dischetto
dump(elencooggetti, fileout=”A:/miofile.s”)
F.M.Stefanini - Facoltà di Agraria 2003
11
Per leggere un file di oggetti in S
source(”A:/miofile.s”)
Per leggere un file contenente soli numeri in ASCII
mieidati <− scan(”c:/miacartella/dati2000.txt”)
Il risultato è il vettore mieidati.
Se i dati sono in formato vario, numerico e carattere, è possibile importarli
da ASCII direttamente, inclusa la prima riga contenente i nomi delle variabili
datitot <− read.table(”c:/miacartella/tuttidati2000.txt”, header=T)
Il risultato è un data frame. Per salvare un dataframe
write.table(mydataframe, ”c:/miacartella/outanalisi2000.txt”)
1.8
Finestre grafiche
Per costruire un grafico occorre aprire una finestra grafica. Un modo semplice
è
windows()
Per dividere lo schermo in aree, in modo di avere più grafici in una stessa
pagina, ad esempio 3 per 2
par(mfrow=c(3,2))
alternativamente
par(mcol=c(3,2))
La funzione par() consente di variare moltissimi parametri del grafico (si
veda l’help in linea).
Il grafico di dispersione si ottiene con
vx <− c(2.1, 7.2, 5.5 ,4.5)
vy <− c(12.0, 17.2, 15 ,12.5)
plot(vx,vy)
È possibile scegliere il valore di moltissime opzioni, ad esempio per avere
linee, punti, linee e punti, linee e punti sovrapposti bisogna assegnare a type
il valore, rispettivamente, l, p, b, o.
plot(vx,vy, xlab=”Dose”,ylab=”resa”,xlim=c(0,10),ylim=c(10,25), type=”b”)
Il titolo si ottiene con
title(”Grafico resa dato dose”)
Per salvare su file si usa il menu dopo avere cliccato sulla finestra grafica.
F.M.Stefanini - Facoltà di Agraria 2003
12
Perc chiudere una finestra si usa
dev.off()
Per salvare su file in formato PDF copiando dalla finestra attiva
dev.copy(pdf.graph,”a:/tmp.pdf”)
dev.off()
È possibile aprire e tenere aperte più finestre grafiche in una sessione di
lavoro.
2
Statistica descrittiva
Con il linguaggio S è possibile implementare funzioni che effettuano specificamente i calcoli desiderati. Tuttava vi sono molte funzioni già pronte da
impiegare per ottenere gli indici riassuntivi più comuni.
Per impiegare tali funzioni, si raccomanda di preparare un dataframe che
contenga le informazioni accessorie sulla natura delle variabili incluse (scala
nominale, per rapporti, ecc...).
Sia dato un dataframe contenente due fattori (variabile su scala nominale
e ordinale) ed due variabili su scala per rapporti.
miodf <− data.frame(peso = c(51,87,87,66,66,78),eta=c(14,34,28,34,29,29)
istruzione = as.factor( c(”Bassa”,”Media”,”Alta”,”Bassa”,”Media”,”Bassa”)),
occhi=as.factor(c(”Blu”,”Mar”,”Ver”,”Blu”,”Mar”,”Blu”)) )
ordered(miodf$istruzione ) <− c(”Bassa”,”Media”,”Alta”)
I principali indici sono calcolati, per ogni variabile, con la funzione
summary(miodf)
In particolare, la variabile occhi è presentata in distribuzione di frequenze
assolute, mentre associate alle variabili continue vi sono numerosi indici.
Analiziamo ora gli indici in dettaglio.
Per ordinare le osservazioni ed ottenere la distribuzione di frequenze
sort(miodf$occhi)
table(miodf$occhi)
La distribuzione di frequenze relative si calcola con
table(miodf$occhi) / sum(table(miodf$occhi))
Per valutare se due osservazioni sono uguali, ad esempio la 1 e la 4
miodf$occhi[1] == miodf$occhi[4]
Per contare quanti individui hanno occhi ”Blu”
F.M.Stefanini - Facoltà di Agraria 2003
13
sum(miodf$occhi == ”Blu”)
Se il fattore è su scala ordinale, possiamo cumulare le frequenze ed, eventualmente ottenere quelle cumulate relative
cumsum(table(miodf$istruzione))
cumsum(table(miodf$istruzione))/sum(table(miodf$istruzione))
Possiamo, inoltre, chiedere se un valore precede (segue, è uguale) l’altro,
miodf$istruzione[2] < miodf$istruzione[6]
Le variabili quantitative possono essere riassunte nella stessa maniera di
quelle non quantitative, ed in più è lecito applicare indici che ricorrono ad
operazioni algebriche sui valori della variabile.
Minimo, massimo, media aritmetica, mediana si calcolano, oltre che con
la funzione summary() anche
min(miodf$peso)
max(miodf$peso)
mean(miodf$peso)
median(miodf$peso)
Il minimo ed il massimo si ottengono anche con
range(miodf$peso)
mentre il campo di variazione con
diff(range(miodf$peso))
I quantili si ottengono con la funzione quantile()
quantile(miodf$peso)
con possibilità di specificare un vettore di probabilità
quantile(miodf$peso, prob=seq(0,1,0.1))
La varianza e la deviazione standard (scostamento quadratico medio) si
calcolano con
var(miodf$peso)
sqrt(var(miodf$peso))
È IMPORTANTE RICORDARE CHE IL DENOMINATORE USATO
DALLA FUNZIONE VAR È n − 1, non il totale delle osservazioni, n. Per
avere la varianza non corretta
var(miodf$peso) * (length(miodf$peso)-1)/length(miodf$peso)
A tal fine è possibile creare una funzione
F.M.Stefanini - Facoltà di Agraria 2003
14
varpop <− function(vettore) var(vettore) * (length(vettore)-1)/length(vettore)
da cui l’equivalente
varpop(miodf$peso)
Per calcolare la correlazione lineare tra due variabili
cor(miodf$peso,miodf$eta)
Gli indici appena visti costituiscono solo una piccola parte, seppur importante, dei possibili indici riassuntivi. Il linguaggio S è stato concepito
per permettere facili estensioni mediante funzioni create dall’utente. Ad
esempio, si consideri la differenza semplice media di un elenco di n valori
P
|xi −xj |
∆1 = n,n
i=1,j=1 n2 −n . Raffigurando il calcolo come la differenza tra ogni elemento e tutti i rimanenti, è possibile visualizzare un elenco a due dimensioni
in cui il generico elemento in posizione i, j è costituito dalla differenza tra
i valori xi e xj del medesimo vettore x. Questa operazione può essere effettuata con la funzione outer(). Avendo intenzione di usare ripetutatmente
questo indice, conviene creare una funzione
DifSemMed <− function(vetto){
dimensione <− length(vetto)
differenze <− outer(vetto,vetto,”-”)
valoreassoluto <− abs(differenze)
somma <− sum(c(valoreassoluto))
denominatore <− dimensioneˆ2 - dimensione
risultato <− somma /denominatore
risultato
}
Un’altro esempio di funzione assai utile è quella che fornisce le mode di
una distribuzione di frequenze
FModa <− function(vect){
freq <− table(vect)
flag <− max(freq)
flag2 <− (freq == flag[1])
table(vect)[flag2]/length(vect)
}
2.1
Metodi grafici
Si consideri un dataframe ottenuto per simulazione (si vedano i paragrafi
successivi).
F.M.Stefanini - Facoltà di Agraria 2003
15
miodf2 <− data.frame(forza1 = rnorm(100,12,3), forza2 = rnorm(100,15,3),
durata= rnorm(100,99,5))
In primo luogo occorre aprire una finestra grafica.
windows()
A questo punto è possibile effettuare rappresentazioni grafiche. Per disegnare un boxplot() per la variabile peso
boxplot(miodf2[c(1,2)])
L’istogramma di frequenze relative è richiamato con
hist(miodf2$forza1, nclass=10, plot=T, probability = T)
L’opzione plot = F restituisce la distribuzione di frequenze con raggruppamento in classi. Le classi possono essere specificate dall’utente, ad esempio,
con
hist(miodf2$forza1, breaks=seq(4,25,3), probability = T)
È possibile sovraimporre una stima kernel della densità con
lines(density(miodf2$forza1))
in cui la funzione density() effettuata la stima che può essere regolata
attraverso la scelta di vari parametri.
Essendo interessati alla verifica della normalità dei dati, si procede con il
quantile-quantile plot. Per forza1 si ottiene con
qqnorm(miodf2$forza1)
qqline(miodf2$forza1)
Se la variabile è un fattore del dataframe allora il grafico risulta adeguatamente modificato
hist(miodf$occhi)
Un diagramma a barre si ottiene con
barplot(c(0.7,0.05,0.1,0.1,0.05),names=c(”Assente”,”Basso”,”Medio”,”Alto”,”Altissimo”))
e riquadro e titolo sono disegnati da
box()
title(”Grado di istruzione statistica in Italia”)
Una rappresentazione di uso similare è il diagramma a punti
dotchart(c(0.7,0.05,0.1,0.1,0.05),labels=c(”Assente”,”Basso”,”Medio”,”Alto”,”Altissimo”))
I diagrammi a settori circolari si ottengono da
F.M.Stefanini - Facoltà di Agraria 2003
16
piechart(c(0.7,0.05,0.1,0.05,0.1),names=c(”Assente”,”Basso”,”Medio”,”Alto”,”Altissimo”),
angle=seq(45,135,len=5), density=10, sub=”Grado di istruzione statistica in
Italia”)
La funzione split() risulta assai utile per dividere una variabile quantitativa in componente definite dal valore di una variabile nominale
split(miodf,miodf$occhi)
Si possono ottenereriassunti grafici per più variabili, ad esempio il diagramma di dispersione per coppied id variabili
pairs(miodf2)
Per rappresentare un diagramma di dispersione in cui si possano distinguere punti appartenenti a gruppi diversi si usa la funzione matplot().
Altre rappresentazioni multivariate per dati appartenenti a gruppi differenti includono i grafici a stella, con stars(), e quelli a facce di Chernoff, con
faces().
Sono accessibili anche rappresentazioni tridimensionali. Ad esempio, si
assuma che durata sia funzione di forza1 e di forza2. Si procede con la stima
della superfice
supgra <− interp(miodf2$forza1,miodf2$forza2,miodf2$durata)
quindi alla rappresentazione 3D
persp(supgra)
oppure a linee di livello (contour)
contour(supgra,nlevels=7)
oppure ad immagine
image(supgra)
3
Modelli statistici e simulazioni Monte Carlo
Il linguaggio S offre una ricca collezione di funzioni per adattare modelli di
vario tipo.
È conveniente iniziare sempre lo studio di un modello lineare con la specifica della parametrizzazione.
La funzione
options(contrasts)
fornisce la scelta della parametrizzazione corrente.
Per cambiare parametrizzazione, la funzione
options(constrasts = c(”contr.sum”, ”contr.poly”))
F.M.Stefanini - Facoltà di Agraria 2003
17
specifica che la sommatoria degli effetti è vincolata a zero, scelta comune
nei modelli di Analisi della Varianza per dati gaussiani.
Un’altra possibilità è di scegliere la parametrizzazione a trattamento, cioè
quella in cui un livello ha effetto pari a zero e di fatto sostituisce la media
generale nel ruolo di parametro di riferimento
options(constrasts = c(”contr.treatment”, ”contr.poly”))
Il passo successivo consiste nella costruzione di un dataframe contenente
i dati. In questa fase si specificheranno anche gli attributi delle variabili
incluse nel dataframe, ad esempio quali sono fattori, quali fattori casuali,
quali fattori ordinati (si veda il paragrafo sui dataframe).
I passi successivi sono specifici per il tipo di modello.
3.1
ANOVA
Si consideri uno studio che riguarda gli effetti di tre fertilizzanti sulla crescita
del mais.
Sia definito il dataframe
stud1 <− data.frame(alte=c(12.3,21.6,17.6,19.1,24.3,28.9),
fert=as.factor(c(”A”,”A”,”B”,”B”,”C”,”C”,)))
L’esplorazione grafica dei dati di partenza si può eseguire con
windows()
par(mfrow=c(1,2))
plot.design(stud1)
plot.design(stud1, fun= var)
a destra il valore delle medie, a sinistra delle varianze per livello del
trattamento.
È anche possibile fare un boxplot dei valori entro livello
plot.factor(stud1)
Un oggetto modello ANOVA è creato con
m1 <− aov(y ˜ fert,data=stud1)
che equivale a yi,j = µ + αi + εi,j se la parametrizzazione è sum oppure a
yi,j = µ1 + αi + εi,j con α1 = 0 se la parametrizzazione è treatment.
Per mostrare un riassunto del modello si procede con
summary(m1)
I valori adattati ed i residui si ottengono, rispettivamente, con
fitted.values(m1)
F.M.Stefanini - Facoltà di Agraria 2003
18
resid(m1)
Una prima analisi grafica dei residui può essere fatta con
hist(resid(m1))
qqnorm(resid(m1))
Il grafico dei valori adattati verso i residui è prodotto con
plot(fitted(m1), resid(m1))
La tavola delle medie di livello si ottiene con
model.tables(m1, type=”means”)
Vi sono molte altre funzioni utili, ad esempio interaction.plot() che disegna
graficamente i valori di interazione per disegni a due o più trattamenti (vie).
La funzione fac.design() produce la matrice del disegno di un fattoriale.
Per richiedere un atavola ANOVA Type III per disegni non bilanciati
library(car)
Anova(modello, type= ”III” )
Infine le funzioni p.adjust(p, method=p.adjust.methods, n=length(p)) e
pairwise.t.test() sono utili per i confronti multipli.
3.2
Regressione
Si consideri uno studio che riguarda gli effetti di tre dosi di fertilizzante sulla
crescita del mais.
Sia definito il dataframe
stud2 <− data.frame(alte=c(12.3,21.6,17.6,19.1,24.3,28.9),
dose=c(1.2,1.3,1.1,1.2,1.3,1.4) )
L’esplorazione grafica dei dati di partenza si può eseguire con
summary(stud2)
Il modello di regressione semplice si adatta con
m2 <− lm(alte ∼ dose,stud2)
summary(m2)
L’esame grafico dei residui si effettua con
windows()
par(mfrow=c(2,3))
plot(m2)
F.M.Stefanini - Facoltà di Agraria 2003
19
Un modello di regresione multipla si ottiene espandendo la formula in
lm(). Ammettiamo di avere una variabile relativa alla quantità di acqua ed
una alla temperatura, rispettivamente
stud2$H2O <− c(13.1, 11.3, 11.5, 9.1, 11.8, 10.4 )
stud2$temp <− c(27.2, 29.8, 27.4, 23.1, 28.3, 26.9 )
Il modello è adattato con
m3 <− lm(alte ∼ dose + H2O + temp,stud2)
summary(m3)
ed i residui sono diagrammati con
plot(m3)
La tavola F del modello si ottiene con
anova(m3)
Le funzioni drop1() e add1() consentono di togliere o aggiungere un termine al modello.
La funzione step() effettua la ricerca stepwise.
La funzione
predict(m3, se.fit=T)
calcola i valori attesi ed il reltivo errore standard.
3.3
Simulazione
Il linguaggio S possiede funzioni per generare numeri pseudocasuali estratti
da distribuzioni quali, la normale, l’esponenziale, l’uniforme, la binomiale, la
Poisson, la t, la F, il χ2 , ed altre ancora.
Le relative funzioni iniziano con la lettera r, per indicare random, generatore di numeri casuali, seguito dal nome della distribuzione, ad esempio per
la Gaussiana
rnorm(10, 3.14, 2)
genera un campione di 10 osservazioni dalla Gaussiana (normale) con
media 3.14 e varianza 22 = 4. In altri termini il terzo argomento di rnorm()
è la deviazione standard.
Anologamente si può generare un campioni con n = 100 dall’esponenziale
con tasso 7.1
rexp(100, 7.1)
Tra le altre funzioni si trovano rbinom per la binomiale rchisq per il chiquadro, rf per la F, rgamma per la gamma, rgeom per la geometrica, rhyper
F.M.Stefanini - Facoltà di Agraria 2003
20
per l’ipergeometrica, rlnorm per la log-normale, rpois per la Poisson, rt per
il t, runif per l’uniforme.
È anche possibile specificare un vettore di valori e relativa probabilità per
ogni elemento
sample(c(1,3,7),size= 10,replace=T,prob=c(0.2,0.5,0.3))
che restituisce un vettore di 10 elementi estratti con reimmissione (replace=T)
dalla tripla (1, 3, 7) con probabilità date dal vettore (0.2, 0.5, 0.3).
La funzione sample() è la base per costruire procedure Bootstrap di ricampionamento. Ad esempio, dato un campione (2.1, 4.5, 5.5, 3.2), si ricava
un campione Bootstrap con
sample(c(2.1,4.5,5.5,3.2),size= 4,replace=T)
3.4
Funzioni di distribuzione e quantili
Per calcolare la probabilità dell’evento X ≤ x, con x = 1.16 e X distributita
normalmente con media µ = 1.4 e deviazione standard σ = 5, si procede
mediante
pnorm(1.16,1.4,5)
Per calcolare il quantile x corrispondente ad un valore p = 0.05 dell’evento
X ≤ x, quando X è distributita normalmente con media µ = 1.4 e deviazione
standard σ = 5,
qnorm(0.05,1.4,5)
Il valore di densità di probabilità in x, per X distribuita normalmente
con media µ = 1.4 e deviazione standard σ = 5 si calcola
dnorm(x,1.4,5)
Anolagamente si procede con le altre distribuzioni, in cui i prefissi p,q,d
rimangono inalterati.
4
Risorse web e bibliografia
contiene risorse per la statistica. Si raccomanda la lettura del manuale di R
in PDF. Il programma R si trova a http://www.r-project.org/
Un ipertesto di statistica è il Laboratorio Virtuale di Prob. & Stat.
http://www.ds.unifi.it/VL/
21
F.M.Stefanini - Facoltà di Agraria 2003
Il sito STATLIB
http://lib.stat.cmu.edu
Tra i libri principali:
• Iacus, Masarotto, (2003). Introduzione alla statistica in R. McGrawHill.
• Venables, W.N. and Ripley, B.D. (2003). Modern Applied Statistics
with S-PLUS, Third Edition. Springer-Verlag, New York.
• Selvin, S. (1998). Modern Applied Biostatistical Methods: Using SPLUS. Oxford University Press.
• Krause, A. and Olson, M. (1997).
Springer- Verlag, New York.
The Basics of S and S-PLUS.
• Burns, Patrick (1998). S Poetry. Download for free from:
http://www.seanet.com/ pburns/Spoetry.