Employing stochastic differential equations to model wildlife motion

24 Pages·2003·0.31 MB·English

BULLETIN
BRAZILIAN
MATHEMATICAL
SOCIETY BullBrazMathSoc,NewSeries33(3),385-408
©2002,SociedadeBrasileiradeMatemática
Employing stochastic differential equations to
model wildlife motion
David R. Brillinger, Haiganoush K. Preisler,
AlanA.Ager, John G. Kie and Brent S. Stewart
Abstract. Theconcerniswiththepropertiesofstochasticdifferentialequations(SDEs)
describingthemotionofparticlesin3dimensionalspace,onthesphereorintheplane.
Thereisconsiderationofthecasewherethedriftfunctioncomesfromapotentialfunc-
tion. ThereisstudyofSDEswhoseparametersareperiodicintime. Theseareuseful
forincorporatingcircadianrhythminthebehavior. Thecasesofasealinafrozenlake
inAlaska,anelephantsealmigratingagreatdistanceinthePacificOceanandofagroup
of“free-ranging”elkinareserveinOregonarereferredto. Fortheelknonparametric
estimatesofthedriftandvariancetermsofanSDEmodelarediscussedandevaluated
and the fit of the model assessed. One issue is how to include explanatories, beyond
locationandtime,inthemodel.Anumberofquestionsmotivatedbythewildlifemotion
concerningdiffusionprocessesofthetypeconsideredareposedattheendofthepaper.
Keywords: Circadianrhythm,Diffusionmodel,Elk,Elephantseal,Nonparametricre-
gression,Potentialfunction,Ringed-seal,Stochasticdifferentialequation,Vectorfield.
Mathematicalsubjectclassification: 60J60,62G08,62M10,70F99.
1 Introduction
Theconcernistheuseofvector-valuedstochasticdifferentialequations(SDEs)
to describe the motion of wildlife, with the particular cases of seals and elk
focused on. To illustrate the character of the data Figure 1 presents the 3D
trajectory of a ringed-seal in a frozen lake inAlaska, while Figure 2 shows the
pathofanelephantsealmigratingoutintoandthenbackfromtheNorthPacific
Ocean,andlastlyFigure3showsanelkmovingaroundinareserveinOregon.
The trajectories of the animals are sampled about every 30 seconds, every 24
hoursandevery2hoursrespectively.
Received25April2002.
386 D.R.BRILLINGER, H.K.PREISLER, A.A.AGER, J.G.KIE AND B.S.STEWART
Z-120-100-80-60-40-20 010080Y604020 0 0 20 (cid:127)40 60 X80 100 120 140 Z00000002468 1-----021- 020406X08010(cid:127)0120140 0 20 40 6Y0 80 100
(cid:127) (cid:127)
Z0 00000020864211------ 0204Y06080100140120100 80 X60 40 20 0 Z0 00000002246811------14012010080X60 40 20 0 100806040Y20 0
Figure 1: The trajectory of a ringed-seal diving into a frozen lake dot as seen
from 4 viewpoints. The seal was released at the surface point indicated by the
dot.
InthecaseofthesealinAlaskathelakewasfrozenandtheconcernwashow
didthesealfinditswaytoairholes. Thedotinthefigureindicatestheicehole
intowhichtheanimalwasreleased. Thesealisseentoswimtothebottomfairly
directly,toswimaroundthereforawhileandthentograduallycomebacktothe
originalhole. Inthecaseoftheelephantsealsonewondershowtheynavigate?
Where do they go? Is their speed approximately constant? Figure 2 indicates
thatsometimestheyfollowclosetogreatcircleroutes. Inthecaseoftheelk,the
reservemanagerswereconcernedwithquestionslike: Howtoallocateforage?
Ischangetakingplace? Whataretheeffectsoftraffic? Whatisthesequenceof
habitatuse? Onecanconjecturewhatthingsareimportanttotheanimals: The
locationofcover? Thetypeofforage? Thepresenceofroads? ...
An important aspect of the use of the SDEs to model such paths is that the
drift and variance terms are meant to include phenomena such as: attraction,
repulsion, barriers, time of day, ... There are interesting statistical questions:
Howtoincludeexplanatoryvariables? Howtoassessthefitofthemodel? How
to predict future motion? ... A number of interesting questions specifically
addressedtoprobabilistsarelistedinthelastsectionofthepaper.
BullBrazMathSoc,Vol.33,N.3,2002
STOCHASTIC DIFFERENTIAL EQUATIONS TO MODEL WILDLIFE MOTION 387
45 (cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)
(cid:127) (cid:127) (cid:127)
latitude (cid:127)(cid:127) (cid:127) (cid:127)(cid:127) (cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127) (cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127) (cid:127)(cid:127) (cid:127)(cid:127)
(cid:127) (cid:127)
40 (cid:127) (cid:127)(cid:127) (cid:127)
(cid:127)(cid:127)(cid:127)(cid:127) (cid:127)(cid:127)
(cid:127)
(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)(cid:127)
(cid:127)
(cid:127)
(cid:127)
(cid:127)
(cid:127)
35 (cid:127)(cid:127)(cid:127)(cid:127) (cid:127)(cid:127)(cid:127)
(cid:127) (cid:127)(cid:127)(cid:127)(cid:127)(cid:127)
(cid:127) (cid:127)(cid:127)
-160 -150 -140 -130 -120
longitude
Figure 2: Points along the trajectory of an elephant seal. The curved line is a
greatcircleroutepresentedforcomparison. Thesealstartsitsjourneyfroman
islandofSantaBarbara,California.
Thepaperbeginswithadescriptionofdeterministicandstochasticmethodsfor
describingthepathsfollowedbyparticlesundertheinfluenceofafield. Section
3providestwospecificmodels. Section4presentssomedetailsofthestatistical
methods employed. In Section 5 the experiment in which the elk data were
collected is described in some detail. Section 6 presents the results obtained.
Section 7 is Discussion and Summary. The final section poses some specific
questionsfordiffusionprocessspecialists.
Referencespresentingmodelsforanimalmovementinclude: DunnandGip-
son [12], Dunn and Brisbin [11], Preisler and Akers [23], White and Garrott
[33],BrillingerandStewart[7],Turchin[32]. ThereferencesWhiteandGarrott
[33],Turchin[32])setdownFokker-Planckdeterministicdifferentialequations
(DDEs)fordensityfunctionsdescribingtheexpectedpatternofspaceuse. Like-
lihoods may be set down using conditional densities and may be used to make
inferences. Thisdeterministicapproachistobecontrastedwiththatinthepresent
BullBrazMathSoc,Vol.33,N.3,2002
388 D.R.BRILLINGER, H.K.PREISLER, A.A.AGER, J.G.KIE AND B.S.STEWART
14
12
10
8
m
K
6
4
2
0
0 2 4 6 8
Km
Figure3: PointsalongthetrajectoryofoneoftheelkmovingwithintheStarkey
ExperimentalForest.
BullBrazMathSoc,Vol.33,N.3,2002
STOCHASTIC DIFFERENTIAL EQUATIONS TO MODEL WILDLIFE MOTION 389
and other papers where stochastic equations are set down describing the indi-
vidualpaths. Fokker-Planckequationsofdesiredordermaybederivedfromthe
SDEmodels. TheSDEapproachisappropriateforleadingdirectlytoresiduals,
simulationandlikelihood.
2 Approachestothedescriptionofmovingparticles
The analytic formulation of the motion of particles is a traditional problem of
physics and applied mathematics. Classical approaches have been developed
offeringpropertiesofsolutionsoftheequationsofmotionandinterpretationsof
theparametersinvolved. Thissubjectmatterisusefulformotivatingtheresults
ofthepresentwork. Firstconsiderthedeterministicapproach.
2.1 Deterministiccase
Motion in Newtonian dynamics may be described by a potential function,
H(r,t), see Nelson [21]. Here r = (x,y)τ is location and t is time. The
equationsofmotiontaketheform
dr(t) = v(t)dt (2.1)
dv(t) = −βv(t)dt −β∇H(r(t),t)dt
withr(t)theparticle’slocationattimet,v(t)theparticle’svelocityand−β∇H
theexternalforcefieldactingontheparticle,β beingthecoefficientoffriction.
Here ∇ = (∂/∂x,∂/∂y)τ is the gradient operator. The function H is seen to
control the particle’s direction and velocity. For example H(r) = |r − a|2
corresponds to motion with a point of attraction at a and H(r) = 1/|r − a|2
correspondstomotionwithapointofrepulsionata.
Inthecasethattherelaxationtimeβ−1issmall(frictionishigh),theequations
areapproximately
dr(t) = −∇H(r(t),t)dt (2.2)
andthevelocity,v(t),isnolongerinvolveddirectly,seeNelson[21].
Therehasbeenconsiderablemathematicaldevelopmentofthismateriale.g.,
Goldstein [14]. An interesting question given a force field, F, is whether there
exists a real function H, such that F = ∇H? When it does exist, the field is
calledconservative,seeStewart(1991)[31]. Thisquestionisaddressedforthe
StarkeyelkinBrillingeretal. [6].
BullBrazMathSoc,Vol.33,N.3,2002
390 D.R.BRILLINGER, H.K.PREISLER, A.A.AGER, J.G.KIE AND B.S.STEWART
2.2 StochasticCase
Apertinentprobabilisticconceptfordynamicmotionofaparticleisastochastic
differentialequation(SDE),e.g.,seeNelson[21],KarlinandTaylor[17],Bhat-
tacharyaandWaymire[3]. SuchequationsoftenleadtoMarkovprocessesand
takethefollowingform
dr(t) = µµ(r(t),t)dt+(cid:4)(cid:4)(r(t),t)dB(t) (2.3)
Herer,µµ,Barevectorswhile(cid:4)(cid:4)isamatrix. Thedriftµµmaybeinterpretedasa
velocityfieldandanexampleofanestimatewillbeprovidedlaterinthepaper.
(cid:4)(cid:4) is the variance or diffusion parameter while B is a random function such as
aBrownianorLevyprocess. TheparametersandtheBrownianprocesscontrol
thedirectionandspeedofmotion.
Theparametershaveinterpretationsprovidedby
E{dr(t)|Ht} = µµ(r(t),t)dt
and
var{dr(t)|Ht} =(cid:4)(cid:4)(r(t),t)dt
withdtissmallandHtrepresentingthetimehistoryoftheprocess,{r(u),u ≤ t}.
The driving process B leads to variability around deterministic motion. This
process might correspond to explanatories omitted from the equations. The
vectorµµ(r(t),t)isseentorepresenttheinstantaneousvelocityoftheparticleat
time t and position r. Since the process is Markov, these conditional moments
dependonlyontheprecedingposition,r(t).
AvarietyofpropertiesareknownconcerningsolutionsofSDEs,forexample
whenH doesnotdependont and(cid:4)(cid:4)(cid:4)(cid:4)τ = σ2I,thereisoftenaninvariantdensity
0
π(r) = cexp{−2H(r)/σ2} (2.4)
0
representing the longrun density of locations that the process visits; see Bhat-
tacharyaandWaymire[3]. Thusinsuchastationarycase,byanalyzingthepaths,
populationdensitiesmaybeestimated.
The situation of a particle being affected by the force field of a potential
function may be visualized by picturing a ball rolling around in the interior
of a perspective plot of the potential function. In the stochastic case there is
a convenient approximation for viewing the situation. First consider a one-
dimensionalprocessdx = µ(x,t)dt+σ(x,t)dB(t). Supposethatattimet the
BullBrazMathSoc,Vol.33,N.3,2002
STOCHASTIC DIFFERENTIAL EQUATIONS TO MODEL WILDLIFE MOTION 391
particle is at location x(t) = x. Then for the particle’s location at time t +dt
take
√ 1 µ(x,t) √
x(t +dt) = x ±σ(x,t) dt withprob ± dt
2 2σ(x,t)
seeKloedenandPlaten[19];ProhorovandRozanov[26]. Asdt → 0thepaths
ofthisprocessapproachthoseoftheSDEwithparametersµ,σ. Inthebivariate
case one generates two such processes. Examples are given in Brillinger et al.
[6].
A random potential/environment may be a useful model. It might be used to
describesealschasingfish,elkattractedorrepelledbyotherelk. Therandomness
is introduced by the motion of these other particles. For the potential function
onemightwrite
(cid:1)J
H(r,t) = α(r,t) |r−Xj(t)|2
1
withtheXj(t)locationsofmovingattractors/repellors. Thiscouldbethesource
ofthedB(t)termin(2.3).
The components Hx,Hy of the gradient ∇H correspond to the components
ofµµof(2.3). Oneadvantageofthepotentialfunctionapproachisthatindepen-
dentpotentialfunctionsfromavarietyofsourcesaddasin
(cid:1)
Hl(r,t)
l
There may be points, lines or regions of attraction or repulsion and barriers
tobeincluded. Thebarrierscanrepresentactualphysicalobjects(e.g.,fences).
The process B includes impulses due to the presence of objects such as trees,
mounds,dipsandnaturalvariabilitycorrespondingtoindividualanimalbehavior.
The fitted SDEs may be used to produce estimates of other parameters, (e.g.
expectedspeedanddistancestoroad),topredictspatialandtemporalpatternsof
animaldistributionandhabitatpreferences,tosimulatetrajectoriesandtostudy
thedirectionalityofthemovementforexample.
3 Somespecificmodels
3.1 Ornstein–Uhlenbeck
A particular case of an SDE that has been employed in describing animal mo-
tion is the mean-reverting Ornstein–Uhlenbeck (O–U) process; see Dunn and
BullBrazMathSoc,Vol.33,N.3,2002
392 D.R.BRILLINGER, H.K.PREISLER, A.A.AGER, J.G.KIE AND B.S.STEWART
Gipson[12];DunnandBrisbin[11]. Here
µµ(r,t) = A(a−r)
and
(cid:4)(cid:4)(r,t) =(cid:4)(cid:4)
whilethemeanisa. TheO–UprocessbecomestherandomwalkwhenA = 0,
i.e.,whenthedriftterm,µµ(r,t),is0.
If A is symmetric and positive definite, the corresponding potential function
is
H(r,t) = (a−r)τA(a−r)/2
and the particle is wandering but being pulled towards the location a. The
invariantdistributionismultivariatenormal,N(a,(cid:8)),with
(cid:2)
∞
(cid:8) = e−Au(cid:4)(cid:4)(cid:4)(cid:4)τe−Audu;
0
seep.597inBhattacharyaandWaymire[3]. If(cid:4)(cid:4)(cid:4)(cid:4)τ = σ2I,then(cid:8) = σ2A−1/2.
0 0
3.2 Diffusiononasphere
Inthecaseoftheelephantsealsonecanconsiderthecaseofadiffusionprocess
with drift on the sphere, see Brillinger [4]. With φ(t) denoting latitude and
θ(t)denotingcolatitudetheequationsoftheprocessare
σ2
dθ = σdU +( −δ)dt
t t 2tanθt
σ
dφ = dV
t sinθt t
Hereδ isthedriftparameterand(Ut,Vt)denotesabivariateBrownianmotion.
4 Thestatisticalmethodsused
One can consider both parametric and nonparametric methods of inference.
DunnandGipson[12]useapproximatemaximumlikelihoodprocedurestoesti-
matetheparametersa,(cid:8)ofthemultivariateO–Uprocessfromsampledtrajecto-
rieswithconstantsamplingintervals. DunnandBrisbane[11]giveextensionsof
themaximumlikelihoodestimatetothecasewhereobservationsareunequally
spacedovertime.
BullBrazMathSoc,Vol.33,N.3,2002
STOCHASTIC DIFFERENTIAL EQUATIONS TO MODEL WILDLIFE MOTION 393
Turning to a more general situation, the Radon-Nikodyn derivative of the
process(2.3)withrespecttothemeasureofr(t) = σ(r,t)dB(t)is
(cid:3) (cid:2) (cid:2) (cid:4)
T T
exp − σ(r,t)−1µ(r,t)·σ(r,t)−1dr+ |σ(r,t)−1µ(r,t)|2dr/2
0 0
[22],andestimatesofafinitedimensionalparametermaybecomputedbymax-
imization. In the sampled time case the SDE (2.3) may be approximated as
follows
√
(r(tl+1)−r(tl))/(tl+1−tl) ≈ µµ(r(tl),tl)+(cid:4)(cid:4)(r(tl),tl)Zl/ tl+1−tl (4.1)
l = 1,2,... witht < t < t < ··· samplingtimesandwiththeZ independent
1 2 3 l
bivariatestandardnormals. Thisfollowsfrom(2.3)directly. Thevalidityofthe
approximationisinvestigatedinKloedenandPlaten[19].
Intermsoftheindividualcomponents(X,Y)ofronehasfrom(4.1)
(cid:12)X
(cid:12)t = µx(X,Y,t)+noise
(cid:12)Y
(cid:12)t = µy(X,Y,t)+noise (4.2)
Ifthedriftfunctioncomponents,µx,µy,aresmooth,oneseesthatonehasanon-
parametricregressionestimationproblem. Principalstatisticaltoolsemployedin
theworkaresmoothingmethodsandresidualplots. Thesmoothingornonpara-
metricregressionprocedureemployedtoestimateµx,µy isthefunctionloess()
of Cleveland et al. [9] within gam() of Hastie [15]. It involves the local fitting
ofquadraticsintheexplanatories.
Intheexampleoftheelk,timeofdayprovesimportantandconditionaldensity
estimatesareevaluatedforselectedtimesofday. Estimateswillbeprovidedof
conditionallongrunpopulationdensitiesandoffunctionsdependingsmoothlyon
timeandlocationandexplanatories. Thereareavarietyofsuchestimatesinclud-
ing: kernel-based,spline-based,localpolynomial,HastieandTibshirani[16].
References to inferential methods for diffusion processes, both parametric
and nonparametric, include Banon and Nguyen [1], Bertrand [2], Burgière [8],
Dohnal[10],Genot-Catalotetal. [13],Sorensen[29],PrakasaRao[22].
Aquestionishowtoincludeexplanatoryvariablescorrespondingtophenom-
enasuchas: fences,roads,streams,forage,cover,terrain,timeofday,timesof
sunrise and sunset, periods of hunting. Attraction/repulson phenomena might
be described by a potential term H(r,t) = αd(r)β where d(r) is the shortest
distancetoafeature.
BullBrazMathSoc,Vol.33,N.3,2002
394 D.R.BRILLINGER, H.K.PREISLER, A.A.AGER, J.G.KIE AND B.S.STEWART
5 TheElkexperiment
ThemainstudyareaatStarkeyExperimentalForestandRangeis7,762halarge
andlocatedintheBlueMountainsofnortheasternOregon,Rowlandetal. [27].
The forest was enclosed with a 2.4-m tall, net-wire fence in 1988 and radio-
telemetrystudieswerebegun. Eachspringanumberofelk,deer,andcattleare
fitted with collars containing Loran-C receivers. The collars are instructed to
interceptLoran-Cbroadcastsatregularintervalsandthentorelaythosesignals
toacentralreceiver. Locationsarethenestimatedfromthetimedelays.
Thestudyareaisalsomanagedforavarietyofpublicusessuchasrecreation,
hunting, forest management, cattle grazing, and other activities. The shape of
theareaisshowninFigure3. Thewhiteareasinthefigurescorrespondtotwo
smallelk-proofexclosureswithinthestudyarea.
The data are spatial-temporal. The locations of M = 53 elk, (labelled by
m = 1,... ,M, and recorded at times, tmk,k = 1,... ,Km for the m-th an-
imal) are given as well as various explanatory variables describing vegetation
andtopography. Otherhabitatfeatures(e.g.,distancetoroad,distancetowater)
suspectedtoinfluenceelkmovement,arealsoavailable. Thelocationsarewrit-
ten as a column vector rmk = (Xm(tmk),Ym(tmk))τ, corresponding to the UTM
(Universal Transverse Mercator) coordinates of the k-th time measurement of
them-thelk.
The data used in the work of this paper were collected betweenApril 7 and
November15,1994andinvolve53femaleelk. Observationswereomittedfor30
daysduringtheautumnof1994whenhuntingwasconductedwithintheproject
area. Preliminary analyses of these data often showed erratic movements that
werenottypicaloftherestoftheyear. Alsoomittedweremovementswherethe
timeintervalexceeded1.5hoursorwaslessthan0.1hours. Velocitiescalculated
from short time intervals are strongly influenced by telemetry error while long
timeintervalsgenerateuncertaintyregardingthetruetrajectoriesofelk.
Figure 3 showed the successive estimated locations for one of the elk. The
trajectories plotted are a sequence of straight-line segments and therefore are
jagged. This discreteness results from the fact that location estimates are only
available once every 0.1-4 hours. Figure 5, to follow, will show the estimated
distributionofall53animalsatfourdifferenttimesoftheday.
BullBrazMathSoc,Vol.33,N.3,2002