1 Vettori Gaussiani

1.1 Vettori Gaussiani nel piano

Generiamo una nube casuale di punti nel piano con distribuzione Gaussiana.

X=rnorm(5000)
Y=rnorm(5000)
plot(X,Y,pch=".")

Come prima ma rendendo le varianze differenti

X=rnorm(5000,0,6)
Y=rnorm(5000,0,2)
plot(X,Y,pch=".")

Non è cambiato nulla nella figura precedente. Per apprezzare la differenza richiediamo che le unità di misura sui due assi siano le stesse.

plot(X,Y,pch=".",asp=1)

Cambiamo anche le medie

X=rnorm(5000,10,6)
Y=rnorm(5000,-8,1)
plot(X,Y,pch=".",asp=1,xlim=c(-10,10),ylim=c(-10,10))
plot(X,Y,pch=".",asp=1)

Applichiamo una rotazione

a=pi/6
U=matrix(c(cos(a),-sin(a),sin(a),cos(a)),2,2)
U%*%t(U)
t(U)%*%U
P=matrix(c(X,Y),5000,2)
rotP=P%*%U
plot(P,pch=".",col="blue",xlim=c(-10,30),ylim=c(-15,15))
points(rotP,col="red",pch=".")

1.2 Vettori Gaussiani nello spazio

Generiamo una nube di punti con distribuzione Gaussiana nello spazio, usando un campionamento Gaussiano con componenti indipendenti e ugual varianza

Per visualizzare i grafici tridimensionali, useremo la libreria rgl.

library(rgl)
X=rnorm(5000)
Y=rnorm(5000)
Z=rnorm(5000)
open3d()
plot3d(X,Y,Z)
axes3d(c('x','y','z'),pos=c(0,0,0),labels=FALSE,tick=FALSE,col=2,lwd=3)

Esaminiamo il campionamento di dati con distribuzione Gaussiana con componenti indipendenti ma con varianze differenti.

X=rnorm(5000,0,6)
Y=rnorm(5000,0,3)
Z=rnorm(5000,0,1)
open3d()
plot3d(X,Y,Z)
plot3d(X,Y,Z,xlim=c(-20,20),ylim=c(-20,20),zlim=c(-20,20))
axes3d(c('x','y','z'),pos=c(0,0,0),labels=FALSE,tick=FALSE,col=2,lwd=3)

Determiniamo la covarianza

A=matrix(c(X,Y,Z),5000,3)
Q=cov(A)
Q
cor(A)
eigen(Q)
eigen(Q)$values
eigen(Q)$vectors
round(eigen(Q)$vectors,2)

Calcoliamo la radice quadrata della matrice Q.

U=eigen(Q)$vectors
U%*%t(U)
round(U%*%t(U),2)
round(t(U)%*%U,2)
D=diag(eigen(Q)$values)
D
U%*%D%*%t(U) - Q
sqQ=U%*%sqrt(D)%*%t(U)
round(sqQ%*%sqQ - Q,6)

Generiamo una matrice simmetrica semidefinita positiva mediante numeri casuali.

A=matrix(rnorm(9),3,3)
A
Q=A%*%t(A)
eigen(Q)$values

Generiamo dei dati casuali che abbiano la matrice trovata come matrice di covarianza. La corrispondenza è naturalmente solo approssimativa.

X=rnorm(5000)
Y=rnorm(5000)
Z=rnorm(5000)
B=matrix(c(X,Y,Z),5000,3)
G=t(A%*%t(B))
dim(G)
round(Q-cov(G),2)
open3d()
plot3d(G)

Si può ottenere una matrice casuale usando la covarianza di dati casuali. Il risultato però sarà vicino alla matrice identica.

D=matrix(rnorm(15000),5000,3)
Q=cov(D)
A=eigen(Q)$vectors%*%sqrt(diag(eigen(Q)$values))%*%t(eigen(Q)$vectors)
round(A%*%A-Q,2)

2 Analisi delle componenti principali

2.1 Un primo esempio artificioso di analisi

Primo esempio con matrice artificiosamente 2-dimensionale

X=rnorm(50)
Y=rnorm(50)
A=matrix(nrow=50,ncol=5)
A[,1]=X
A[,2]=X+0.1*Y
A[,3]=(X+Y)/2
A[,4]=Y+0.1*X
A[,5]=Y
round(cor(A),2)
plot(data.frame(A))

Le prime due componenti principali catturano praticamente tutta la proporzione di varianza.

pca=princomp(A)
summary(pca)
plot(pca)
biplot(pca)
loadings(pca)

2.2 Analisi su caratteristiche di auto

Consideriamo il dataset mtcars sulle caratteristiche di alcune auto (dataset standard di R). Prima di procedere con l’analisi, rimuoviamo due variabili categoriche).

mt=mtcars[,-(8:9)]
head(mt)
str(mt)

Le variabili rappresentano le seguenti caratteristiche:

  • mpg: miglia per gallone
  • cyl: numero di cilindri
  • disp: cilindrata
  • hp: cavalli
  • drat: quoziente volante/asse
  • wt: peso (in migliaia di libbre)
  • qsec: tempo necessario a percorrere 1/4 di miglio
  • gear: numero delle marce
  • carb: numero di carburatori

La tabella richiede chiaramente standardizzazione.

biplot(princomp(mt),cex=0.3)

Una prima valutazione della bontà dell’analisi.

mt.pca=princomp(scale(mt))
summary(mt.pca)
screeplot(mt.pca)
plot(cumsum(mt.pca$sdev^2)/sum(mt.pca$sdev^2),type="b",ylim=c(0,1))
segments(1,0.8,9,0.8,col="red")

Proviamo ad associare i fattori alle nuove variabili. Dal grafico si potrebbe ipotizzare:

  • 1ma comp.: mpg, cyl, disp, wt
  • 2da comp.: drat, qsec, gear, carb
  • incerti: hp
biplot(mt.pca,col=c("gray","red"))

Approfondiamo l’associazione per mezzo dei loadings.

  • 1ma comp.: mpg, cyl, disp, wt
  • 2da comp.: qsec, gear, carb
  • incerti: hp, drat
mt.ld=loadings(mt.pca)
mt.ld

Integriamo con il contributo delle rotazioni.

  • 1ma comp.: mpg, cyl, disp, drat, wt
  • 2da comp.: hp, qsec, gear, carb
varimax(mt.ld[,1:2])

Rotazioni successive in questo caso non aiutano.

varimax(mt.ld[,1:3])

Concludiamo con una possibile interpretazione:

  • 1ma comp. (mpg, cyl, disp, hp, wt): caratteristiche legate alla potenza, o al consumo.
  • 2da comp. (drat, qsec, gear, carb): caratteristiche dotazionali (e prestazionali).

Alla luce di questa interpretazione, osserviamo dal biplot che mpg è concorde alla prima coomponente, mentre cyl, disp, hp e wt sono discordi, in accordo con l’idea che l’asse rappresenti un indice di consumo (a sinistra le auto meno prestazionali dal punto di vista del consumo, a destra le auto più prestazioni dal punto di vista del consumo).

Analogamente, tutte le frecce dei fattori associati alla seconda componente sono concordi alla seconda componente (inclusa, in effetti, qsec, in quanto minore è il valore di questo fattore, maggiore la prestazione). Il fattore carb, al contrario degli altri fattori associati alla seconda componente, ha una proiezione negativa sulla prima componente, in accordo con l’interpretazione data. Infine, carb e qsec sono opposti.

Scopriamo che la rappresentazione sul piano delle componenti principali può rivelare aspetti inaspettati. Indichiamo con colori differenti auto di provenienza differente.

orig=rep(2,32)
orig[c(1:3,19:21)]=3
orig[c(4:7,15:17,22:25,29)]=1
plot(mt.pca$scores,pch=20,col=orig)
text(mt.pca$scores,labels=as.character(row.names(mt)),col=orig,pos=3,cex=.6)

Chiediamoci se questo aspetto poteva emergere senza l’analisi delle componenti principali.

plot(mt,col=orig)

Proviamo a indagare sulla stabilità del risultato.

mt_s=data.frame(scale(mt))
res=rep(0,32)
for(i in 1:32){
  mt_r=mt_s[-i,]
  mt_r.pca=princomp(mt_r)
  mt_p=predict(mt_r.pca,newdata=mt_s[i,])[1:2]
  res[i]=mean((mt_p-predict(mt.pca)[i,1:2])^2)
}
sqrt(mean(res))
hist(res)
mt[31,]

2.3 Analisi sulla produzione agricola

Esempio con i dati sulla produzione agricola delle regioni.

pa<-read.csv("08prodagri.csv",row.names=1)
head(pa)

Svolgiamo l’analisi delle componenti principali.

pa.pca=princomp(pa)
summary(pa.pca)

Il grafico delle varianze e delle varianze cumulate:

plot(pa.pca)
plot(cumsum(pa.pca$sdev^2)/sum(pa.pca$sdev^2),type="b",ylim=c(0,1))
segments(1,0.8,7,0.8,col="red")

Rappresentiamo il piano principale.

biplot(pa.pca)

Osserviamo la presenza di frecce più corte e più lunghe, standardizziamo i dati e verifichiamo la validità dell’analisi.

pa_s.pca=princomp(scale(pa))
summary(pa_s.pca)
plot(cumsum(pa_s.pca$sdev^2)/sum(pa_s.pca$sdev^2),type="b",ylim=c(0,1))
segments(1,0.8,7,0.8,col="red")

Confrontiamo le analisi con/senza standardizzazione

layout(t(1:2))
biplot(pa.pca)
biplot(pa_s.pca)
layout(1)

Possiamo associare alla prima componente la produttività, la seconda permette di distinguere coltivazioni che hanno necessità di climi più caldi/freddi.

biplot(pa_s.pca)

Verifichiamo l’interpretazione con le informazioni quantitative. In effetti tutti i fattori contribuiscono in diverse misure alle componenti principali.

pa_s.ld=loadings(pa_s.pca)
pa_s.ld
varimax(pa_s.ld[,1:2])

Usiamo l’asse principale per classificare le regioni rispetto alla produzione agricola.

pa_s.pt=predict(pa_s.pca)
sort(-pa_s.pt[,1])

Esplorazione dei piani. 

layout(t(1:2))
biplot(pa_s.pca,choices=c(1,3))
biplot(pa_s.pca,choices=c(2,3))
layout(1)

Avendo in mente il significato delle prime componenti principali, si può usare l’analisi delle componenti principali in presenza di nuovi dati per capire che produttività e tipo di coltivazioni hanno i nuovi individui. Per esempio:

pa_s=data.frame(scale(pa))
pa_sm=pa_s[-4,]
pa_sm.pca=princomp(pa_sm)
predict(pa_sm.pca,newdata=pa_s[4,])[c(1,2)]
predict(pa_s.pca)[4,c(1,2)]

2.4 Analisi del dataset relativo agli Iris

Un altro esempio di analisi delle componenti principali usando uno dei dataset standard di R. I dati si riferiscono alle misure (in centimetri) di lunghezza e larghezza dei sepali e dei petali di 50 fiori da tre diverse specie di iris.

library(datasets)
head(iris)
summary(iris)

Per ora dimentichiamo la variabile categoriale. La salviamo però in un vettore per uso successivo.

ir=iris[,1:4]
species=factor(t(iris[5]))

Valutiamo numericamente e graficamente le correlazioni.

round(cor(ir),2)
plot(ir,pch=19)

o meglio ancora, per evidenziare le differenze tra le specie,

plot(ir,pch=19, col=as.numeric(species)+1)

Analisi delle componenti principali.

ir.pca=princomp(ir)
summary(ir.pca)

Le proporzioni di varianza spiegata.

plot(ir.pca)
plot(cumsum(ir.pca$sdev^2)/sum(ir.pca$sdev^2),type="b",ylim=c(0,1))

Il biplot mostra frecce di lunghezza disomogenea (l’opzione cex è il character expansion factor, qui usato per rendere i numeri più piccoli e poter osservare le frecce). Risulta opportuno (come sempre nell’analisi PCA) standardizzare.

biplot(ir.pca,cex=.4)

L’analisi delle componenti principali continua ad essere buona.

ir_s.pca=princomp(scale(ir))
summary(ir_s.pca)
plot(cumsum(ir_s.pca$sdev^2)/sum(ir_s.pca$sdev^2),type="b",ylim=c(0,1))

Notiamo il miglioramento con la standardizzazione.

layout(t(1:2))
biplot(ir.pca,main="iris",cex=0.6)
biplot(ir_s.pca,main="iris std.",cex=0.6)
layout(1)

Interpretazione delle componenti: la prima componente sembra venire dal contributo (approssimativamente) simile di tre dei quattro fattori, con il ruolo di separare il quarto fattore dagli altri. La seconda componente principale ha tutti i valori negativi e si potrebbe interpretare come una misura della dimensione (principalmente dei sepali).

biplot(ir_s.pca)
loadings(ir_s.pca)
varimax(loadings(ir_s.pca)[,1:2])

Studiamo gli assembramenti, utilizzando la variabile categoriale.

levels(species)<-c("red","blue","green3")
ir_s.pt <- predict(ir_s.pca)
layout(t(1:2))
plot(ir_s.pt[,c(1,2)],pch=20,asp=1,col=species)
plot(ir_s.pt[,c(1,3)],pch=20,asp=1,col=species)
layout(1)

Eventualmente anche…

plot(ir_s.pt[,c(2,3)],pch=20,asp=1,col=species)

2.5 Analisi dei dati socio-economici di alcuni paesi europei.

Consideriamo i dati contenuti in un dataset relativo a dati socio-economici di alcuni paesi europei.

eu<-read.csv("10euecosoc.csv",header=TRUE,row.names=1)
head(eu)
str(eu)
  • Education la percentuale della popolazione che risulta diplomata o laureata;
  • Health la percentuale della popolazione che ritiene di essere in buono stato di salute;
  • Accidents il numero di incidenti sul lavoro che hanno richiesto almeno 4 giorni di assenza dal posto di lavoro, calcolati sulla percentuale di popolazione lavoratrice;
  • Employment lo stipendio annuale medio;
  • Export le quote di mercato sulle esportazioni;
  • HealthyYears la speranza di vita media (media tra uomini e donne);
  • GenderPaymentGap una sorta di indicatore della diversità di retribuzione tra uomini e donne;
  • IndustrialProduction un indicatore della produzione industriale del Paese;
  • TaxesLabour e TaxesTotal, due indicatori riguardanti la pressione fiscale.

Una prima analisi delle correlazioni.

library(corrplot)
corrplot(cor(eu),"number")
plot(eu)

L’analisi delle componenti principali.

eu.pca=princomp(eu)
summary(eu.pca)

Valutazione della bontà dell’analisi. Per arrivare a circa l’80% è necessario considerare tre componenti principali.

plot(eu.pca)
plot(cumsum(eu.pca$sdev^2)/sum(eu.pca$sdev^2),type="b",ylim=c(0,1))

Ancora una volta il biplot rivela la necessità della standardizzazione.

biplot(eu.pca,cex=.2)

Ripetiamo l’analisi delle componenti principali.

eu.pca=princomp(scale(eu))
summary(eu.pca)

Valutazione della bontà dell’analisi. Il problema si rivela più complesso dimensionalmente.

plot(eu.pca)
plot(cumsum(eu.pca$sdev^2)/sum(eu.pca$sdev^2),type="b",ylim=c(0,1))
segments(1,0.8,9,0.8,col="red")

L’attribuzione di significato risulta più complessa.

biplot(eu.pca,cex=.4)

Anche guardando alle informazioni quantitative.

eu.ld=loadings(eu.pca)
eu.ld

L’assegnazione di Health, HealthyYears, IndustrialProduction, TaxesTotal alla prima componente sembra essere abbastanza stabile.

Export e TaxesLabour sembrano essere associati alla seconda componente.

varimax(eu.ld[,1:2])
varimax(eu.ld[,1:3])
varimax(eu.ld[,1:4])

Si può tentare di esplorare anche i piani principali successivi.

biplot(eu.pca,choice=c(1,3),cex=.6)
biplot(eu.pca,choice=c(2,3),cex=.6)

3 Esercizi

Generare una matrice A simmetrica e semidefinita positiva con 3 righe e 3 colonne. Implementare il calcolo della radice quarta di A, ovvero di una matrice Q tale che Q^4=A.

Generare 100 coppie di dati Gaussiani che abbiano la bisettrice del primo quadrante come asse principale e varianza 9 lungo tale direzione. Dire quale deve essere l’altro asse principale.

Si generi una tabella di 3 colonne e 25 righe. Le prime due colonne siano Gaussiane (media e varianza a piacere), la terza colonna sia la somma delle prime due colonne più un piccolo rumore. Ottenere le varianze spiegate delle componenti principali della tabella, rappresentarle in un grafico illustrativo e visualizzare il grafico delle prime due componenti principali.

Generare un vettore Gaussiano di 5000 punti con componenti principali di varianze 9, 4 e 1 ed in modo che i suoi assi principali siano le bisettrici dei quadranti del piano xy e l’asse delle quote. Ottenere il grafico della proiezione dei punti sul piano xz.

Generare un campione Gaussiano quadri-dimensionale di numerosità 100, in modo che questo risulti approssimativamente bidimensionale e che la prima componente principale sia la somma delle prime due componenti.

Creare una tabella con 5 colonne e 350 righe, e popolarla di valori in modo che il biplot dell’analisi delle componenti principali per la tabella creata presenti alcune frecce molto corte.

Generare una tabella di dati in modo che, analizzata attraverso l’analisi delle componenti principali, riveli un accumulo di dati sul piano principale che però si rivela ingannevole se visto attraverso gli altri piani.

Considerare nuovamente i dati (non standardizzati) relativi alla produzione agricola delle regioni.

  1. Normalizzare i dati in termini di distribuzione pro-capite usando i dati sulle popolazioni contenuti nella tabella 09prodagri.csv.
  2. Ricavare l’analisi delle componenti principali del nuovo set di dati e osservare che è necessario standardizzare la tabella.
  3. Ricavare l’analisi delle componenti principali del set di dati della produzione agricola pro-capite standardizzata, eventualmente analizzando la posizione delle regioni nei primi piani principali

Considerare la tabella di dati iris. Studiare la corrispondente tabella ottenuta attraverso i logaritmi dei dati originari (ovvero, indagare rispetto a una dipendenza di tipo potenza tra i fattori).

  1. Valutare l’opportunità di standardizzare la tabella ottenuta.
  2. Ricavare l’analisi delle componenti principali per la nuova tabella di dati.
  3. Verificare, colorando i punti rappresentati sui piani principali in accordo con la corrispondente specie, la collocazione dei dati in dipendenza dalla specie.

Ripetere una completa analisi delle componenti principali della tabella 10euecosoc, alla quale però è stata rimossa la colonna IndustrialProduction. Confrontare i risultati con l’analisi della tabella completa.