Modelli autoregressivi per la stima della connettività tra segnali...

Modelli autoregressivi per la stima della connettivitàtra segnali biologici

Ing. Laura Astolfi

Dip. Informatica e Sistemistica

Fondazione S.Lucia IRCCS

E-mail: [email protected]

Corso di Modelli di Sistemi Biologici – Prof. S. Salinari

Università di Roma “La Sapienza” - Facoltà di Ingegneria

L. Astolfi - Modelli autoregressivi per la connettività tra segnali biologici - Corso di Modelli di Sistemi Biologici 2

Piano delle lezioni

•Definizione di connettività tra segnali biologici• Connettività anatomica, funzionale ed effettiva• Definizioni di causalità

•Modelli autoregressivi (AR)• Identificazione dei parametri di un AR: metodo di Yule-Walker• Impiego degli AR come predittori lineari

•Modelli bivariati autoregressivi•Metodi bivariati per la stima della connettività

• Test di causalità di Granger• Coerenza Ordinaria• Coerenza Diretta

•Modelli multivariati Autoregressivi (MVAR)•Metodi multivariati per la stima della connettività

• Directed Transfer Function• Partial Directed Coherence• Direct Directed Transfer Function

•Metodiche non parametriche per la stima della connettività•Confronto tra metodi bivariati e metodi multivariati


Connettività tra segnali biologici

• SEGNALI UNIVARIATI:•L’informazione biologica può essere contenuta in un unico segnale -> ANALISI UNIVARIATA(es: temperatura)

• SEGNALI MULTIVARIATI:•Più spesso, numerosi segnali biologici vengono registrati simultaneamente, e l’analisi dell’interdipendenza tra di essi può consentire di comprendere il funzionamento dei sistemi biologici che li hanno prodotti (ANALISI MULTIVARIATA) (es: segnali elettroencefalografici, elettrocardiografici, elettromiografici…)



• L’osservazione separata dell’andamento temporale di alcuni indici dell’attività biologica fornisce solo informazioni parziali sul funzionamento dei sistemi in esame

• Per comprenderne realmente il funzionamento è necessario capire anche come le diverse parti del sistema stesso interagiscano tra loro -> connettività funzionale



Es: Attività cerebrale• Può essere rilevata mediante

Elettroencefalografia (EEG), Magnetoencefalografia (MEG), Risonanza Magnetica Funzionale (fMRI), Local FieldPotentials (LFP),….


Applicazione: connettivitàcorticale

• Tutte queste metodiche restituiscono informazione sull’andamento temporale dell’attività nei vari punti della corteccia cerebrale (Mappe di attivazione)

• Ma l’informazione multicanale, da sola, non dice abbastanza sul funzionamento della corteccia

• Connettività corticale: descrive come le diverse aree della corteccia cooperano tra loro per l’esecuzione di un determinatocompito sperimentale


Diverse definizioni di connettività•Connettività anatomica= esistenza di legami anatomici che permettono la trasmissione di informazione da un distretto biologico ad un’altro (ad es: fasci nervosi)

•Connettività funzionale= esistenza di legami funzionali tra le attività rilevate in diversi distretti anatomici, indipendentemente dall’esistenza di una effettiva connessione anatomica diretta.

Definizione intuitiva: descrizione di come diversi sottosistemi di un sistema biologico (ad es: diverse aree della corteccia cerebrale, il cervello e i muscoli, il cervello e il cuore) comunichino tra loro. Può essere espressa dal concetto di CAUSALITA’ (un evento ne provoca un altro)

Ma per poterla stimare occorre una definizione matematica!


Definizione di causalità

Norbert Wiener (1956). Prima definizione di causalità in senso statistico:

Dati due segnali, misurati simultaneamente, se la predizione del primo è migliorata dall’aver incorporato le informazioni relative al passato del secondo segnale, piuttosto che usando solo l’informazione relativa al passato del primo, allora si può dire che il secondo segnale è causa (in senso statistico) del primo.


Causalità secondo Granger

Clive Granger (Premio Nobel per l’economia2003), 1969: Formulazione matematica del concetto espresso da Wiener.

Dati due segnali X e Y, si può dire che X influenza Y se includendo i valori passati della prima variabile nella stima regressiva della seconda, questo migliora la qualità della stima, cioè riduce l’errore di predizione.

nn--33 nntt

X(t)X(t)

Y(t)Y(t)


Modelli Autoregressivi (AR)

[ ] [ ] [ ] [ ]∑=

+−−=p

k

nuknxkanx1

x[n]= serie temporale

a[n]= parametro autoregressivo relativo al lag k

p= ordine del modello

u[n]= residuo del modello

Ipotesi:

x[n] processo stazionario in senso lato

u[n] bianco, a media nulla


Modello AR come filtro generatore

Hp:Filtro causale (h[n]=0 per n<0)

h[n]u[n] x[n]

Input: processo bianco

Output: processo x[n]

[ ] [ ] [ ]∑∞

=

−=0k

knukhnx


Stima dei parametri AR

Ipotesi:

1. x[n] processo stazionario in senso lato

2. la sequenza di autocorrelazione rxx[m] si suppone nota (in realtà per segnali finiti si può solo stimare)

3. u[n] bianco, a media nulla


Moltiplicando per x*[n-m]:

Calcolando il valore atteso:

O, equivalentemente:


[ ] [ ] [ ] [ ]∑=

+−−=p

k

nuknxkanx1

[ ] [ ] [ ] [ ] ][*][*][*1

mnxnumnxknxkamnxnxp

k

−+−−−=− ∑=

[ ] [ ] [ ] [ ] ]}[*{]}[*{]}[*{1

mnxnuEmnxknxEkamnxnxEp

k

−+−−−=− ∑=

[ ] [ ] [ ]mrkmrkamrp

kuxxxxx ∑

=

+−−=1

][



[ ] [ ] [ ] [ ] [ ] [ ]∑∑∞

==

−=+−−=01 k

p

k

knukhnuknxkanx


kuxxxxx ∑

=

+−−=1

][

[ ] [ ] [ ]{ } [ ] [ ] [ ] =⎭⎬⎫

⎩⎨⎧

⎥⎦

⎤⎢⎣

⎡−−=−= ∑

∞

=0

***

kux kmnukhnuEmnxnuEmr

[ ] [ ] [ ] [ ] [ ]∑∞

=

−∗+=++=1k

uu*

uuuu*

uu mrmh]m[rkmrkhmr

[ ] [ ] [ ] [ ] =⎭⎬⎫

⎩⎨⎧

⎥⎦

⎤⎢⎣

⎡−−+−= ∑

∞

=1

***

k

kmnukhmnunuE

Poiché:

Sostituendo nella equazione precedente :



Poiché u è per ipotesi bianco e a media nulla, indicando con σ2 la sua

varianza si ha:[ ][ ] 20

00

σ=

≠∀=

uu

uu

r

mmr

[ ][ ]

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧

>=<−

=00

0

02

*2

m

m

mmh

mrux σσ

[ ] [ ] [ ] [ ] [ ] [ ]∑∞

=

−∗+=++=1k

uu*

uuuu*

uuux mrmh]m[rkmrkhmrmr

Inoltre, poiché il filtro è causale, h[m]=0 per m<0 (i campioni futuri

delle serie in ingresso sono incorrelati con il campione attuale della

serie in uscita)

diventa:

[ ] 0mkruu ≠+ 2σe = solo se k=-m

se m>0 allora m+k diverso da 0



[ ][ ]

[ ] [ ]

[ ] [ ] ⎪⎪⎪

⎭

⎪⎪⎪

⎬

⎫

⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

>−−

=+−−

<−

=

∑

∑

=

=

0

0

0

1

2

1

*

mkmrka

mkrka

mmr

mr

p

kxx

p

kxx

xx

xx σ


kuxxxxx ∑

=

+−−=1

][ [ ][ ]

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧

>=<−

=00

0

02

*2

m

m

mmh

mrux σσ

Inoltre, per la stazionarietà in senso lato di x si ha che

rxx[m]= rxx[-m]

EQUAZIONI DI YULE-WALKER


Equazioni di Yule-Walker

• Sono p+1 equazioni-> consentono di ricavare p parametri a[n] e la varianza del rumore in ingresso σ2:

[ ] [ ] [ ][ ] [ ] [ ][ ] [ ] [ ]

[ ] [ ] [ ]

[ ][ ]

[ ] ⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

=

⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

∗

⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

−

+−+−

−−

0

0

0

2

1

1

01

212

101

10 2

MM

LL

MMM

L

L

L σ

pa

a

a

rprpr

prrr

prrr

prrr

xxxxxx

xxxxxx

xxxxxx

xxxxxx

var =⋅In forma matriciale:


Modelli Autoregressivi (AR)

Il modello AR può essere impiegato:

• come un filtro generatore, ponendo in ingresso una sequenza pilota (rumore bianco) e ottenendo in uscita un processo con le stesse proprietà del segnale da stimare;

• come un predittore lineare, che consente di predire il campione al passo n-esimo a partire dai campioni ai passi precedenti


Predizione lineare autoregressiva

• Un filtro AR può essere impiegato come predittore lineare (è questo l’uso nella connettività):

[ ] [ ] [ ]∑=

−−=p

k

knxnanx1

ˆ

[ ] [ ] [ ]nxnxne −=

X[n-6]X[n-5]

X[n-4]

X[n-3]X[n-2]

X[n-1]X[n][ ]nX• L’errore di

predizione sarà:

• Lo scopo della predizione lineare èdeterminare i coefficienti a[k], minimizzando la potenza dell’errore e[n]



Ipotesi:

1.x[n] processo stazionario in senso lato

2.Ortogonalità tra le coppie di serie e[n] e x[n]. Questo per rispondere alla condizione di minimo della potenza dell’errore. Infatti se queste non fossero scorrelate significherebbe che l’errore contiene in sé ancora delle informazioni utili sul segnale.



• Si dimostra che i coefficienti del vettore di predizione lineare che minimizzano la varianzaσ2 sono dati dalla soluzione delle equazioni:

[ ]⎟⎟⎠

⎞⎜⎜⎝

⎛=⎟⎟

⎠

⎞⎜⎜⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛

− 0

10 2

1

σaRr

rr

ll

Hlxx Di nuovo le eq.

di Yule- Walker!

var =⋅In forma matriciale:


Modello AR e predittore AR

• Differenza tra un modello AR e un predittore AR:

h[n]u[n] x[n]

h[n]x[n] e[n]

Input: processo bianco

Output: processo x[n]

Input: processo x[n]

Output: errore e[n]

•e[n] è scorrelato dal segnale stimato

•In generale, non è detto che l’errore e[n] sia un processo bianco

•Se però il processo x[n] è autoregressivo nella sua natura, e modellizzabile con un filtro AR dello stesso ordine scelto per il predittore, allora e[n] è bianco.

][ˆ nx


Applicazione alla connettività

• Come visto, un modello predittore di tipo AR consente, una volta che i suoi parametri siano stati stimati, di predire il campione all’istante n, x[n], in funzione dei campioni nei p istanti precedenti x[n-1],x[n-2],…,x[n-p].

X[n-3]X[n-2]

X[n-1]X[n][ ]nX

Y[n-3]Y[n-2]

Y[n-1]X[n][ ]nY

•Se si potesse estendere questo modello a 2 serie temporali, includendo i campioni agli istanti precedenti di una serie nella predizione dell’altra…?


Modelli bivariati

• Il campione di x è predetto inserendo anche i campioni agli istanti precedenti di y, e viceversa.

• Anche l’errore di predizione è quindi influenzato da entrambi i segnali

• I parametri del modello vengono stimati con un analogo bivariato delle equazioni di Yule-Walkerin cui compaiono le funzioni di cross-correlazione rxy e ryx

[ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ]∑ ∑

∑∑

= =

==

+−+−=

+−+−=

p

1k

p

1kyxyxyx

p

1kxyxy

p

1kxy

neknykbknxkany

neknykbknxkanx


Test di causalità di Granger

• Ora si può applicare la definizione di Granger: X influenza Y se includendo i valori passati della prima variabile nella stima regressiva della seconda, questo riduce l’errore di predizione.

• Confrontando i modelli univariati e bivariati stimati su x e y:

[ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ]∑ ∑

∑∑

= =

==

+−+−=

+−+−=

p

k

p

kyxyxyk

p

kxyxy

p

kxy

neknykbknxkany

neknykbknxkanx

1 1

11

[ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ]∑

∑

=

=

+−=

+−=

p

kyy

x

p

kx

neknykany

neknxkanx

1

1Caso univariato

Caso bivariato



• La capacità predittiva dei due modelli può essere espressa dalla varianza degli errori di predizione:

( )( )yyy

xxx

eV

eV

var

var

=

= ( )( )yxx,yy

xyy,xx

evarV

evarV

=

=Per i modelli univariati

Per i modelli bivariati

Una misura dell’influenza di y su x può allora essere espressa da:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛=→

yxx

xx

xy V

VG

,

ln

Se y non migliora la predizione di x,

Se la migliora,

0, ≈⇒≈ GVV xxyxx

↑⇒↓ GV yxx ,



•Applicazione:Analisi della connettività tra segnali ottenuti su 5 elettrodi da registrazioni EEG.

X1

X2

X3

X4

X5

1. Su tutte le possibili coppie all’interno del set di segnali viene calcolato il coefficiente di Grangerin entrambe le direzioni (1->2 e 2->1)

2. Per ogni coppia può esistere un legame di tipo causale in un verso, nell’altro, in entrambi o in nessuno dei due

3. Risultato: un pattern di connettività che lega i 5 siti sullo scalpo



• Vantaggi del test di causalità di Granger:• Indica la DIREZIONE del legame di causalità:

• Limiti:• Dà indicazioni nel solo dominio del tempo

(nell’intervallo temporale analizzato)

• Spesso invece l’informazione associata ad un segnale biologico è codificata nel dominio della frequenza -> occorre una misura spettrale

yxxy GG →→ ≠


Nozioni di analisi spettrale per segnali tempo-discreti

• Trasformata di Fourier di un segnale tempo-discreto:

• Densità spettrale di energia:

• Funzione di autocorrelazione:

( ) ( )2xx fXfS =

( ) [ ] ( )fkT2jexpkxfXk

π−= ∑∞

−∞=

[ ] [ ] [ ]nkxkxnrk

*xx += ∑

∞

−∞=


Teorema di Wiener- Khinchin

• Lo Spettro Densità di Potenza (PSD) è calcolato come trasformata di Fourier della Funzione di Autocorrelazione:

• Nella pratica, sono disponibili sequenze discrete ottenute per campionamento a passo T:

( ) ( )[ ] ( )∫∞

∞−

−=ℑ= ττ τπ dertrfS f2jxxxxx

( ) [ ] ( )fkT2jexpkrfSk

xxxx π−= ∑∞

−∞=


Nozioni di analisi spettrale per segnali tempo-discreti

• Funzione di cross-correlazione:

• Densità spettrale di mutua potenza:

• Proprietà: Sxy = Syx

( ) ( ) ( )( ) ( ) ( )fXfYfS

fYfXfS*

yx

*xy

=

=

[ ] [ ] [ ]nkykxnrk

*xy += ∑

∞

−∞=


Matrice spettrale dei segnali

• Matrice spettrale dei segnali:

• In cui:

Sii= Densità spettrale di potenza di xi

Sij= Densità spettrale di mutua potenza tra xi e xj

( ) ( ) ( )( ) ( )⎥⎦

⎤⎢⎣

⎡=

fSfS

fSfSfS

2221

1211 E’ una matrice SIMMETRICA


Coerenza Ordinaria

• La Coerenza Ordinaria indica la correlazione lineare tra due segnali in funzione della frequenza.

• E’ data dalla mutua densità spettrale dei due segnali normalizzata sul prodotto delle loro densità spettrali individuali.

( )( )

( ) ( )fSfS

fSfC

yyxx

2

xy

xy =


Proprietà della coerenza ordinaria

• La coerenza assume valori compresi nell’intervallo [0, 1]

• Per una data frequenza f0:•Cxy(f0)=0 => l’attività dei segnali a questa

frequenza è indipendente.

•Cxy(f0)=1 => indica la massima correlazione a questa frequenza.


Coerenza ordinaria

• LIMITE: non indica il verso dell’interazione

• Infatti la coerenza ordinaria misura l’esistenza di un sincronismo tra i segnali, ma non è legata al concetto di causalità => non esprime la direzionalità

xx yy??Cxy

( )( )

( ) ( )( )fC

fSfS

fSfC yx

yyxx

2

xy

xy ==


Coerenza diretta

• Le grandezze che sostituiscono la matrice spettrale S possono essere ricavate mediante il teorema di Wiener- Khinchin :

• oppure con un modello di filtro generatore che dia informazioni sulla struttura dei segnali => AR bivariato

( ) [ ] ( )

( ) [ ] ( )fkT2jexpkrfS

fkT2jexpkrfS

kyyyy

kxxxx

π

π

−=

−=

∑

∑∞

−∞=

∞

−∞=


Coerenza diretta

• Infatti, dato un modello AR bivariato si ha:

• Da cui, portando le x al primo membro:

[ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ]∑ ∑

∑∑

= =

==

+−−−−=

+−−−−=

p

1k

p

1kjjjjjijij

p

1kijjij

p

1kiiii

neknxkaknxkanx

neknxkaknxkanx

[ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]∑ ∑

∑ ∑

= =

= =

=−+−

=−+−

p

1k

p

0kjjjjjiji

ij

p

0k

p

1kjijiii

neknxkaknxka

neknxkaknxka Con

aii[0]=1

ajj[0]=1.


Coerenza diretta

• Passando al dominio della frequenza otteniamo:

• Dove:

( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )fEfXfAfXfA

fEfXfAfXfA

jijjjiji

ijjijiii

=+

=+

( ) [ ]∑=

−=p

0k

fTk2jijij ekafA π

è la trasformata di Fourier dei coefficienti del modello.


Coerenza diretta

• In forma matriciale si ha:

( ) ( ) ( )fEfXfA =

( ) ( ) ( )( ) ( )⎥⎦

⎤⎢⎣

⎡=

fAfA

fAfAfA

jjji

ijii

( ) ( )( )⎥⎦⎤

⎢⎣

⎡=

fX

fXfX

j

i

( ) ( )( )⎥⎦⎤

⎢⎣

⎡=

fE

fEfE

ji

ij

Con:


Coerenza diretta

• L’equazione precedente può essere riscritta come segue:

dove H(f) è detta MATRICE DI TRASFERIMENTO del fitroAR:

( ) ( ) ( ) ( ) ( )fEfHfEfAfX 1 == −

( ) ( ) ( ) ( )( ) ( )⎥⎦

⎤⎢⎣

⎡== −

fHfH

fHfHfAfH

jjji

ijii1


Matrice di trasferimento

L’elemento ij della matrice H è la funzione di trasferimento tra l’i-esimo ingresso e la j-esima uscita del filtro AR bivariato. In generale si ha che H ij (f) è diverso da H ji (f)

E1(f)

X2(f)E2(f)

X1(f)

H22(f)

H11(f)

H12(f) H

21(f)

FILTRO AR bivariato


Coerenza diretta

• Passando alle densità spettrali di potenza, si ottiene:

• Poiché la matrice spettrale S espressa in funzione dei parametri del modello AR è quindi data da:

( ) ( ) ( ) ( ) ( ) ( ) ( )fHfEfEfHfXfXfS HHH ==

( ) ( ) ( )( ) ( ) ( )[ ]

( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ( )( ) ( ) ( )fS

fSfS

fSfS

fXfXfXfX

fXfXfXfX

fXfXfX

fXfXfX

jjji

ijii

*jj

*ij

*ji

*ii

*j

*i

j

iH

=⎥⎦

⎤⎢⎣

⎡=

⎥⎥⎦

⎤

⎢⎢⎣

⎡

=⎥⎦

⎤⎢⎣

⎡=

( ) ( ) ( )fEfHfX =


Coerenza diretta

• è la matrice spettrale di e1[n], e2[n]

• Ma, per Hp, e1[n], e2[n] sono bianchi, quindi:

1) la loro densità spettrale di potenza è costante e pari alla loro potenza σ11

2 , σ112

2) La loro densità di mutua potenza spettrale è costante e pari a σ12

2= σ212

• La matrice spettrale di E è pertanto

indipendente dalla frequenza e pari a:( )

⎥⎥⎦

⎤

⎢⎢⎣

⎡== 2

jj2ji

2ij

2iifE

σσσσ

Σ

⎥⎥⎦

⎤

⎢⎢⎣

⎡= 2

jj

2ii

0

0

σσ

ΣInoltre, se il modello fornisce una buona stima del processo X, Eij e Eii sono scorrelati=> la loro mutua correlazione è zero =>

( ) ( ) ( )fSfEfE EH =


Coerenza diretta

• La matrice spettrale S espressa in funzione dei parametri del modello AR è quindi data da:

( ) ( ) ( ) ( ) ( ) ( ) ( )fHfHfHfEfEfHfS HHH Σ==

( ) ( ) ( )fHfHfS HΣ=

⎥⎥⎦

⎤

⎢⎢⎣

⎡= 2

jj

2ii

0

0

σσ

ΣCon:


Coerenza diretta

• Sulla base della fattorizzazione della matrice spettrale S fornita dalla equazione precedente, si può definire la COERENZA DIRETTA da j verso i:

( ) ( )( )

( )( ) ( )2ij

2jj

2

ii2ii

ijjj

ii

ijjjij

fHfH

fH

fS

fHf

σσ

σσγ

+==


Coerenza diretta

• Poiché , si ha che:

A differenza della Coerenza Ordinaria, la Coerenza Diretta fornisce indicazioni sulla direzione dell’interazione tra i due segnali

( ) ( )ff jiij γγ ≠( ) ( )fHfH jiij ≠


Coerenza diretta

• Quindi per la Coerenza Diretta si ha:

xjxj xixi

( )fijγ

xjxj xixi

( )fjiγ

Come per il coefficiente di Granger => si dimostra che la Coerenza Diretta è un equivalente in frequenza del test di Granger


METODI BIVARIATI

• Test di causalitàdi Granger:

• Coerenza Ordinaria:

• Coerenza Diretta:

ijGxjxj xixi

jiG

xjxj xixi

ijγ

jiγ

xjxj xixiijC

Nel tempo

In frequenza

In frequenza

Modelli autoregressivi per la stima della connettività tra segnali...

Documents

Transcript of Modelli autoregressivi per la stima della connettività tra segnali...