########################################################################## # Laboratorio di "Metodi Statistici avanzati" con R # Prof.ssa Lea Petrella # Facolta' di Economia # Dipartimento Metodi e Modelli per l'Economia il Territorio e la Finanza # Sapienza Universita' di Roma # a.a. 2015-2016 # alberto.diiorio@uniroma1.it ########################################################################## ########################################################################## # Lezione 2: INTRODUZIONE ALL'USO DI R (Parte 2) ########################################################################## # # Indice: # # * importazione dei dati # * statistiche descrittive # * rappresentazione grafica dei dati # * fine 2a parte ---> salviamo tutto! # ########################################################################## #----------------------------------------------------------------- # * importazione dei dati #----------------------------------------------------------------- # Leggere dati da 'console' o linea di comando # funzione scan() x=scan() 1 2 3 4 7 32 1 21 45 3 x y=scan() 1 1 1 2 2 2 1 2 1 1 y # nota 1: scan() accetta in imput solo numeri (sarà necessaria la ricodifica per la variabili qualitative) # nota 2: lasciate uno spazio dopo la sequenza di numeri... # nota 3: doppio invio per uscire dati=data.frame(Eta=x,Sesso=factor(y)) #In questo modo stiamo dicendo di leggere la seconda variabile come qualitativa str(dati) #Comando che fa vedere la struttura del dataset levels(dati$Sesso) levels(dati$Sesso)=c("M","F") str(dati) # Importare dati da file esterni # funzione read.table(...) e suoi derivati come # read.csv(...) # read.xls(...) # read.fwf(...) ?read.table # IMPORTANTE: prima di importare un data set aprire il file che lo contiene # per capire come settare gli argomenti della funzione read.table!! # ad esempio vedere se e' presente l'intestazione col nome delle variabili # oppure se ci sono (e quante!) righe iniziali da saltare ecc... # IMPORTANTE: settare la working directory nel percorso dove si trovano i dati getwd() # setwd("C:/Users/alberto/Documents/Tutoraggi/MSA/Lezione2) w=read.table("provastringhe.txt") w w=read.table("provastringhe.txt",stringsAsFactors=F) #Con questa opzione gli stiamo dicendo che le stringhe sono parole. w str(w) gelati=read.table(file="icecream.txt",header=T) str(gelati) # per importare un file .txt possiamo usare il seguente comando dati_school<-read.table("caschool.txt",sep="\t",header=T) dati_school dim(dati_school) str(dati_school) # per importare un file .dat possiamo usare il seguente comando dati_goal<-read.table("goal.dat",header=T) dati_goal dim(dati_goal) str(dati_goal) # uno dei formati piu' "scambiabili" e' il .csv [Comma Separated Values] ?read.csv dati.2 = read.csv("goal.csv") dati.2 dim(dati.2) str(dati.2) #Provare a importare dati indagineUSA.csv dati_USA=read.csv("indagineUSA.csv",head=T,sep=",") #Con R Studio Funzione import dataset #----------------------------------------------------------------- # * Statistiche descrittive #----------------------------------------------------------------- # Variabili qualitative # vediamo qualche rappresentazione per variabili "qualitative" # riprendiamo il data.frame dati_USA str(dati_USA) dati_USA$sesso # Distr. di frequenza assolute table(dati_USA$sesso) sesso.ass=table(dati_USA$sesso) sum(sesso.ass) # Distr. di frequenza relative table(dati_USA$sesso)/length(dati_USA$sesso) sesso.rel=sesso.ass/length(dati_USA$sesso) # verifica sum(sesso.rel) # Distr. di frequenza percentuali table(dati_USA$sesso)/length(dati_USA$sesso)*100 sesso.perc=sesso.rel*100 # verifica sum(sesso.perc) # per avere delle sintesi su tutte le variabili in un colpo solo? summary(dati_USA[,-1]) #Escludo la prima colonna library(epicalc) codebook(dati_USA[,-1]) # Tabelle a doppia entrata Frequenze assolute table(dati_USA$razza,dati_USA$sesso) razza_sesso.tabASS=table(dati_USA$razza,dati_USA$sesso) sum(razza_sesso.tab) # Tabelle a doppia entrata Frequenze relative table(dati_USA$razza,dati_USA$sesso)/sum(table(dati_USA$razza,dati_USA$sesso)) razza_sesso.tabREL=table(dati_USA$razza,dati_USA$sesso)/sum(table(dati_USA$razza,dati_USA$sesso)) sum(razza_sesso.tabREL) # comandi molto efficienti ma SOLO se le variabili sono definite opportunamente # nella struttura del data frame altrimenti si possono ottenere # sintesi inappropriate! # ... e la moda? # posso guardare alle frequenze ma anche determinarla cosi' names(table(dati_USA$sesso))[which.max(table(dati_USA$sesso))] # passiamo a qualche rappresentazione grafica # Barplot barplot(dati_USA$razza) # OPS! barplot(table(dati_USA$razza)) barplot(table(dati_school$gr_span)) # Impariamo ad usare qualche parametro grafico barplot(table(dati_USA$razza), main="Distribuzione per razza",xlab="Razza",ylab="Frequenza") # Barplot per valutare più variabili insieme barplot(table(dati_USA$razza,dati_USA$sesso)) #Rendiamolo più leggibile barplot(table(dati_USA$razza,dati_USA$sesso), main="Distribuzione per razza e sesso",xlab="Sesso",ylab="Frequenza",col=c("red","orange","coral")) #Risulta ancora poco leggibile, abbiamo bisogno di una legenda barplot(razza_sesso.tabASS,col=c("red","orange","coral"),legend=TRUE) #OPS #Muoviamo la legenda barplot(razza_sesso.tabASS,col=c("red","orange","coral"),legend=TRUE,args.legend = list(x = ncol(razza_sesso.tabASS)+0.5 , y=max(colSums(razza_sesso.tabASS))+100, bty = "n")) # Grafico a torta pie(table(dati_USA$razza)) pie(table(dati_school$gr_span),col=c("green","red")) # ...mmmmmhhh aggiungiamo qualcosa... razza.perc=table(dati_USA$razza)/length(dati_USA$razza)*100 pie(table(dati_USA$razza), main="Distribuzione per razza",col=c("green","violet","blue")) legend("topright",legend=paste(razza.perc,"%"),col=c("green","violet","blue"),pch=15,bty="n",cex=1.5,horiz=F) sposato.perc=table(dati_USA$sposato)/length(dati_USA$sposato)*100 par(xpd=TRUE) pie(table(dati_USA$sposato), main="Distribuzione per stato civile",col=1:nlevels(dati_USA$sposato)) legend("topright",legend=paste(sposato.perc,"%"),col=1:nlevels(dati_USA$sposato),pch=15,bty="n",horiz=F) par(xpd=TRUE) pie(table(dati_school$gr_span), main="Distribuzione grade span of district",col=c("green","red")) legend("topright",legend=table(dati_school$gr_span),col=c("green","red"),pch=15,bty="n",horiz=F) # provate con le altre variabili del data set a costruire grafici e sintesi opportune! # Variabili quantitative # Riprendiamo il dataset "dati_goal" str(dati_goal) # chi sono le unita' statistiche? # le colonne fanno riferimento alla medesima variabile (# goal) misurate in anni diversi # lavorando sulla singola variabile... mean(dati_goal$g9697) # per i quantili specifichiamo type=1 visto che si tratta di una v.a. discreta quantile(dati_goal$g9697,type=1) quantile(dati_goal$g9697,probs=c(.1,.9),type=1) quantile(dati_goal$g9697,probs=c(.05,.95),type=1) range(dati_goal$g9697) sd(dati_goal$g9697) # in che unita' di misura e' espressa la s.d.? round(sd(dati_goal$g9697),2) var(dati_goal$g9697) # in che unita' di misura e' espressa la varianza? library(labstatR) sigma2(dati_goal$g9697) #Varianza non corretta. Come la potrei ottenere senza la funzione sigma2? sqrt(sigma2(dati_goal$g9697)) # e la devianza? sigma2(dati_goal$g9697)*length(dati_goal$g9697) # ...oppure... var(dati_goal$g9697)*(length(dati_goal$g9697)-1) summary(dati_goal$g9697) table(dati_goal$g9697) # la moda e' dunque 17 names(table(dati_goal$g9697))[which.max(table(dati_goal$g9697))] # possiamo calcolare le freq. , relative e percentuali cumulate cumsum(table(dati_goal$g9697)) # assolute # arrotondiamo round(cumsum(table(dati_goal$g9697)/length(dati_goal$g9697)),2) # relative arrotondate a 2 cifre decimali round(cumsum(table(dati_goal$g9697)/length(dati_goal$g9697)*100),2) # percentuali # passiamo a qualche rappresentazione grafica # Diagramma ad aste plot(table(dati_goal$g9697)) plot(ecdf(dati_goal$g9697)) plot(ecdf(dati_goal$g9697),main="Funz. di ripartizione empirica per il \n numero di goal nell'anno '96/'97",xlab="Numero di goal",ylab="Freq. relativa cumulata") # come si individua la mediana dal grafico della funzione di ripartizione empirica? abline(h=.5,col="red") # aggiungiamo anche le rette il primo e terzo quartile! abline(h=c(.25,.75),col=c("green","blue")) # lavorando sull'intero data set summary(dati_goal) codebook(dati_goal) # Riprendiamo il dataset "dati_school" str(dati_school) head(dati_school) # soffermiamoci sul reddito medio dati_school$avginc summary(dati_school$avginc) mean(dati_school$avginc) # ora per i quantili possiamo omettere l'opzione type=1 median(dati_school$avginc) quantile(dati_school$avginc) quantile(dati_school$avginc,probs=c(.1,.9)) quantile(dati_school$avginc,probs=c(.05,.95)) # coefficiente di variazione (numero puro) per confrontare la variabilita' # di distribuzioni relative a variabili diverse library(labstatR) cv(dati_school$avginc) cv(dati_school$str) # dunque la variabilita' delle scuole rispetto al reddito medio e' quasi 5 volte # la variabilità rispetto al rapporto numero studenti per insegnante boxplot(dati_school$avginc,main="Boxplot per il reddito medio" ,xlab="Reddito medio") ?hist hist(dati_school$avginc) hist(dati_school$avginc,breaks=c(seq(5,25,5),seq(30,60,10))) hist(dati_school$avginc,breaks=20) hist(dati_school$avginc,xlab="Reddito medio",ylab="Densità di frequenza",xlim=c(0,60),ylim=c(0,200), main="Reddito medio nei distretti delle 420 scuole americane", col="green") hist(dati_school$avginc,xlab="Reddito medio",ylab="Densità di frequenza",xlim=c(0,60),ylim=c(0,.1),prob=T, main="Reddito medio nei distretti delle 420 scuole americane", col="green") plot(ecdf(dati_school$avginc),main="Funzione di ripartizione empirica \n per il reddito medio" ,xlab="Reddito medio") rip.inc=ecdf(dati_school$avginc) rip.inc(median(dati_school$avginc)) # quante scuole si trovano in distretti con reddito medio inferiore o uguale a 20? rip.inc(20) # posso mettere su una finestra piu' grafici? # ...vedi l'argomento mfrow della funzione par() ?par par(mfrow=c(2,1)) hist(dati_school$avginc,xlab="Reddito medio",ylab="Densita' di frequenza",xlim=c(0,60),ylim=c(0,.1),prob=T, main="Reddito medio nei distretti delle 420 scuole americane", col="green") hist(dati_school$comp_stu,xlab="Numero computer per studente",ylab="Densita' di frequenza",prob=T, main="Numero di computer per studente nelle 420 scuole americane", col="lightblue") #----------------------------------------------------------------- # * Rappresentazione grafica dei dati #----------------------------------------------------------------- # Abbiamo visto # per variabili qualitative # 1. Barplot # 2. Grafico a torta # per variabili quantitative # 1. Diagramma ad aste (per dati discreti) # 2. Funzione di ripartizione empirica # 3. Boxplot # 4. Istogramma (per dati continui) # 5. Stima della densita' kernel (la vedremo tra poco) # funzione plot() ?plot vec=c(1,2,3,3,9,20) plot(vec) # specificare una distribuzione (discreta) e disegnarla k = 0:4 p = c(1,2,3,2,1) p=p/sum(p) # normalizzo per avere una distribuzione di probabilità # ---> sum(p) plot(k,p) plot(k, p, type = "h", xlab = "k", ylab = "Probabilita'", ylim = c(0,max(p)) ) # aggiungiamo i pallini points(k,p) # aggiungiamo una retta ?abline abline(0,0.05,lty=2) # aggiungiamo retta verticale abline(v=0.5) # aggiungiamo retta orizzontale abline(h=0.09) # aggiungiamo linea di forma qualsiasi lines(0:4,c(0.1,0.04,0.2,0.1,0),col="red") # tipi di rette e pallini e loro grandezza ??pch plot(1:25,1:25,pch=1:25,main="tutti i simboli", xlab="",ylab="",cex=1:3) # pch: tipo di pallino # cex: grandezza del pallino abline(0,0.1,lty=1,lwd=1,col=1) abline(0,0.2,lty=2,lwd=2,col=2) abline(0,0.3,lty=3,lwd=3,col=3) abline(0,0.4,lty=4,lwd=4,col=4) abline(0,0.5,lty=5,lwd=5,col=5) abline(0,0.6,lty=6,lwd=6,col=6) # lty: tipo di linea # lwd: spessore della linea # col: colore della linea plot(density(dati_school$avginc),main="Stima della densita' kernel per il reddito medio") # come sovrappongo la densità kernel al corrispondente istogramma? hist(dati_school$avginc,xlab="Reddito medio",ylab="Densità di frequenza",xlim=c(0,60),ylim=c(0,.1),prob=T, main="Reddito medio nei distretti delle 420 scuole americane", col="green") lines(density(dati_school$avginc),main="Stima della densita' kernel per il reddito medio",col="red") # funzione curve() ?curve curve(x^2,xlim=c(-5,5),ylim=c(-25,25),col="green") curve(x^3 - 5*x, -4, 4,col="purple",add=T) #--------------------------------------- # * fine 2a lezione ---> salviamo tutto! #--------------------------------------- # Per elencare gli oggetti presenti in memoria ls() # dove salviamo? getwd() # Salvare la lista dei comandi savehistory() loadhistory() # Salvare e recuperare il contenuto del workspace save.image("seconda_parte.RData") load("seconda_parte.RData") # Per uscire : q() # ...e rispondiamo Si/No alla richiesta di salvataggio del workspace