readLines("http://macosa.dima.unige.it/R/visitleva.txt",n=11)
[1] "# tabella (dati da visite arruolamento in Marina; 1997; primi contingenti)"
[2] "# (qualche dato relativo al torace e' poco attendibile: la statistica"
[3] "# deve tener conto anche di questo)"
[4] "# torace (troncato ai cm)"
[5] "# peso (troncato ai kg)"
[6] "# regione (numerata in ordine di latitudine del capoluogo: 0ao 1tr 2fr"
[7] "# 3lomb 4ven 5piem 6em-rom 7lig 8tosc 9mar 10um 11abr 12laz 13pugl"
[8] "# 14camp 15sard 16cal 17sic; le altre regioni non sono presenti)"
[9] "# statura (troncata ai cm e diminuita di 100: 178 cm -> 78)"
[10] "torac; peso; regio; statu"
[11] "91;62;15;73"
# Vedo che si sono righe di commento inzianti con # e 1 riga di intestazioni
# I dati sono separati da ";"
# Le righe inizianti con # il programma le salta automaticamente; se ci fossero
# in testa ad es. altre 7 righe da saltare metterei skip=7. C'è anche una riga
# (la 10) di intestazioni: metto header=TRUE in modo che chiami le variabili
# con i nomi torac, peso, regio, statu; altrimenti le chiamerebbe V1,...,V4
dati <- read.table("http://macosa.dima.unige.it/R/visitleva.txt",sep=";",header=TRUE)
str(dati)
'data.frame': 4170 obs. of 4 variables:
$ torac: int 91 86 89 88 92 95 91 83 87 88 ...
$ peso : int 62 64 64 68 70 77 65 62 55 60 ...
$ regio: int 15 15 15 15 15 15 15 15 15 15 ...
$ statu: int 73 70 74 79 66 69 66 79 75 65 ...
# Ho visto che i dati sono molti (non avrebbe avuto senso visiualizzarli con R; e' per
# questo che nel comando readLines ho messo n=11
# Potrei perņ visualizzarli tutti in una tabella esterna, col comando:
# data.entry(dati)
summary(dati)
torac peso regio statu
Min. : 48.00 Min. : 48.00 Min. : 0.00 Min. :47.00
1st Qu.: 86.00 1st Qu.: 63.00 1st Qu.: 7.00 1st Qu.:70.00
Median : 90.00 Median : 70.00 Median :12.00 Median :74.00
Mean : 91.19 Mean : 70.48 Mean :10.18 Mean :74.42
3rd Qu.: 95.00 3rd Qu.: 76.00 3rd Qu.:14.00 3rd Qu.:79.00
Max. :179.00 Max. :126.00 Max. :17.00 Max. :96.00
# Aggiungo 1/2 ai dati troncati agli interi, in modo da non dover poi
# aggiungere 1/2 alle medie
dati$torac <- dati$torac+1/2; dati$peso <- dati$peso+1/2; dati$statu <- dati$statu+1/2
summary(dati)
torac peso regio statu
Min. : 48.50 Min. : 48.50 Min. : 0.00 Min. :47.50
1st Qu.: 86.50 1st Qu.: 63.50 1st Qu.: 7.00 1st Qu.:70.50
Median : 90.50 Median : 70.50 Median :12.00 Median :74.50
Mean : 91.69 Mean : 70.98 Mean :10.18 Mean :74.92
3rd Qu.: 95.50 3rd Qu.: 76.50 3rd Qu.:14.00 3rd Qu.:79.50
Max. :179.50 Max. :126.50 Max. :17.00 Max. :96.50
# i coefficienti di correlazione lineare tra le diverse variabili
cor(dati)
torac peso regio statu
torac 1.0000000 0.69749442 0.01582720 0.2553230
peso 0.6974944 1.00000000 0.02331792 0.4581852
regio 0.0158272 0.02331792 1.00000000 -0.1523683
statu 0.2553230 0.45818521 -0.15236827 1.0000000
# Vediamo, ad es., la relazione tra torace e peso
plot(c(45,180),c(45,130),type="n",xlab="", ylab="")
abline(h=seq(50,130,10),v=seq(40,180,20),lty=3)
# Con densCols posso colorare diversamente le uscite in base alle frequenze
dcols <- densCols(dati$torac,dati$peso)
points(dati$torac,dati$peso,col=dcols,pch=".", cex=5)
# cex dimensiona i punti
# la retta di regressione ("y in funz. di x": rossa)
mod = lm(dati$peso ~ dati$torac); mod$coefficients; abline(mod,col="red")
(Intercept) dati$torac
-15.300034 0.940987
lm(dati$torac ~ dati$peso)$coefficients
(Intercept) dati$peso
54.9919381 0.5170087
# e quella di "x in funzione di y" (verdone)
# x = a y + b -> ay = x-b -> y = (x-b)/a
g <- function(x) (x-54.9919381)/0.5170087
plot(g,45,180,add=TRUE,col="green4")
#
# Elaborazioni simili per la relazione tra "statura-100" e peso:
# Posso, invece che usare diverse intensità di uno stesso colore, usare
# colori diversi. Ecco come:
plot(c(45,180),c(45,130),type="n",xlab="", ylab="")
abline(h=seq(50,130,10),v=seq(40,180,20),lty=3,col="blue")
colori <- c("red","green","blue")
dcols <- densCols(dati$torac,dati$peso,colramp=colorRampPalette(colori))
points(dati$torac,dati$peso,col=dcols,pch=".", cex=4)
#
# Per l'istogramma di "peso" posso usare:
hist(dati$peso,right=FALSE)
# Per avere l'istogramma delle freq. relative posso usare:
hist(dati$peso,right=FALSE,freq=FALSE)
# Sulla scala vert. ho le freq. rel. divise per l'ampiezza dei rettangolini
# I dati sono:
length(dati$peso)
[1] 4170
# La sintesi delle informazioni statistiche:
summary(dati$peso)
Min. 1st Qu. Median Mean 3rd Qu. Max.
48.50 63.50 70.50 70.98 76.50 126.50
#
# Vediamo, infine, un altro strumento spesso usato per studiare le relazioni
# tra due diverse variabili aleatorie relative allo stesso fenomeno:
# i QQ-plot (quantili di una variabile contro i quantili dell'altra).
# Vediamolo nel caso del confronto tra torace e peso:
torace <- dati$torac; peso <- dati$peso
qqplot(peso, torace)
abline(h=axTicks(2),v=axTicks(1),lty=3,col="blue")
# sotto a sinistra la rappresentazione grafica
# sopra a destra è stata sovrapposto un tratteggio
qqplot(peso, torace,col="red")
abline(h=axTicks(2),v=axTicks(1),lty=3,col="blue")
i <- (0:1e4)/1e4; x <- quantile(peso,i); y <- quantile(torace,i)
lines(x,y,lty=3)
# È evidente la tendenza delle due variabili ad essere proporzionali,
# a parte pochi valori estremi.
# Come sappiamo, ciò non accade per peso e altezza:
stat <- dati$statu
qqplot(peso, stat)
abline(h=axTicks(2),v=axTicks(1),lty=3,col="blue")
