Alcuni esempi artefatti.
Ci poniamo in questa sezione di quantificare quanto sia una correlazione bassa, ovvero una correlazione alta, in termini di un confronto con un campione casuale di numerosità n=50.
N=1000
n=50
X=matrix(rnorm(n*N),n,N)
Y=matrix(rnorm(n*N),n,N)
res=diag(cor(X,Y))
hist(res,30,xlim=c(-1,1))Se N è più grande questo non funziona bene perché cor(X,Y) calcola le mutue covarianze tra tutte le colonne, mentre a noi servono solo i termini sulla diagonale. In alternativa,
N=100000
res=rep(0,N)
for(i in 1:N){
X=rnorm(n)
Y=rnorm(n)
res[i]=cor(X,Y)
}
hist(res,100,xlim=c(-1,1),freq=FALSE)Cerchiamo una banda di confidenza empirica al 95%
Teoricamente potremmo cercare una distribuzione che si adatta e cercare intervalli di confidenza attraverso la distribuzione, che in questo caso è ovviamente Gaussiana.
I quantili teorici sono dati da
Si importa una tabella contenente i dati di peso e altezza di 500 individui.
Esploriamo la tabella, scoprendone innanzitutto la classe,
come sono fatti i primi dati (e le etichette delle colonne),
un breve sommario del suo contenuto.
Valutiamo le correlazioni.
Diamo una rappresentazione grafica. Ricordiamo che la rappresentazione grafica di un data frame è l’insieme dei grafici di dispersione delle colonne a due a due.
Diamo una rappresentazione grafica delle prime due colonne con il diagramma di dispersione,
Valutiamo il modello di regressione lineare semplice.
Possiamo ottenere intervalli di confidenza per i coefficienti (nell’ipotesi di modello Gaussiano!)
Rappresentiamo il diagramma di dispersione includendo la retta di regressione.
Rappresentiamo i residui.
La rappresentaziomne dei residui rispetto alle predizioni è più efficace nel rivelare struttura.
Si standardizza la tabella (solo le colonne numeriche).
oppure con
La standardizzazione non ha variato il grafico di dispersione.
La stessa operazione può essere svolta usando comandi più raffinati di R,
Valutiamo le correlazioni separate per sesso.
PAo<-PA[order(PA$sesso),]
plot(PAo$altezza[1:271],PAo$peso[1:271],col="red",xlim=c(min(PA$altezza),max(PA$altezza)),
ylim=c(min(PA$peso),max(PA$peso)))
points(PAo$altezza[272:500],PAo$peso[272:500],col="blue")
abline(pa.lm)
abline(lm(peso~altezza,data=PAo[1:271,]),col="red")
abline(lm(peso~altezza,data=PAo[272:500,]),col="blue")Spesso le dipendenze lineari sono meno ovvie di quanto uno si aspetti.
Si opera una valutazione della variabilità.
Rappresentazione delle bande di confidenza empiriche (linee blu) e parametriche (linee verdi).
La tabella esaminata presentava una correlazione pari a 0.92 tra altezza e peso. Per capire se tale correlazione sia alta o bassa, confrontiamone il valore con quello di un campione casuale. Ricordiamo che avevamo prodotto un campione di correlazioni con il codice
La versione parametrica permette di stimare il p-value, ovvero la probabilità di avere una correlazione tanto estrema quanto la correlazione tra altezza e peso.
Recuperiamo la tabella su altezza e peso di individui.
Valutiamo la varianza spiegata.
Rappresentiamo graficamente i residui.
Una misura quantitativa di asimmetria è la skewness, \[ \frac1n\sum_{k=1}^n\Bigl(\frac{x_k-\bar x}\sigma\Bigr)^3. \]
Una misura quantitativa di presenza più rilevante di valori estremi è la kurtosi (il valore 3 è il valore di riferimento per una distribuzione normale), \[ \frac1n\sum_{k=1}^n\Bigl(\frac{x_k-\bar x}\sigma\Bigr)^4. \]
Valutiamo graficamente le frequenze e le confrontiamo con la densità di una normale avente le stesse media e deviazione standard campionarie. In questo caso il numero di dati è basso e non ci aspettiamo informazioni troppo rilevanti.
hist(pa.lm.res,10,freq=F)
lines(sort(pa.lm.res),dnorm(sort(pa.lm.res),mean(pa.lm.res),sd(pa.lm.res)),col="red")Una analisi quantitativa di vicinanza tra distribuzioni può essere effettuata attraverso il grafico quantile-quantile qqplot. Nel caso del confronto con i quantili teorici di una normale, c’è il comando abbreviato qqnorm. Il comando qqline disegna la retta di riferimento per quantili uguali.
Infine, esistono test di normalità, quali ad esempio il test di Shapiro-Wilk, un test che ha per ipotesi nulla la normalità dei dati. Nell’output il valore interessante è il p-value, che ci indica (se piccolo) se rigettare l’ipotesi nulla.
Per rendere sensata la previsione, estraiamo le ultime 10 osservazioni dalla tabella peso/altezza, e calibriamo un modello regressivo sulle restanti osservazioni.
Confrontiamo le predizioni con i valori effettivi, valutando l’affidabilità della previsione.
plot(PAts[,2:3],col="blue",pch=16)
pa.pt <- predict(patr.lm,newdata=PAts)
abline(patr.lm,col="gray")
points(cbind(PAts[,2],pa.pt),col="red",pch=16)
# Valutiamo intervalli di confidenza e predizione.
# intervalli non parametrici
cf=coef(patr.lm)
rs=resid(patr.lm)
abline(cf-c(quantile(rs,0.025)),col="darkgreen")
abline(cf+c(quantile(rs,0.975)),col="darkgreen")
# bande di confidenza
pa.ptc <- predict(patr.lm,newdata=PAts,interval="confidence",level=0.95)
lines(cbind(PAts$altezza,pa.ptc[,2]),col="green")
lines(cbind(PAts$altezza,pa.ptc[,3]),col="green")
# bande di predizione
pa.ptp <- predict(patr.lm,newdata=PAts,interval="prediction",level=0.95)
lines(cbind(PAts$altezza,pa.ptp[,2]),col="blue")
lines(cbind(PAts$altezza,pa.ptp[,3]),col="blue")Si considera un esempio artefatto di regressione non-lineare (un esempio più realistico sarà esaminato nel seguito).
Struttura nei residui…
Residui più soddisfacenti.
Carichiamo un set di dati standard nell’installazione di R.
Per maggiori informazioni sui dataset disponibili, si può leggere il risultato del comando library(help=datasets).
Usiamo il dataset EuStockMarkets, che riporta il valore degli indici azionari del mercato tedesco (DAX), svizzero (SMI), francese (CAC) e britannico (FTSE)
Valutiamo un modello lineare per l’andamento dell’indice SMI rispetto al tempo.
I residui presentano una chiara struttura.
Alla luce dei residui ottenuti, analizziamo un modello non-lineare, passando ai logaritmi dei dati.
I residui presentano ancora una chiara struttura (non dimentichiamo che è una serie temporale).
Il modello non-lineare si è rivelato più efficace. Non è una conclusione ovvia. Per capirlo, proviamo a studiare l’andamento di \(R^2\) nei due modelli lineare e non-lineare.
Per capire il ruolo della skewness, generare dei dati casuali artificiosamente simmetrici e asimmetrici, e valutare i corrispondenti valori di skewness.
Per capire il ruolo della kurtosi, generare dei dati casuali che presentino artificiosamente valori estremi, e valutare i corrispondenti valori di kurtosi.
Per valutare il risultato di qqnorm, generare serie di dati casuali con distribuzione normale, e valutare l’aspetto tipico dell’output del grafico quantile-quantile.
Generare un vettore V di numeri casuali di lunghezza 1000. Implementare il calcolo diretto (ovvero senza usare il comando quantile) dei quantili relativi ai dati contenuti in V.
Generare 1000 serie di 10000 numeri casuali con distribuzione \(\chi^2(4)\). Calcolare per ogni serie la frequenza dei valori superiori a 7. Trovare poi un intervallo di confidenza empirico al 95% per la probabilità che una variabile aleatoria \(\chi^2(4)\) sia maggiore di 7.
Generare numeri aleatori con distribuzione esponenziale (e parametro a scelta). Esaminare i dati ottenuti per mezzo dell’analisi dei residui per valutarne l’aderenza ad una distribuzione Gaussiana. Ripetere l’analisi con numeri con distribuzione binomiale, Cauchy, \(\chi^2\), geometrica, lognormale, Poisson, Student e uniforme.
Generare una tabella A con 2 colonne e 10000 righe, in modo che le due colonne siano correlate. Implementare i comandi per tracciare un grafico dei punti con coordinate in A, e per un modello regressivo con fattore di uscita la seconda colonna. Analizzare i residui del modello. Tracciare sul grafico la retta di regressione e due bande di confidenza empiriche al 90%.
Generare una tabella con 2 colonne e 700 righe, in modo tale che la seconda colonna contenga numeri casuali esponenziali di tasso 3, e la prima colonna sia tale che la regressione lineare applicata alla tabella con fattore di uscita la seconda colonna restituisca l’intercetta pari (circa) a -3, la pendenza pari (circa) a 2, e in modo tale che i residui siano non nulli.
Effettuare l’analisi di regressione lineare semplice per la tabella 06pesoaltezza differenziando il campione per indice corporeo.
Analizzare i residui del modello lineare e del modello logaritmico per l’andamento dell’indice SMI.
Scaricare dal sito http://seriestoriche.istat.it la Tavola 2.3 sulla popolazione residente e la Tavola 1.20 sulle emissioni di inquinanti in atmosfera. Comporre una tabella di due colonne contenente nella prima colonna la popolazione dal 1980 al 2010, e nella seconda colonna le corrispondenti emissioni di ammoniaca. Effettuare l’analisi di regressione lineare semplice, comprensiva dell’analisi dei residui. Utilizzare i dati sulla popolazione dal 2011 al 2014 per una previsione delle emissioni inquinanti in atmosfera. Corredare la previsione con bande di confidenza e previsione.