Statistica multivariata (principal component analysis?)

wedge
ciao, vorrei capire empiricamente qual è la strategia per risolvere un determinato problema.

Ho 4 set di dati. Li chiamerò y, $x_1$, $x_2$ e $x_3$
y ha una blanda correlazione coni vari $x_i$.

Vorrei vedere se la correlazione è più forte tra y e una combinazione lineare degli $x_i$.

In altre parle, vorrei trovare i valori di a, b e c

$y = a * x_1 + b * x_2 + c * x_3$

che migliorano la correlazione, se essa esiste. (è una regressione lineare a più dimensioni, alla fine)
qual è la strategia da seguire? centra qualosa con la principal component analysis?

brutalmente, mi interessa un algoritmo operativo ;)

Risposte
olaxgabry
Ciao,
con l'analisi in componenti principali puoi ricavarti tre componenti del tipo

$z_{i}=a_{i}x_{1}+b_{i}x_{2}+c_{i}x_{3}$ $i=1,2,3$

però secondo me le cose non si miglioreranno di molto. Se ti interessa c'è un software opens source che te le crea in pochissimo tempo, si chiama gretl (in realtà anche R lo fa ed anche assai meglio).

wedge
grazie ad entrambi.

"Sergio":

Lascio perdere considerazioni più sofistiche e assumo che si tratti di un esperimento: hai un fenomeno che misuri in $y$ e tre sue possibili "spiegazioni" $x_1,x_2,x_3$ (fattori che possono influenzare i valori di $y$).
Le $x_i$ sono "blandamente" correlate con $y$. Si tratta di correlazioni "a coppia": $y$ con $x_1$, $y$ con $x_2$, $y$ con $x_3$.
Può la correlazione aumentare se si considera una combinazione lineare delle $x_i$?

esatto, il problema è proprio quello



Se esegui una regressione di $y$ su $x1$ e $x2$ ottieni:
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.3750     4.7405   0.079 0.940016    
x1            5.3750     0.6638   8.097 0.000466 ***
x2            9.2500     1.3276   6.968 0.000937 ***



interpreto bene? qui hai fatto (con R?) una regressione del tipo
$y = a \cdot x_1 + b \cdot x_2 + c$ (1)
esatto? nello spazio 3d dei valori osservati (1) geometricamente rappresenta un piano

ora ti dico cosa credo centri la PCA. (forse ho interpretato tutto male dall'inizio, può anche darsi).
pensiamo ancora ai dati distribuiti in uno spazio 3D.
trovare gli autovettori della matrice di covarianza, ridurre la dimensionalità da 3 a 2, non equivale a definire il piano (1) ?

in altre parole direi che l'autovalore più piccolo corrisponde alla direzione dove la varianza è più piccola, e quindi alla normale del piano che meglio fitta i dati, piano che posso ricavare facilmente con un conto di geometria a partire dalla normale.

chiaramente, se per contrario avessi tre autovalori uguali, ognuno dei quali pesasse per il 33% della varianza, ì dati sarebbero distribuiti come una sfera.

mi segui su questo ragionamento?

wedge
visto che mi piace il metodo sperimentale, ho provato a far girare sugli stessi tuoi dati la procedura che ho descritto per la PCA,
implementata in IDL

il risultato è
Eigenvectors: 
Mode#1    0.1275    0.0056    0.0131
Mode#2    0.0845    0.5980   -1.0734
Mode#3   -0.6113    6.0925    3.3459

Reconstruction error:   1.95968e-11

     Mode   Eigenvalue  PercentVariance
       1      60.7309      98.8931
       2     0.659229      1.07347
       3    0.0205393    0.0334457


Prendo il terzo eigenvalue e scrivo l'equazione del piano
$ y * -0.611306 + x_1 * 6.09250 + x_2 * 3.34588 = cost$

che equivale a:
$ y = 9.97 * x_1 + 5.47 * x_2 + cost$

(problema, da qui non saprei come valutare l'errore sui parametri)

che mi sembra simile a quella che trovi tu, come pure la correlazione tra $y$ e $k= 9.97 * x_1 + 5.47 * x_2 $
che viene $R(y,k) = 0.978402$

olaxgabry
"Sergio":

PS @olaxgabry:
a) davvero gretl si può usare anche per la PCA? Non mi pare. Dipende dal fatto che lo sto usando... solo da ieri?
b) la PCA con R non è poi così semplice. Non per altro, ma le due funzioni "standard", prcomp() e princomp(), non sono equivalenti e non fanno tutto. Alla fine, ho preferito usare una mia funzione derivata dalla principal() della libreria psych.
c) Se, partendo da tre variabili, arrivi a tre componenti principali, il beneficio è nullo e ottieni dimensioni più difficilmente interpretabili ;-)


Su gretl: se vai su visualizza trovi il tutto. L'unica cosa che non mi convincono sono i pesi.
Su R: vero, ma in fondo R non è semplicissimo. Io usao factoMineR per fare la ACP.
Sul punto c. Sono d'accordo in parte e ti faccio un esempio con cui mi imbatto ogni giorno. Quando si hanno due serie storiche non stazionarie ed $I(1)$ si vuole sapere se esiste una combinazione lineare stazionaria; in letteratura ci sono tanti test per vedere se tale relazione esista o meno: a mio parere, però, una cosa molta vantaggiosa sarebbe quelle di costruirsi le due ACP e vedere se una di queste è stazionaria. Chiaramente se le due componenti sono non stazionarie non vuol dire che non esista una relazione di cointegrazione. Fatto sta che ogni volta, nei casi analizzati da me, che i test mi dicono che tale relazione esiste guarda caso una componente ACF è stazionaria. Questo mi risolve molti problemi perché lavorando spesso con i modelli autoregressivi vettoriali (VAR), nel caso che verifichi l'esistenza di tale relazione, posso lavorare direttamente con le serie originali (non con le differenze) perché mantengo la consistenza dei parametri (non la normalità, ma del rapporto t sui VAR me ne frego).
Poi un'altra applicazione, su dati cross section, sarebbe quella di evitare il problema della multicollinearità nella regressione: insomma, dipende da quello che si vuole fare.
Scusa se mi so dilungato, ma appena sento odore di serie storiche mi dilungo spesso :).

wedge
ok, la mia grezzezza in statistica inizia a farsi sentire.

"Sergio":
[quote="wedge"]qui hai fatto (con R?) una regressione del tipo
$y = a \cdot x_1 + b \cdot x_2 + c$ (1)
esatto? nello spazio 3d dei valori osservati (1) geometricamente rappresenta un piano

Prima nota: se hai $n$ osservazioni (8 nell'esempio che ti proponevo), $Y$ è un vettore di $n$ elementi, quindi appartiene a $RR^n$. Non ti stai muovendo in uno spazio di dimensione 3; cerchi, invece, di arrivarci (v. dopo).
[/quote]

si, hai ragione. avendo 8 osservazioni di 3 variabili (l,m,n) sono molto abituato a pensare la cosa come 8 punti disposti in uno spazio (l,m,n), è questo quello che intendevo.


Poi... capiamoci. Cosa intendi per $c$? Potresti intendere un'intercetta, oppure un errore. Ovvero, il modello è:
$y=beta_0 + beta_1 x_1 + beta_2 x_2 + epsilon$


intendevo l'intercetta.


La regressione calcola stime dei coefficienti $beta_i$ e ti consente di calcolare valori "teorici" di $y$, che sono dati da $\hat y=b_0 + b_1 x_1 + b_2 x_2$, senza $epsilon$ (l'errore), dove i $b_i$ sono stimatori dei $beta_i$.


bene, fino a qui ci sono.


.Nella PCA (in cui $Y$ non compare), provi a proiettare lo spazio generato dalle colonne di $X$ su uno spazio di dimensione minore. Ovviamente lo fai lavorando sugli autovalori:


invece io qui tratterei la Y al pari di $X_1$ e $X_2$
per me sono tre misure che effettuo indipendentemente, e solo per scelta arbitraria quando invece penso alla regressione, scrivo Y in funzione delle $X_i$
e quindi trovo un valore teorico per la Y.

quando con la PCA cerco di ridurre le dimensioni, è perchè cerco di passare dallo spazio delle osservazioni [Y,$X_1$,$X_2$] a due vettori che mi definiscano il "piano più favorevole" per la descrizione dei dati (attraverso l'autovettore con autovalore più piccolo, che è chiaramente la normale a tale piano)

forse non mi sono spiegato bene sul significato che do ai diversi set, che vedo tutti equivalenti tra loro.
è come se facessi osservazioni indipendenti. esempio: misuro peso P, lunghezza L e velocità massima V di 8 autovetture.
poi cerco se e come sono correlate.
in una regressione arbitrariamente potrei scrivere il peso in funzione delle altre due P= a*L +b*V + c, mentre con le PCA riduco di uno la dimensionalità di [P,L,V].

in definitiva, per un problema del genere tu favoriresti l'approccio regressione a quello PCA?

ingenuamente mi sembrano equivalenti entrambe a cercare il piano $k_1*P+k_2*L+k_3*V = cost$ che minimizza la somma delle distanze dei punti osservati $[P_n, L_n, V_n]$ dal piano stesso.

grazie per la pazienza.

wedge
"Sergio":

Onestamente, non ho capito cosa hai fatto. Ma una cosa è certa: se prendo una matrice di dati le cui tre colonne sono $y,x_1,x_2$, e se chiedo tre componenti principali, ottengo tre nuovi vettori tra loro incorrelati.
Una cosa un po' diversa dall'obiettivo, che era: trovare una correlazione meno blanda tra $y$ e una combinazione lineare delle variabili esplicative.


certo, sono 3 vettori tra loro linearmente indipendenti.
ma non vedo perchè sia diverso dall'obiettivo, visto che di quei tre vettori i primi 2 mi rappresenteranno le direzioni lungo le quali è massima la varianza, ossia sono sufficienti a definirmi il piano ove i dati tendono a disporsi.
ricavata l'equazione di quel piano, essa mi definisce la combinazione lineare delle variabili che minimizzano la correlazione con y.

olaxgabry
"Sergio":
[quote="wedge"]

La PCA ti dice "solo" che, se un fenomeno è descritto (non spiegato, descritto) da tre variabili, e se trovi 2 componenti (combinazioni di quelle variabili) che danno conto di un'alta percentuale della variabilità, puoi descrivere quel fenomeno ricorrendo a 2 componenti (entrambe correlate a tutte e tre le variabili di partenza) invece che a 3 variabili. Ti rimane poi il problema -- da non sottovalutare -- di interpretare le componenti [inciso: in ambito sperimentale, la regressione soffre della cosiddetta multicollinearità, cioè di una correlazione non perfetta ma non trascurabile tra le variabili esplicative; ciò nonostante, si ricorre raramente ad una PCA preliminare, per poi eseguire la regressione su variabili incorrelate, proprio perché risulta spesso difficile interpretare le componenti].
[/quote]

Se ti può interessare, io spesso la uso così in ambito di regressione: magari mi dai dei suggerimenti o discutiamo delle conseguenze visto che sei come al solito ferratissimo in materia :).
Sull'uso dell'ACP per evitare la multicollinearità: se devo essere onesto, io la uso spesso perché voglio avere un gruppo di modelli plausibili per cui non vedo perché non metterci anche uno o più modelli con regressori ACP. Poi valuto i modelli con dati fuori dal campione e vedo come si comportano. Ora, se me ne fregassi dei test sui parametri (con il rischio di modelli sovraidentificati che causano overfitting) probabilmente non dovrei pormi troppi problemi sulla multicollinearità che in fondo incide sulla varianza degli stimatori. Con l'ACP risolverei i problemi della multicollinearità e potrei senza problemi effettuare test sui parametri.
Tutto dipende dallo scopo dell'analisi, io comunque tendo sempre a considerare la ACP nelle regressioni. Ovvio che se uno volesse valutare come incidono dei regressori sulla variabile dipendente, non so se utilizzando le componenti principali si possa soddisfare la richiesta nonostante io possa sempre definire il modello in base ai regressori iniziali: però i pesi che avrei alla fine per ogni regressore di partenza sono effettivamente interpretabili per valutarne l'incidenza? Qui non ti nascondo delle perplessità.
Sarà deformazione professione, ma è una delle parti della statistica che utilizzo di più soprattutto nelle serie storiche.

olaxgabry
Interpretare le componenti è la cosa sicuramente più difficile. Supponi che abbia tre regressori $x_{i}$ e che ci costruisca le 3 componenti principali $PC_{i}$ e definisca la regressione

$y_{i}=aPC_{1}+bPC_{2}+cPC_{3}+u_{i}$

Stimati i parametri del modello, supponendo per comodità che questi siano tutti significativi, posso ricondurmi ad un modello

$y_{i}=c_{1}x_{1}+c_{2}x_{2}+c_{3}x_{3}+u_{i}$

in quanto le $PC_{i}$ sono combinazioni lineari delle $x_{i}$. Qui ho il dubbio se adesso questi pesi $c_{i}$ siano effettivamente utilizzabili per fini di interpretazione del modello ovvero come incidono i regressori sulla $y$.

Sulle serie storiche: la cointegrazione è molto interessante. Però i test hanno spesso delle ipotesi forti (johansen): se ad esempio venisse meno l'ipotesi di normalità il test di Johansen ha seri problemi. Poi qui si aprono tanti discorsi: ad esempio, in caso di cointegrazione, è utile usare i modelli a correzione dell'errore vettoriale VECM? Gli economisti gli danno molta importanza a questi modelli (funzioni impulso risposta), a me però interessano di più le previsioni e qui un classico VAR sulle serie nei livelli è perfetto (c'è un bellissimo articolo di Stock-Sim-Watson a riguardo, quando vuoi te lo giro).
Non vedo l'ora di discutere su questi argomenti allora ;-).

olaxgabry
Hai ragione sullìOT e mi scuso con wedge.
Concludo: io penso che se mettono molti regressori crea non pochi problemi ai modelli, quindi mi dissocio da questa corrente di pensiero tipica di alcuni econometrici. Io sarò limitato, ma per il momento mi bastano i VAR.
Su Granger: qui c'è molto, diciamo che parliamo di colui che ha sviluppoto gran parte dell'econometria come cointegrazione o casualità secondo Granger. Io ho usato la Granger casuality, intesa come miglioramento nella previsione: dato che dico previsione, la applico con dati out of sample. I test che si trovano in quasi tutti i libri per me non ci dicono niente sulla previsione, magari sulla stima che è ben diverso.
Dai, la chiudo qui altrimenti mi dilungo troppo. Appena vuoi ne riparlerò volentieri magari aprendo dei topic apposta.

wedge
rieccomi.

riguardo agli autovettori diversi, è possibile che abbia sbagliato a implementare il mio codice.
"Sergio":

Allora: sei partito da tre vettori tra loro non ortogonali. Hai poi ottenuto altri tre vettori tra loro ortogonali. Guardi i relativi autovalori e vedi che i primi due spiegano quasi tutta la variabilità. Concludi quindi che il sottospazio (un piano) generato dai primi due autovettori racchiude il "segnale", mentre il terzo, ortogonale al piano, è "rumore". Detta per immagini, hai trovato che quella che poteva sembrare una generica nuvola di punti (un ellissoide), in realtà è molto piatta, come un osso di seppia; se ti muovi lungo gli assi di questo osso di seppia trovi quasi tutta l'informazione racchiusa nei tuoi dati.
A questo punto, che senso ha costruire l'equazione di un piano usando gli elementi del terzo autovettore?


un piano è definito a partire da un punto $(x_0,y_0,z_0)$ e dal vettore ad esso normale $(nx, ny, nz)$, come
$(x-x_0) * n_x + (y-y_0) * n_y + (z-z_0) * n_z = 0$
credo sia il modo più semplice per scriverlo in questo caso.


Dici poi che ottieni un'equzione $y=9.97*x_1+5.47*x_2+cost$ molto simile a quella ottenuta con la regressione. In realtà, però, questa era:
$y=5.375*x_1+9.25*x_2+cost$
Mica tanto simile! I coefficienti si somigliano, ma sono scambiati di posto.
Un conto è dire che l'aumento di una unità di $x_1$ porta ad un aumento di $9.97$ unità di $y$, tutt'altro è dire che porta un aumento di $5.375$.


hai ragione sullo scambio di posto (sto facendo troppe cose insieme :oops: ), chiedo perdono!

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.