Detection of Change--Points in the Spectral Density. With Applications to ECG Data

2.2 MB·English

by  Pierre R. Bertrand

Detection of Change–Points in the Spectral Density.
With Applications to ECG Data
PierreR.Bertrand,GillesTeyssière,GilBoudet,AlainChamoux
9
INRIASaclay,APISTeam
0
0 UMRCNRS6620,Clermont-FerrandUniversity,France.
2 E-mail:[email protected],
,CREATES,AarhusUniversity,Denmark
n
a E-mail:[email protected]
J ,InstitutdeMédecinedutravail,UFRMédecine,Clermont-FerrandUniversity

0 [email protected]
2 ,InstitutdeMédecinedutravail,UFRMédecine,Clermont-FerrandUniversity

CHUClermont-Ferrand
]
E E-mail:[email protected]
M
. Résumé. Nous proposonsune nouvelleméthode pour estimer les rupturesdu rythme cardiaquedans les
t
a bandes parasympathiques et orthosympathiques basée sur la transformée en ondelettes dans le domaine
t
s complexeetl’étudedesrupturesdesmomentsdesmodulesdecestransforméesenondelettes.Nousobser-
[

vonsdesrupturesdansladistributionpourlesdeuxbandesdefréquence.
1
v
6
9
9
2 1 Introduction
.
1
0
ECGsignalanalysishasalongstoryaftertheimplementationoftheambulatorymonitoringbyHolteratthebeginning
9
0 ofthefifties.Recentmeasurementmethods,duetothesizereductionofthemeasurementdevices,seeChamoux(1984),
:
v allowustorecordECGdataforhealthypeopleoveralongperiodoftime:longdistance(marathon)runners,individuals
i
X daily(24hours)records,etc.Wethenobtainlargedatasetsthatallowustocharacterizethevariationsoftheheartbeatrate
r inthetwocomponentsofthenervousautonomoussystem:theparasympatheticandtheorthosympatheticones.
a
1500
1000
500
50h50 9h43’ 13h38’ 16h32’ 20h07’ 0h14’ 4h52’
FIG. 1– RRintervalforahealthysubjectduringaperiodof24hours
DetectionofChange–PointsintheSpectralDensity
The data studied in this paper have been recorded the 29th of January 2008, at the Clermont-Ferrand University
Hospital (CHU) by G. Boudet. Unlike Boudet et al. (2004) and Cottin et al (2006), these data are not recorded in the
frameworkofalaboratoryexperiment,butduringthereallifeoftheindividualunderinvestigation.
Wehavedailydatasets,over100,000observations,andinthenearfuturewewillhavedataoverseveraldays,i.e.,over
thatsamplesize.WearemeasuringtheintervalbetweentwoRRpeaks,i.e.,ifwedenoteby(ti)i=1,...,N,thesequenceof
peaktimesmeasuredwithaprecisionof1.E-03second,weconsidertheseriesX(ti)= ti−ti−1 measuredinseconds.
Severalindicatorsarerelatedto heartbeatdata.Themostpopularisthe instantaneous(cid:16)averagefr(cid:17)equency,i.e.,X(t )−1,
i
whichisdisplayedbyrunnerswatches.Thisquantityisinformativeondailyactivityofobservedindividual:sleepingand
wakinguptimes, physicalactivity,e.g.,sport,manualwork,butdoesnotsummarizeallrelevantinformation.Notethat
heartbeatdatadisplaylargevariations,clustering,etc,asonlyindividualswithseriousdiseasedisplayaregularheartrate.
Cardiologistsareinterestedinthestudyofthissignalintwofrequencybands:theorthosympatheticandparasympa-
thetic bands,i.e., the frequencybands(0.04Hz,0.15Hz)and (0.15Hz, 0.5Hz)respectively.The definition of these
bandsistheoutcomeofresearchworks,seee.g.,TaskforceoftheEuropeanSocietyofCardiologyandtheNorthAmeri-
canSocietyofPacingandElectrophysiology(1996),andisbasedonthefactthattheenergycontainedinsidethesebands
wouldbearelevantindicatoronthelevelofthestressofanindividual.
Indeed, for the heart rate, the parasympathetic system is often compared to the brake while the orthosympathetic
systemwouldbeaniceaccelerator;seee.g.Goldberger(2001).Atrestthereisapermanentbrakingeffectontheheart
rate. Anysolicitation of the cardiovascularsystem, anyactivity initially producesa reductionof parasympatheticbrake
followedbyagradualinvolvementofthesympatheticsystem.Thesemechanismsareveryinterestingtowatchinmany
diseasesincludingheartfailure,butalsorhythmdisordersthatmayfallunderoneorotherofthesetwoeffects,monitoring
thetherapeuticeffectofseveralMedicinesincludingsomepsychotropic.Inthefieldofphysiologysuchdataarecrucial
formeasuringthelevelofvigilanceandparticularlytheriskoffallingasleepwhiledrivingavehicle,thelevelofstress
inducedbyphysicalactivityorlevelofperceivedstress,whichcanbeconsideredasacriterionofovertraininginsport.
Fractals modelshave been used in cardiologyafter the worksby Ivanovet al. (1999),who appliedthe multifractal
spectrumanalysisadvocatedbyFrisch(1996),formodelingRRseriesandclassifyingindividualsaccordingtothismul-
tifractalspectrum,asthisspectrumdiscriminatesbetweenindividualswhoexperiencedhearthtrouble,andthosewhodid
not. However,this toolhas some shortcomingsas it requireshuge samplescommonin turbulenceanalysis, andis then
inappropriateforstudyingphenomenaoccurringataresolutionlowerthanthedailytimeinterval,suchasthevariations
oftheparasympatheticandorthosympatheticsystemsinsidetheday.
Wavelets basedmethodshavebeen used in biostatistics byDiab etal. (2007)for uterineEMG signalanalysis. Ho-
wever,theyconsiderthat the processstudiedis homogeneous,and used these methodsas a classification tool. Another
significantdifferenceisthefactthattheyusediscretewaveletdecomposition,i.e.,afrequencydecompositiononadyadic
waveletbasis,thechoiceforthefrequencybandsismadewithoutreferencetoabiologicalphenomenon.Inourcase,the
choiceforthefrequencybandisjustifiedbybiologicalconsiderations,andwefitthewaveletsinsidethesebands.
Thisiswhythecontinuousanalysisofbothsystemsandtheirquantificationisaparticularlypromisingresearcharea.
So, the example in Figures 4, 5, 6 shows at the observation 28,220 i.e., 13h40’50”, a simultaneous variation of both
systems.
2 Mathematical formulation
Forbiologicalreasons,theheartrateiswithintheinterval[20,250]bpm(beatsperminute),i.e.,X(t)belongstothe
RRinterval60/250s.< X(t) < 60/20s..Thisleadstomodelingbystationaryorlocallystationaryprocesses.Wewish
to decomposethis signalin a sequenceof homogeneousintervalshavingthe same mean and/orvariance.We makethe
assumptionthattheseriesisGaussian,andthenfindtheoptimalsegmentationoftheprocessusingthesameapproachas
inLavielleandTeyssière(2006)orLavielleandMoulines(2000).
Bertrand,Teyssière,Boudet,Chamoux
We assume thatthe signalis thesum ofa piecewiseconstantfunctionanda Gaussian process,centeredandlocally
stationary.Wethenhavethefollowingrepresentation:
X(t)=µ(t)+ eitξf1/2(t,ξ)dW(ξ), forall t∈R, (1)
Z
IR
where
– ξ 7→f(t,ξ)isanevenandpositivefunction,calledspectraldensitypiecewiseconstant,i.e.,thereexistsapartition
τ1,...,τK suchthatf(t,ξ)=fk(ξ)fort∈[τi,τi+1[
– thefunctiont7→µ(t)isalsopiecewiseconstantforanotherpartitionτ1,...,τL withµ(t)=µℓift∈[τℓ, τℓ+1[
Thewaveletcoefficientassociatedwithψis
e e e e
W (b)= ψ(t−b)X(t)dt unitinsecond2. (2)
ψ
Z
IR
In Bardet and Bertrand (2007) or Bardet et al. (2008), one can find a theoretical study of the wavelet coefficient for
stationary (or with stationary increment) centered Gaussian processes, i.e., for X given by (1) with µ(t) = 0. This
approachcanbegeneralizedtolocallystationaryGaussianprocesses,thiswillbedetailedinaforthcomingpaper.
AccordingtorecommendationsofTaskforceoftheEuropeanSoc.CardiologyandtheNorthAmericanSoc.ofPacing
andElectrophysiology(1996),weusethefollowingnotations:
– [ω1,ω2]=(0.04Hz, 0.15Hz)denotestheorthosympatheticfrequencyband,
– [ω2,ω3]=(0.15Hz, 0.5Hz)denotestheparasympatheticband
Theenergyassociatedwitheachofthesefrequencybandsandlocalizedaroundthetimeb,ismeasuredbythemodulusof
thecomplexwaveletcoefficients|W (b)|2 fori=1,2,with
i
b7−→W (b) = ψ (t−b)X(t)dt, i=1,2.
i i
Z
IR
Inthiswork,thesewaveletscoefficientsarecomputedateachsecond,i.e.,thedifferencebetweentwoconsecutivevalues
forbisequalto1second.
Howtochoosethewaveletsψ1andψ2?
Intheidealisticcase,onewouldusetwofiltersψ1andψ2withcompactsupport,theFouriertransformsofwhichhave
supportinsidetheorthosympatheticandparasympatheticbands.
Unfortunately, it does not exist a non null function ψ with compact time domain support and compact frequency
support, see for instance Mallat (1998)Th 2.6 p.34. Therefore,the best we can do is to choose between a filter with a
compactfrequencysupportandafilterwithacompacttimedomainsupport.Thefirstchoiceiswellsuitedforstationary
models, see Bardet and Bertrand (2007). But, in this work, we are interested by locally stationary models, thus our
specificationsareafilterwithcompacttimedomainsupportasaDaubechieswavelet.Thepricetopayforthecompactness
of the time domainsupportis the loss in the compactnessof the frequencysupport.However,the frequencysupportis
"almostcompact"inthefollowingsense:
Definition2.1(ρpseudosupport) Let0<ρ<1,beamapg ∈L2(IR)thatadmitsthecompactintervalIasaρpseudo
|g(t)|2dt
supportif I =ρ.
R |g(t)|2dt
IR
R
FouriertransformofDaubechieswavelethaveareasonablysmallρpseudosupportwitharatioρcloseto1.Moreover,the
largerthenumberoftheDaubechieswaveletis,thecloserto1theratioρis;seetheexamplebelowwiththeDaubechies
waveletD6.
By scaling and modulation, one can adjust the pseudo support inside a specified frequency band as stated by the
followingproposition
DetectionofChange–PointsintheSpectralDensity
Proposition2.1 Let ψ be a filter with compact support [L1,L2] and a frequency ρ pseudo support [Λ1,Λ2], for any
frequenciesband[ω1,ω2],themap ψ1(t)=µ×eiηtψ(λt)with
ω2−ω1 ω1+ω2 Λ2+Λ1
µ>0, λ= and η = −(ω2−ω1)
Λ2−Λ1 2 Λ2−Λ1
Λ2−Λ1 Λ2−Λ1
hasaρpseudosupport[ω1,ω2]andatimedomainsupport L1, L2 .
(cid:20)ω2−ω1 ω2−ω1 (cid:21)
ξ−η
Proof. Sinceψ1(ξ) = µ×ψ ,onecandeduce ρ pseudosuppψ1 = η+λ×ρ pseudosuppψ andthenthe
(cid:18) λ (cid:19)
proposition. c b
Differentchoicesforthewaveletsψ1andψ2
FromProposition2.1,onecandeducethedifferentpossiblechoicesofthefiltersψ1andψ2
– DaubechieswaveletD6:Inthiscase,wehaveΛ1 =0.08,Λ2 =1.75,ρ=0.9999andonecanset
ψ1(t)=µ×eiη1tD6(λ1t) and ψ2(t)=µ×eiη2tD6(λ2t)
with η1 =−0.0255,λ1=0.0659,η2=−0.0585, andλ2 =0.2096.
Thelengthofthetimesupportare|Suppψ1|and|Suppψ2|.OnecanseeonFig.2(a)thattheFouriertransforms
ψ1(x)andψ2(x)havealmostdisjointsupports.
b b
ω2−ω1 ω1+ω2 Λ2+Λ1
λ1 = and η1 = −(ω2−ω1)
Λ2−Λ1 2 Λ2−Λ1
ω3−ω2 ω2+ω3 Λ2+Λ1
λ2 = and η2 = −(ω3−ω2)
Λ2−Λ1 2 Λ2−Λ1
Λ2−Λ1 Λ2−Λ1
|Suppψ1|= ×|SuppD6| and |Suppψ2|= ×|SuppD6|
ω2−ω1 ω3−ω2
– Gaborwavelet:Forcomputationalreasons,wewillusetheGaborwaveletdefinedas
ψ(t)=eiηtg(t), g(t)= 1 e−2tσ22 (3)
(σ2π)1/4
see,e.g.,Mallat(1998).Thiswavelethasthesametimeandfrequencyρpseudosupport[−L,L]=[−3.5,3.5]with
ρ=0.9995.Inthespectraldomain,wehave
ψ(t)=gˆ(ξ−η), gˆ(ξ)=(4πσ2)1/4e−σ22ξ2 (4)
b
WefittheGaborwaveletinsidetheorthosympatheticortheparasympatheticfrequencybandsbyusingProp.2.1or
directcalculations.OneobtainsthetwoGaborwaveletsdefinedby(3)withthefollowingchoicefortheparameters:
ω1+ω2 2L
η1 = and σ1 =
2 ω2−ω1
ω2+ω3 2L
η2 = and σ2 =
2 ω3−ω2
4L2 4L2
moreover |ρpseudoSuppψ1|= and |ρpseudoSuppψ2|= withρ=0.9995
ω2−ω1 ω3−ω2
OnecanseeonFig.2thattheFouriertransformsψ1(x)andψ2(x)stillhavealmostdisjointsupports.
UsingtheGaborwaveletismoreefficientintermsofcomputingtime,asitisatleast8timesfaster.Figure3below
b b
displaystheGaborwaveletscoefficientsintheorthosympatheticandparasympatheticbandsrespectivelyforthesample
plottedinFigure1.
Bertrand,Teyssière,Boudet,Chamoux
12
ψ1
ψc2
c
10
8
6
4
2
−00 ,2 −0,1 −0 0,04 0,1 0.15 0,2 0,3 0,4 0,5 0,6 0,7
FIG. 2– TheFouriertransformsψ1(x)(left)andψ2(x)(right)
b b
1000 600
400
500
art art
x p x p 200
e e
pl 0 pl
m m
o o 0
c c
1 2
W −500 W
−200
−1000 −400
1000 500
500 10 10
8 8
0 6 0 6
−500 4 x 104 4 x 104
2 2
W real part −1000 0 W real part −500 0
1 time 2 time
FIG.3 – Thewaveletcoefficientsintheorthosympatheticband(left)andintheparasympatheticband(right)ofthesame
healthysubjectduringaperiodof24hours
3 Segmentation analysis
Weassumethattheprocess{X }isabruptlychangingandischaracterizedbyaparameterθ ∈Θthatremainsconstant
t
between two changes. We use this assumptionto define our contrastfunction J(τ,X). Let K be some integerand let
τ ={τ1,τ2,...,τK−1}beanorderedsequenceofintegerssatisfying0<τ1 <τ2 <...<τK−1 <n.
Forthedetectionofchangesinthemeanandvarianceofasequenceofrandomvariables,i.e.,achangeindistribution,
DetectionofChange–PointsintheSpectralDensity
weusethefollowingcontrastfunction,basedonaGaussianlog–likelihoodfunction:
1 K kX −X¯ k2
J (τ,X)= τk τk +n log(σˆ2). (5)
n n σˆ2 k k
Xk=1 k
where nk = τk − τk−1 is the length of segment k, σˆk2 is the empirical variance computed on that segment k, σˆk2 =
n−k1 τi=kτk−1+1(Xi−X¯)2,andX¯ istheempiricalmeanofX1,...,Xn.
InPfact,weminimizeapenalizedversionofthiscontrastfunction,thepenaltyparameterbeingselectedinanadaptive
waysothattheobtainedsegmentationdoesnotdependtoomuchonthepenaltyparameter;seealsoBirgéandMassart
(2007)for another choice for the penalty parameter. Althoughthis method has been devised for Gaussian processes, it
empirically providesrelevant results for non-Gaussian processes, e.g., financial time series, see Lavielle and Teyssière
(2006)forfurtherdetails.
3.1 Orthosympathetic band
Forthatfrequencyband,weobtainthefollowingsegmentation:
τ ={28220,33366,71048},
i.e.,13h40’50”,15h06’36”,1h34’38”.
500
450
400
350
300
250
200
150
100
50
0
5h50’30" 13h40’50" 15h06’36" 1h34’38"
FIG. 4– Segmentationinthemeanandvariance(changeinthedistribution)ofthemodulusofthewaveletcoefficientsin
theorthosympatheticband
3.2 Parasympatheticband
Forthatfrequencyband,weobtainthefollowingsegmentation:
τ ={11620,21912,28054,31540,36022,40172,52622,70550},
i.e.,9h04’10”,11h45’42”,13h38’04”,14h36’10”,15h50’22”,17h00’02”,20h27’32”,1h26’40”.
4 Conclusion
TheexampleillustratedbyFigures4,5and6showsataroundt = 28,220asimultaneousvariationofbothsystems.
But,onecanalsoobservechange-pointsintheorthosympatheticandparasympatheticbandsoccurringatdifferenttimes.
Bertrand,Teyssière,Boudet,Chamoux
400
350
300
250
200
150
100
50
0
5h50’30" 9h04’10" 11h45’42" 13h38’14h36’ 15h50’17h00’02" 20h27’32" 1h26’40"
FIG. 5– Segmentationinthemeanandvariance(changeinthedistribution)ofthemodulusofthewaveletcoefficientsin
theparasympatheticband
Inthefuture,wewillstudytheexistenceofapossiblecausalityorsequentialitybetweenthesechange–pointsindifferent
bands.
900
800
Orthosympathetic band
700 Parasympathetic band
Hearthbeat process
600
500
400
300
200
100
5h050 ’30" 9h04’10" 11h45’42" 13h38’14h36’ 15h50’17h00’02" 20h27’32" 1h26’40"
FIG.6– Meanofthemodulusifthewaveletscoefficientsfortheorthosympatheticandparasympatheticbands,andofthe
RRprocessX(t)
Références
[1] Bardet,J.M.andBertrand,P.R.(2007).IdentificationofthemultiscalefractionalBrownianmotionwithbiomecha-
nicalapplications.J.TimeSer.Anal.28,1–52.
[2] Bardet,J.M.,Bertrand,P.R.andBillat,V.(2008).Estimationnon-paramétriquedeladensitéspectraled’unprocessus
gaussienobservéàdesinstantsaléatoires.Pub.ISUP52,123–138.
[3] Birgé, L. and Massart, P. (2007).Minimalpenalties for Gaussian modelselection. Probab.Theory Related Fields
138,33–73.
DetectionofChange–PointsintheSpectralDensity
[4] Boudet,G.,Albuisson,E.,Bedu,M.andChamoux,A.(2004).Heartraterunningspeedrelationshipsduringexhaus-
tiveboutsinthelaboratory.Can.J.Appl.Physiol.29(6),731–742.
[5] Chamoux,A.(1984).LesystèmeHolterenpratique.MédecineduSport58,43-273,54–284.
[6] Cottin,F.,Leprêtre,PM.,Lopes,P.,Papelier,Y.,Médigue,C.andBillatV.(2006).Assessmentofventilatorythre-
sholdsfromheartratevariabilityinwell-trainedsubjectsduringcycling,Int.J.SportsMed.27,959–967.
[7] Diab,M.O.,Marque,C.andKhalil,M.A.(2007).ClassificationforuterineEMGSignals:ComparisonbetweenAR
modelandstatisticalclassificationmethod.Int.J.ComputationalCognition5,8–14,
[8] Frisch,U.(1996).Turbulence.CambridgeUniversityPress.
[9] Goldberger,A.L.(2001).Heartbeats,hormonesandhealth:isvariabilitythespiceoflife?Am.J.Crit. CareMed.
163,1289–1290.
[10] Ivanov,P.C., Amaral, L.A.N., Goldberger,A.L., Havlin, S., Rosenblum,M., Struzik,Z, and Stanley,H.E. (1999).
Multifractalityinhumanheartbeatdynamics.Nature399,461–465.
[11] Lavielle,M.andMoulines,E.(2000).Least-squaresestimationofanunknownnumberofshiftsinatimeseries.J.
ofTimeSeriesAnal.21,33–59.
[12] Lavielle, M. andTeyssière,G. (2006).Detection ofmultiple change-pointsin multivariatetime series. Lithuanian
Math. J. 46, 287–306. [Détection de ruptures multiples dans des séries temporelles multivariées.] (Translation in
French)Liet.Mat.Rink.46,351–376.
[13] Mallat,S.(1998).Awavelettourofsignalprocessing.AcademicPress.
[14] Task force of the European Soc. Cardiology and the North American Society of Pacing and Electrophysiology
(1996).Heartratevariability.Standardsofmeasurement,physiologicalinterpretation,andclinicaluse.Circulation
93,1043–1065.
Summary
Weproposeanewmethodforestimatingthechange-pointsofheartrateintheorthosympatheticandparasympathetic
bands,basedonthewavelettransforminthecomplexdomainandthestudyofthechange-pointsinthemomentsofthe
modulusofthesewavelettransforms.Weobservechange-pointsinthedistributionforbothbands.