1. No category

FYS 4780
Analyse av diagnostisk, dynamisk bildeinformasjon
Del 3
Modellering av dynamiske data
Atle Bjørnerud, Rikshospitalet
[email protected]
975 39 499
Modellering av dynamiske data
• Behov for å tilpasse data til modell
• lineær og ikke-lineær kurvetilpasning
• diskrete løsninger på kontinuerlige problemer:
– foldings-integral / dekonvolusjon
Når trenger vi å modellere data?
• Estimering av spesifikke parametre av
interesse fra målbare observasjoner:
– relaksasjonsparametre
– kinetikk-parametre
– perfusjons-parametre
Regresjons-analyse
lineær regresjon:
y=Ax + b
1.2
y = 0.0104x - 0.0452
R² = 0.8074
1
0.8
0.6
0.4
0.2
0
0
-0.2
20
40
60
80
100
120
Regresjons-analyse
Ikke-lineær (eksponensiell) regresjon:
y=A*exp(-TE*R2)
800
SI
700
600
y = 750.53e-0.019x
R² = 0.9819
500
400
300
200
100
0
0
50
100
150
200
250
300
ekkotid (ms)
Regresjons-analyse
Ikke-lineær (eksponensiell) regresjon:
y=|A*(1-B*exp(-TI*R1))|+C
ROI signal intensity
850
g
b
c
d
e
f
b
c
d
e
f
g
800
750
700
650
600
550
500
450
400
350
300
250
200
150
100
50
200
400
600
Inversion Time [ms]
800
ROI #1
T1 relaxation f it
Regresjons-analyse
Ikke-lineær (gamma-variat) regresjon:
ROI signal intensity
18
g
b
c
d
e
f
b
c
d
e
f
g
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0
0
5
10
15
20
25
30
Image Number
35
40
45
ROI #1
Gamma var f it
Datamodellering
• Hva bestemmer om en kurvetilpasning er ’bra’
eller ’dårlig’?
• Trenger å definere en ’merit-funksjon’ (’figure
og merit’) som måler grad av overenstemmelse mellom observasjoner og data
Datamodellering – hva er
problemstillingen?
• Vi har et visst antall (N) målepunkter: (xi,yi) i=1..N
• Vi har en modell som (vi tror) predikerer et funksjonelt forhold
mellom de målte avhengige (y) og uavhengige (x) variabler slik
at: y(x) = y(x;a1,...aM)
• Spørsmålet vi ønsker å besvare er: gitt et sett med parametre;
hva er sannsynligheten for at det målte datasettet kunne
oppstå?
• Merk: vi kan ikke besvare det ’motsatte’ spørsmålet: Hva er
sannsynligheten for at et visst sett med parametre a1..aM er
’korrekt’ gitt observasjonene y(x). Hvorfor? : det finnes ikke et
statistisk univers med modeller; bare en: den korrekte...
Datamodellering – hva er
problemstillingen?
• Mål: identifisere sannsynligheten for
dataobservasjonene gitt modellen
• Hva menes med ’sannsynlighet’ i denne
sammenheng? det er kun en statistisk intuisjon
av sannynligheten (likelihood) for sammenheng
er og har egentlig ikke noe matamatisk
fundament...
• Med dette utganspunkt ønsker vi nå å tilpasse
parameterne a1...aM slik at vi får en ’maximum
likelihood’ for dataobservasjonen gitt modellen
y(x;a1,...aM).
Minste-kvadrat (least-squares) som
merit-funksjon
• Basert på forutsetningen om at hvert
datapunkt yi er assoisert med en uavhengig,
tilfeldig målfeil som er normalfordelt rundt
den ’sanne’ modellen y(x).
• Dersom standard-avviket (σ) i
normalfordelingene for hvert målpunkt er
konstant så er sannsynligheten for et målt
datasett lik produktet av sannsynligheten for
hvert målepunkt:
Minste-kvadrat som merit-funksjon
• Sannsynligheten for et målt datasett:
2

  1  yi  y ( xi )   

P   exp  

y
  
2

  
i 1 
  

N
∆y =måle-usikkerhet pga ikke-kontinuerlig parameter-rom (= konstant)
yi = målt datapunkt
y(xi)=modellert datapunkt, basert på modell-parametre a1,...aM
σ = standardavik i normaldistribusjonen for hvert målpunkt (=konstant)
Vi ønsker nå å maksimere P som er det samme som å minimere -log(P):
  yi  y ( xi )  2

min  
  N log y   min


 i 

N
2


y

y
(
x
;
a
..
a
)
 i
i
1
M
i 1
Chi-kvadrat som merit-funksjon
• Minstekvardats minimering forutsetter at hvert
målpunkt er normalfordelt rundt den ’sanne’
verdien og at denne normalfordelingen er lik
(samme varians) for alle datapunkter
• Dersom hvert datapunkt har sin egen, kjente
varians σi2, er minimeringsproblemet beskrevet
ved å minimere chi-kvadrat uttrykket:
 yi  y ( xi ; a1..aM ) 

   
i
i 1 

N
2
2
Chi-kvadrat minimering for lineær
regressjon
y ( x)  y ( x; a, b)  a  bx;
 yi  a  bxi 
2

 (a, b)   
i
i 1 

N
2
N
N
 2
 2
 yi  a  bxi 
 yi  a  bxi 
 2 

0
;


2


0

2
2
a
 i
b
 i


i 1 
i 1 
2 likninger med 2 ukjente
S 
i 1
aS  bS x  S y
1
 2i
N
S xx  
aS x  bS xx  S xy
a
N
i 1
S xx S y  S x S xy
S xx  S x
2
;b 
N
; Sx  
i 1
x 2i

2
i
xi
 2i
N
; S xy  
i 1
S xy  S x S y
S xx  S x
2
N
;Sy  
xi yi
 2i
i 1
yi
 2i
Lineær regressjon – estimering av
usikkerthet i tilpasning
y ( x)  y ( x; a, b)  a  bx;
 yi  a  bxi 

i
i 1 

N
 2 (a, b)   
2
2  f 
 f    i  
i 1
 yi 
N
’propagation of error’:
2
2
varians i estimerte parametre a,b (rett linje):
S  S x xi
S xx
a
2
 2 xx



a
yi  i ( S xx  S x2 )
S xx  S x2
Sxi  S x
b
S
2
 2
 b 
2
yi  i ( S xx  S x )
S xx  S x2
N
S 
i 1
1
 2i
N
S xx  
i 1
N
; Sx  
x 2i
 2i
i 1
xi
N
 2i
N
; S xy  
i 1
;Sy  
xi yi
 2i
i 1
yi
 2i
regresjons-koeffisient (R2)
(’coefficient of determination’)

y ( x)  y ( x; a, b)  a  bx  y;
Definisjoner:
y = middelverdier for målte data data

y = modelldata
’Total sum of squares’:
’Regression (explained) sum of squares’:
SStot   ( yi  y ) 2
i

SS reg   ( yi  y ) 2
i
’Resdiual sum of squares’:

SSerr   ( yi  yi ) 2
i
Regresjons-koeffisient:
R2  1
S err
S tot
=1-(’unexplained residue’/’total residue’)
R2 kan tolkes som hvor stor del av en respons-variasjon som kan forklares av
regressorene i modellen:
f,eks: R2=0.7 =>70% av respons-variasjonen kan forklares av modellen; resterende 30% er
forårsaket av enten ’ukjente’ regressorer eller av usikkerhet (støy) i datamålingene.
Regresjons-analyse
lineær regresjon:
y=ax + b

SSerr   ( yi  yi ) 2
i
1.2
y = 0.0104x - 0.0452
R² = 0.8074
1
SStot   ( yi  y ) 2
i
yˆ  Ax  b
0.8
0.6
y
0.4
0.2

SS reg   ( yi  y ) 2
0
0
-0.2
y
20
40
60
80
100
120
i
Linearisering av ikke-lineære uttrykk
S (TE )  S o e TE / T 2  ln( S )  ln( S o )  TE / T 2
900
800
700
y = 741.74e-0.019x
R² = 0.9777
500
8
400
7
300
6
200
5
ln(SI)
SI
600
100
0
0
50
100
150
Echo time (ms)
200
250
y = -0.0195x + 6.609
R² = 0.9777
4
3
300
2
1
0
0
50
100
150
Echo time (ms)
200
250
300
Linearisering av ikke-lineære uttrykk
Gamma-variat:
S (t )  A(t  t0 ) e (t T0 ) /  ; t  T0
S (t )  0; t  T0
S (t )  At  e t / 
Dersom T0 kan bli estimert:
45
A=2;alfa=2;beta=2.2
A=2;alfa=2;beta=2.5
A=15;alfa=2.2;beta=2
40
35
30
25
20
15
10
5
0
0
5
10
15
20
25
30
35
40
45
50
Linearisering av ikke-lineære uttrykk
Finne toppunktet,tmax som fn av α,β:
Finne A som fn av S(tmax):
S ' (t )  0    tmax / 
 e

S (t max )  Atmax e   A  S max 
 t max
Uttrykke t som fn av tmax: t’=t/tmax:
S (t ' )  At ' e (1t ')  ln( S (t ' ))  ln( S max )   (1  ln( t ' )  t ' )  ax  b
Dersom T0 og tmax kan estimeres kan gamma-variat funksjonen lineariseres




General Linear Least Squares
n lineære likninger med m ukjente, a1, a2.. am og n>m
m
y( xi )   X j ( xi )a j ; (i  1,2..., n)
j 1
f .eks : y( xi )  a1  a2 xi  a3 xi ; X1 ( xi )  0, X 2 ( xi )  xi , X 3 ( xi )  xi
2
2
Som før, minimere chi-kvadrat:
 y  yˆ ( xi ) 
2   i

i
i 1 

n
2
m
yˆ ( xi )   X j ( xi )a j
j 1
= minimum når deriverte mht alle m parametre, ai , går mot 0:
m

1
0    yi   a j X j ( xi ) X k ( xi ); k  1..m
i 1  i 
j 1

n
Dette kalles normal-likningene, og kan uttrykkes i matriseform som:
( XT X)a  XT y
General Linear Least Squares
Normal-likning:
( XT X)a  XT y
basis-funksjoner: X1(), X2(), ...Xm()
X 2 ( x1 )
.
.
X 2 ( x1 )
2
.
.
.
.
.
.
.
.
.
.
.
1
X M ( x1 ) 

1 
X M ( x2 ) 
2 

.


.
X M ( xN ) 

N 
datapunkter:1..n
 X 1 ( x1 )

 1
 X 1 ( x2 )
 
2
X
.


.
 X 1 ( xN )

 N
 a1 
 
 a2 
a . 
 
 . 
a 
 n
Dersom de m basis-funksjoner er lineært uavhengige er løsningen gitt ved:
a  ( XT X) 1 XT y
(kan løses ved std metoder som Gauss-Jordan elimination)
Se f.eks: Inf-Mat 4350
General Linear Least Squares
Eks: kvadratisk funksjon:
Eks: rett linje:
y  a1  a2 x
1 x1 


1 x2 
X  . . 


. . 
1 x 
n

y  a1  a2 x  a3 x 2
 a1 
a   
 a2 
1 x1

 1 x2
X  . .

. .

1 xn
a  ( XT X) 1 XT y
x12 

2
x2 
. 
. 

xn2 
 a1 
 
a   a2 
a 
 3
General Linear Least Squares
Normal-likning:
( XT X)a  XT y
Dersom de m basis-funksjoner er lineært uavhengige er løsningen gitt ved:
a  ( XT X) 1 XT y
I mange praktiske tilfeller er normal-likningene singulære eller
nær singulære; dvs to eller flere av basisfunskjonene er ikke
lineært uavhengige; => to av funksjonene (eller to forskjellige
kombinasjoner) tilpasser måldata like bra...
Dette betyr at systemet er både overbestemt (flere datapunkter
enn parametre) og underbestemt (ikke uavhengige parametre) på
en gang...
Singulærverdidekomposisjon (SVD)
• For overebestemte systemer gir SVD ’beste’
løsning (i form av minimert minstekvadrat).
• For underbestemte systemer gir SVD den
løsningen som gir de ’minste’ parametere (βi) i
en minstekvadrats forstand: dvs dersom det
finnes kombinasjoner av basisfunksjonene som
er irrelevant mtp løsningen vil denne
kombinasjonen bli redusert til nær null (i
stedet for å bli forstørret opp mot ∞...)
Singulærverdidekomposisjon (SVD)
Ax  y
Dekomponere A til produkt av ortogonale
matriser (UTU=I =>U-1=UT) :
A  UV
T
1
 x  VΣ U y
T
U,  og VT er alltid inverterbare. Diagnonalen i 
er singulærverdiene. Alle ikke-diagonale elementer
i  er null.
SVD
Ax  y
X  UV
 1 / w1
 0

1
Σ  .

 .
 0
.
1 / w2
.
.
.
.
.
.
.
.
T
 A  VΣ 1UT y
.
.
.
1 / wr
.
0 
0 
0 

0 
1 / wM 
w1
wr
rank
wi = systemts singulærverdier
For små singluærverdier setter vi σ=1/w =0 (dvs 1/0 = 0!) siden
disse komponentene er resultat av avrundingsfeil eller støy.
Utfordring: finne riktig cutoff wr
Ikke-lineær kurvetilpasning
•
•
•
•
Minimerings-prosedyre må utføres iterativt
Starte med test-verdier for parametre a1..am
Gjenta interasjonen til χ2(a) er ’liten nok’
Trenger generelt å finne den deriverte av modellfunskjonen
med hensyn på alle modell-parametre
Forutsetning: nær minimum kan χ2(a) approksimeres til en kvadratisk funksjon:
1
2
 2 (a)    d  a  a  D  a
Iterasjons-step:
d=M-vektor; D=MxM matise
D=2. deriverte (Hessian matrise) av χ2(a)
a i 1  a i  k   2 (a i )
en populær algoritme for ikke-lineær kurvetilpasning er Lovenberg-Marquadt metoden.
Ikke-lineær kurvetilpasning
ROI signal intensity
1 300
g
b
c
d
e
f
b
c
d
e
f
g
1 250
1 200
1 150
1 100
1 050
1 000
950
900
850
800
750
700
650
600
550
500
450
400
350
300
250
200
150
100
200
400
600
Inversion Time [ms]
800
ROI #1
T1 relaxation fit
S (TI )  M 0 | 1  e TI / T1 | C
Kurvetilpasning: anvendelser i
perfusjonsanalyse:
Foldingsintegral for perfusjonsanalyse (se del 2):
t
C (t )  f  Ca ( )R(t   )d
0
To mulige metoder for å estimere perfusjon, f:
Parametrisk løsning: forutsetter at residualfunksjonen R(t) kan uttrykkes
parametrisk (dvs som en analytisk, kjent funksjon); f.eks
t
R(t ) 
1
1  t / MTT
R(t )  e
 t / MTT
C (t )  f  Ca ( )e (t  ) d
0
Kan nå løse for perfusjon, f, ved ikke-lineær kurvetilpasning
Ulemper:
• Forutsetter at formen på R(t) er kjent (vanligvis ikke)
• Ikke-lineær kurvetilpasning er data-intensivt og ikke alltid robust (særlig for denne
type uttrykk)
Ikke-parametrisk analyse
Foldings-integral i diskret form:
i
C (ti )  Ft  Ca (t j ) R(ti  t j ); i  1..n
j 1
i matrise-notatsjon:
1
ftR  A C
C  ftAR
0
 Ca(t1 )

 Ca(t 2 ) Ca(t1 )
A .
..

 .
 Ca(t n ) Ca(t n 1 )
0 

..
0 

..
.. 

.. Ca(t1 ) 
..
Arteriell inputfunskjon i konv. matriseform
∆t = samplingintervall
f=perfusjon
 R(t1 ) 


 R(t1 ) 
R  . 


 . 
 R(t n ) 
Diskret residualfunksjon
Ikke-parametrisk analyse
Fordeler med ikke-parametrisk analyse:
•Forutsetter ikke at vi kjenner formen på R(t)
•Kan bruke SVD for å finne f .R
ftAR  C
A  UV
T
1
ftR  A C
ftR  VΣ 1UT C
SVD
FtAR  C
A  UV
Σ 1r
 1 / wmax
 0


.

.

 0
FtR  A 1C
1
FtR  VΣ r U C
T
.
1 / wmax 1
.
.
.
.
.
.
.
.
.
.
.
1 / wr
.
0
0
0
0
0
T
 wmax




 wr

rank
Utfordring: finne korrekt rank for A (# uavhengige kolonner)for
å fjerne støy og beholde mest mulig av sant signal. Bare
beholde de r største singulærverdier
Reelt eksempel
1
FtR  VΣ r UT C
35
250
30
Ca(t)
C(t)
25
200
singlulærverdi, w
20
15
10
Singulærverdier: Σ
150
100
5
0
-5
50
0
5
10
15
20
25
30
35
40
45
50
0
0
AIF of vevsrespons
Σ 1r
5
10
 1 / wmax
 0


.

.

 0
15
20
25
30
Singulærverdi #
.
1 / wmax 1
.
.
.
.
.
.
.
.
35
40
.
.
.
1 / wr
.
45
50
0
0
0
0
0







singulærverdi cutoff
f=175 mL/100 g / min
250
150
100
50
0
0
5
10
15
20
25
30
Singulærverdi #
35
40
45
50
wr= 0.01 wmax
f=128 mL/100 g / min
250
200
singlulærverdi, w
singlulærverdi, w
200
150
100
50
0
0
5
10
15
20
25
30
Singulærverdi #
wr =0.2 wmax
35
40
45
50
singulærverdi cutoff
f=9.6 mL/100 g / min
250
singlulærverdi, w
200
150
100
50
0
0
5
10
15
20
25
30
Singulærverdi #
wr= 0.8 wmax
35
40
45
50
Metoder for å bestemme ’korrekt’ SVD
cutoff:
• Fixed cutoff (f.eks. wr =20% av wmax)
• Iterative estimering:
– øke wr iterativt til oscillering i R(t) er under en gitt
grenseverdi
– Introdusere en dempnings-faktor i 1/w leddene
slik at høye singulærverdier bidrar mer enn lave
(Tikhonov regularisering)
Metoder for å bestemme ’korrekt’ SVD
cutoff:
Iterative estimering: øke σr iterativt til oscillering i R(t) er under en gitt
grenseverdi. Definere en oscillasjons-index (se Wu et al: MRM 50:164–174;2003):
1 1  M 1

OI 
  | R(k )  2 R(k  1)  R(k  2) | 
M  max  k 2

wr=0.02 x wmax
OI=0.32
wr=0.2 x wmax
OI=0.056
wr=0.5 x wmax
OI=0.030
Tikhonov regularisering:
Utgangspunkt for minstekvadratsløsninger er å minimere kvadrated av normen av
forskjellen mellom estimert løsning A.R og målte inputdata C
2
(=’residualnormen’)
R
min || AR  C ||
Liknings-systemet vårt kan (og vil ofte) være ’ill-posed’ og det er problematisk å
få en stabil løsning for R (som fører f.eks til oscillasjoner i R(t)). Vi kan derfor få
en løsning som er ’nøyaktig’ (||AR-C|| er liten) men som ikke nødvendigvis er
korrekt (f.eks ved at R ikke er fysiologisk meningsfull eller oscillerer kraftig)
Tikhonov regularisering: stabiliser systemt ved i stedet å minimere normen for et
alternativt uttrykk:

r  min || AR  C ||  || LR ||
R
2
Matrisen L defineres slik at uttrykket:
2
2 || LR ||2
2

(=’løsningsnormen’)
’straffer’ løsninger som gir en ikke ønskelig løsning for R(t) (f.eks kraftige
oscillasjoner). Valg av L er da basert på apriori forutsetninger om ønskelige
begrensninger i R(t)
Tikhonov regularisering:
I sin enkleste form setter vi L til enhetsmatrisen og λ blir da en enkel
dempnings –faktor som filterer ut små singulærverdier:
 w12
 2
2
 w1  
D=diag. matrise: 
 0
D
.

for små wi: Di ->0

.
for store wi: Di ->1

 0


1
ftR  VDΣ U C
2
i
w
Di  2
wi  2
λ=1
T
λ=50
0
.
w22
w22  2
0
.
.
.
0
λ=200
.
.




0 

.


.

wn2 
wn2  2 
0
Tikhonov regularisering:
Hvoran bestemme riktig verdi av λ ??:
1
1
ftR  A C  VDΣ U C
T
wi2
Di  2
wi  2
Finne en løsning som er en ’optimal’ avveining mellom liten residualnorm og
liten løsningsnorm :
Residualnorm:
Løsningsnorm:
|| ftAR  C ||
|| R ||
mimimert gir ’nøyaktig’ løsning
mimimert gir ’meningsfull’ løsning
Iterativ metode: Øk λ til man når en ’optimal’ tradeoff mellom residualnormen
og løsningsnormen
Tikhonov regularisering:
’L-curve’ optimalisering: Knekkpunktet (pil) kan bestemmes analytisk (P.C. Hansen).
Tikhonov regularisering:
-0.5
-1
-1.5
-2
-2.5
-3
-3.5
-4
-4.5
-0.5
0
0.5
1
1.5
λ=64
Ulemper med iterative Tikhonov metode:
Finne riktig λ (L-kurve knekkpunkt) i reelle data...
Tidkrevende (mange interasjoner nødvendig for å kunne identifisere
knekkpunktet)
2
Kurvetilpasning: anvendelser i
permeabilitetsanalyse
Foldingsintegral for permeabilitetsanalyse (se del 2):
t
C (t )  k1  Ca ( )e k2 (t  ) d  V pCa (t )
0
Tilsvarer parametrisk løsning for perfusjons-foldingsintegralet
D.v.s: kan nå løse for k1 og k2 ved ikke-lineær kurvetilpasning.
Ulemper:
• Ikke-lineær kurvetilpasning er data-intensivt og ikke alltid robust (særlig for denne
type uttrykk)
• Må estimere deriverte av C(t) mhp alle modellparametre (k1,k2,Vp)
Kurvetilpasning: anvendelser i
permeabilitetsanalyse
Kan vi linearisere foldings-integralet??
tar utgangspunkt i transport-likningen (se del 2.2)
d (Ct  V p C p )
dt
 k1C p  k 2 (Ct  V pC p )
Integrerer begge side (og forutsetter Ct(0)=0= gir:
t
t
0
0
Ct (t )  (k1C p  k2V p )  C p ( )d  k2  Ct ( )d  V p C p (t )
I dikret form kan dette uttrykkes som en lineær likning: Ax =C !
SVD: anvendelser i permeabilitetsanalyse
t
t
0
0
Ct (t )  (k1C p  k2V p )  C p ( )d  k2  Ct ( )d  V p C p (t )
=> Ax =C
 t1
 C p ( )d
 0
 t2
 C p ( )d
0

A

.

.
t
N
  C p ( )d
0
t1
  Ct ( )d
0
t2
  Ct ( )d
0
.
.
tN
  Ct ( )d
0

C p (t1 ) 


C p (t 2 ) 

. 

. 

C p (t N ) 

 k1  k 2V p 


x   k2 
 V

p


 Ct (t1 ) 


C
(
t
)
 t 2 
C . 


 . 
 C (t ) 
 t N 
SVD: anvendelser i permeabilitetsanalyse
t
t
0
0
Ct (t )  (k1C p  k2V p )  C p ( )d  k2  Ct ( )d  V p C p (t )
=> Ax =C
kan løses ved SVD!
A  UV
T
 k1  k 2V p 


1
x   k 2   A C  VΣ 1UT C
 V

p


NB: Σ-1 har nå bare tre kolonner (tre singulærverdier) og alle tre kan vanligvis
brukes (ikke behov for å finne noe ’cutoff’ )
Løsning på foldings-integralet ved bruk
av SVD
ROI signal intensity
320
g
b
c
d
e
f
b
c
d
e
f
g
b
c
d
e
f
g
310
300
290
280
270
260
250
240
230
Pixel Intensity [a.u.]
220
210
K1=0.05;
K2=0.157;
Vp=0.2
Ve=0.27
200
190
180
170
160
150
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0
0
50
100
150
200
Time [sec]
250
300
350
ROI #1
ROI #2
Permeability fit