\documentclass{article} % \VignetteIndexEntry{adehabitatLT: Analysis of Animal Movements} \usepackage{graphicx} \usepackage[colorlinks=true,urlcolor=blue]{hyperref} \usepackage{color} \usepackage{url} \usepackage{Sweave} \newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}} \let\pkg=\strong \DeclareTextSymbol{\degre}{T1}{6} \title{ Analysis of Animal Movements in R:\\ the {\tt adehabitatLT} Package } \author{Clement Calenge,\\ Office national de la chasse et de la faune sauvage\\ Saint Benoist -- 78610 Auffargis -- France.} \date{March 2019} \begin{document} \maketitle \tableofcontents <>= owidth <- getOption("width") options("width"=80) ow <- getOption("warn") options("warn"=-1) .PngNo <- 0 wi <- 600 pt <- 15 @ <>= .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".png", sep="") png(file=file, width = wi, height = wi, pointsize = pt) @ <>= dev.null <- dev.off() cat("\\includegraphics[height=10cm,width=10cm]{", file, "}\n\n", sep="") @ <>= dev.null <- dev.off() cat("\\includegraphics[height=14cm,width=14cm]{", file, "}\n\n", sep="") @ \section{History of the package adehabitatLT} The package {\tt adehabitatLT} contains functions dealing with the analysis of animal movements that were originally available in the package {\tt adehabitat} (Calenge, 2006). The data used for such analysis are generally relocation data collected on animals monitored using VHF or GPS collars.\\ I developped the package {\tt adehabitat} during my PhD (Calenge, 2005) to make easier the analysis of habitat selection by animals. The package {\tt adehabitat} was designed to extend the capabilities of the package {\tt ade4} concerning studies of habitat selection by wildlife.\\ Since its first submission to CRAN in September 2004, a lot of work has been done on the management and analysis of spatial data in R, and especially with the release of the package {\tt sp} (Pebesma and Bivand, 2005). The package {\tt sp} provides classes of data that are really useful to deal with spatial data...\\ In addition, with the increase of both the number (more than 250 functions in Oct. 2008) and the diversity of the functions in the package \texttt{adehabitat}, it soon became apparent that a reshaping of the package was needed, to make its content clearer to the users. I decided to ``split'' the package {\tt adehabitat} into four packages:\\ \begin{itemize} \item {\tt adehabitatHR} package provides classes and methods for dealing with home range analysis in R. \item {\tt adehabitatHS} package provides classes and methods for dealing with habitat selection analysis in R. \item {\tt adehabitatLT} package provides classes and methods for dealing with animals trajectory analysis in R. \item {\tt adehabitatMA} package provides classes and methods for dealing with maps in R.\\ \end{itemize} We consider in this document the use of the package {\tt adehabitatLT} to deal with the analysis of animal movements. All the methods available in \texttt{adehabitat} are also available in \texttt{adehabitatLT}. Contrary to the other brother packages, the classes of data returned by the functions of \texttt{adehabitatLT} are the same as those implemented in the original package \texttt{adehabitat}. Indeed, the structure of these classes were described in a paper (Calenge et al. 2009).\\ Package {\tt adehabitatLT} is loaded by <>= library(adehabitatLT) @ <>= suppressWarnings(RNGversion("3.5.0")) set.seed(13431) @ \section{What is a trajectory?} \subsection{Two types of trajectories} We designed the class \texttt{ltraj} to store the movements of animals monitored using radio-tracking, GPS, etc. The rationale underlying the structure of this class is described precisely in Calenge et al. (2009). We summarize this rationale in this vignette.\\ Basically, the trajectory of an animal is the curve described by the animal when it moves. The sampling of the trajectory implies a step of discretization, i.e., the division of this continuous curve into a number of discrete ``steps'' connecting successive relocations of the animal (Turchin, 1998). Two main classes of trajectories can be distinguished:\\ \begin{itemize} \item \textbf{Trajectories of type I} are characterized by the fact that the time is not precisely known or not taken into account for the relocations of the trajectory;\\ \item \textbf{the trajectories of type II} are characterized by the fact that the time is known for each relocation. This type of trajectory may in turn be divided into two subtypes:\\ \begin{itemize} \item \textbf{regular trajectories}: these trajectories are characterized by a constant time lag between successive relocations;\\ \item \textbf{irregular trajectories}: these trajectories are characterized by a variable time lag between successive relocations;\\ \end{itemize} \end{itemize} Note that the functions of \texttt{adehabitatLT} are mainly designed to deal with type I or type II regular trajectories. Irregular trajectories are harder to analyze, as the descriptive parameters of these trajectories (see below) may not be compared when computed on different time lags. \subsection{Descriptive parameters of the trajectory} \label{sec:param} Marsh and Jones (1988) noted that a good description of the trajectory is achieved when the following criteria are fullfilled:\\ \begin{itemize} \item the description is achieved with a minimum set of relatively easily measured parameters; \item the relationships between these parameters are defined precisely (e.g., with the help of a model); \item the parameters and the relationships between them are sufficient to reconstruct characteristic tracks without loosing any of their significant properties.\\ \end{itemize} Based on a literature review (see Calenge et al. 2009), we have chosen to characterize all the trajectories by the following parameters: <>= par(mar=c(0.1,0.1,0.1,0.1)) plot(c(0,1), c(0,1), ty="n", axes=FALSE) x <- c(0.2, 0.3, 0.5, 0.7, 0.4) y <- c(0.5, 0.3, 0.35, 0.6, 0.9) points(x,y, pch=16, cex=2) points(x[1],y[1], pch=16, col="green") lines(x,y, lty=1) arrows(0.4, 0.9, 0.1, 0.8, lty=2, length=0.1) lines(c(0.5, 0.7),c(0.35, 0.6), lwd=6) lines(c(0.5, 0.7),c(0.35, 0.6), lwd=2, col="red") ## dx arrows(0.5, 0.32, 0.7, 0.32, code=3, length=0.15) text(0.6, 0.3, "dx") ## dy arrows(0.75, 0.35, 0.75, 0.6, code=3, length=0.15) text(0.77, 0.5, "dy") ## abs.angle ang <- atan2(0.25, 0.2) ang <- seq(0,ang, length=10) xa <- 0.05*cos(ang) + 0.5 ya <- 0.05*sin(ang) + 0.35 lines(c(0.5, 0.7),c(0.35, 0.35), lty=3, col="red") lines(xa, ya, col="red", lwd=2) text(0.56, 0.38, expression(alpha), col="red") ## rel.angle lines(c(0.3, 0.5),c(0.3, 0.35), col="blue", lwd=2) lines(c(0.5, 0.7),c(0.35, 0.4), lty=3, col="blue") ang1 <- atan2(0.25, 0.2) ang0 <- atan2(0.05, 0.2) ang <- ang0+seq(0, ang1-ang0, length=10) xa <- 0.1*cos(ang) + 0.5 ya <- 0.1*sin(ang) + 0.35 lines(xa, ya, col="blue", lwd=2) xa <- 0.107*cos(ang) + 0.5 ya <- 0.107*sin(ang) + 0.35 lines(xa, ya, col="blue", lwd=2) text(0.61, 0.43,expression(beta), col="blue") ## dist arrows(0.48, 0.37, 0.68, 0.62, code=3, length=0.1) text(0.53, 0.5, "dist, dt") ## R2n arrows(0.21, 0.5, 0.49, 0.35, col="darkgreen", code=3, length=0.1) text(0.35, 0.49, expression(R[n]^2)) ## The legend: text(0.2, 0.2, expression(paste(alpha, ": abs.angle")), col="red") text(0.2, 0.17, expression(paste(beta, ": rel.angle")), col="blue") ## x0, y0, t0 text(0.14, 0.5, expression(paste(x[0], ", ", y[0], ", ", t[0])), col="darkgreen") box() @ \begin{center} <>= <> <> <> @ \end{center} \begin{itemize} \item \texttt{dx, dy, dt}: these parameters measured at relocation $i$ describe the increments of the x and y directions and time between the relocations $i$ and $i+1$. Such parameters are often used in the framework of stochastic differential equation modelling (e.g. Brillinger et al. 2004, Wiktorsson et al. 2004);\\ \item \texttt{dist}: the distance between successive relocations is often used in animal movement analysis (e.g. Root and Kareiva 1984, Marsh and Jones 1988);\\ \item \texttt{abs.angle}: the absolute angle $\alpha_i$ between the $x$ direction and the step built by relocations $i$ and $i+1$ is sometimes used together with the parameter \texttt{dist} to fit movement models (e.g. Marsh and Jones 1988);\\ \item \texttt{rel.angle}: the relative angle $\beta_i$ measures the change of direction between the step built by relocations $i-1$ and $i$ and the step built by relocations $i$ and $i+1$ (often called ``turning angle''). It is often used together with the parameter \texttt{dist} to fit movement models (e.g. Root and Kareiva 1984, Marsh and Jones 1988);\\ \item \texttt{R2n}: the squared distance between the first relocation of the trajectory and the current relocation is often used to test some movements models (e.g. the correlated random walk, see the seminal paper of Kareiva and Shigesada, 1983).\\ \end{itemize} \subsection{Several bursts of relocations} Very often, animal monitoring leads to several ``bursts'' of relocations for each monitored animal. For example, a GPS collar may be programmed to return one relocation every ten minutes during the night and no relocation during the day. Each night corresponds to a burst of relocations for each animal. We designed the class \texttt{ltraj} to take into account this burst structure. \subsection{Understanding the class \texttt{ltraj}} \label{sec:ltraj} An object of class \texttt{ltraj} is created with the function \texttt{as.ltraj} (see the help page of this function). We will take an example to illustrate the creation of an object of class \texttt{ltraj}. First load the dataset \texttt{puechabonsp} from the package \texttt{adehabitatMA}: <<>>= data(puechabonsp) locs <- puechabonsp$relocs locs <- as.data.frame(locs) head(locs) @ The data frame \texttt{locs} contains the relocations of 4 wild boar monitored using radio-tracking at Puechabon (Near Montpellier, South of France). First the date needs to be transformed into an object of the class \texttt{POSIXct}.\\ \noindent \textit{Remark:} The class \texttt{POSIXt} is designed to store time data in R (see the very clear help page of \texttt{POSIXt}). This class extends two sub-classes:\\ \begin{itemize} \item \textbf{the class \texttt{POSIXlt}}: This class stores a date in a list containing several elements related to this date (day of the month, day of the week, day of the year, month, year, time zone, hour, minute, second).\\ \item \textbf{the class \texttt{POSIXct}}: This class stores a date in a vector, as the number of seconds passed since January, 1st, 1970 at 1AM. This class is more convenient for storing dates into a data frame.\\ \end{itemize} We will use the function \texttt{strptime} (see the help page of this function) to convert the date in \texttt{locs} into a \texttt{POSIXlt} object, as then \texttt{as.POSIXct} to convert it into the class \texttt{POSIXct}: <<>>= da <- as.character(locs$Date) head(da) da <- as.POSIXct(strptime(as.character(locs$Date),"%y%m%d", tz="Europe/Paris")) @ We can then create an object of class \texttt{ltraj} to store the wild boar movements: <<>>= puech <- as.ltraj(xy = locs[,c("X","Y")], date = da, id = locs$Name) puech @ The result is a list of class ltraj containing four bursts of relocations corresponding to four animals. The trajectory is of type II and is irregular. There are no missing values. \\ This object is actually a list containing 4 elements (the four bursts). Each element is a data frame. Have a look, for example, at the first rows of the first data frame: <<>>= head(puech[[1]]) @ The function \texttt{as.ltraj} has automatically computed the descriptive parameters described in section \ref{sec:param} from the x and y coordinates, and from the date. Note that \texttt{dx,dy,dist} are expressed in the units of the coordinates \texttt{x,y} (here, metres) and \texttt{abs.angle,rel.angle} are expressed in radians.\\ \noindent A graphical display of the bursts can be obtained simply by: <>= plot(puech) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \subsection{Two points of views: steps (ltraj) or points (data.frame)?} We noted in the previous section that, in \texttt{adehabitatLT}, we consider the trajectory as a collection of successive ``steps'' ordered in time. We will see later in this vignette that most functions of \texttt{adehabitatLT} deal with trajectories considered from this point of view. However, several users (in particular, many thanks to Mathieu Basille and Bram van Moorter) noted that although this point of view may be useful to manage and analyse trajectories, it may be too restrictive to allow an easy management of such data.\\ Actually, the trajectory data may also be considered as a set of successive \textit{points} (the relocations) ordered in time. At first sight, the distinction between these two models may seem trivial, but it is important to consider it in several cases.\\ For example, any slight change in the coordinates/date of a relocation will change the value of all derived statistics (\texttt{dt}, \texttt{dist}, etc.). In the previous versions of adehabitat, it was possible to change directly the values of coordinates/dates in the object, and then to compute again the steps characteristics thanks to the function \texttt{rec} (it is still possible in the present version, but not recommended).\\ For example, consider the object \texttt{puech} created in the previous section. Have a look at the first relocations of the first burst: <<>>= head(puech[[1]]) @ Imagine that we realize that the X coordinate of the second relocation is actually equal to 700146 instead of 700046: <<>>= puech2 <- puech puech2[[1]][2,1] <- 700146 head(puech2[[1]]) @ The coordinate has been changed, but the step characteristics are now incorrect. The function \texttt{rec} recompute these statistics according to these changes: <<>>= head(rec(puech2)[[1]]) @ Although the function \texttt{rec} can be useful for sporadic use, it is limited when a larger number of modifications is required on the relocations (e.g. filtering incorrect relocations when ``cleaning'' GPS monitoring data). This is where the class \texttt{ltraj} does not fit. For such work, it is more convenient to see the trajectory as a set of points located in both space and time. And for such operations, it is sometimes more convenient to work with data frames. Two functions are provided to quickly convert a ltraj to and from data.frames: the functions \texttt{ld} and \texttt{dl}.\\ The function \texttt{ld} allows to quickly convert an object of class \texttt{ltraj} to the class \texttt{data.frame}. Consider for example the object \texttt{puech} created in the previous section. We can quickly convert this object towards the class \texttt{data.frame}: <<>>= puech2 <- ld(puech) head(puech2) @ Note that the data frame contains all the descriptors of the steps. In addition, two variables \texttt{burst} and \texttt{id} allow to quickly convert this object back towards the class \texttt{ltraj}, with the function \texttt{dl}: <<>>= dl(puech2) @ Using \texttt{dl} and \texttt{ld} can be extremey useful during the first steps of the analysis, especially during data ``cleaning''. \section{Managing objects of class \texttt{ltraj}} \subsection{Cutting a burst into several segments} Now, let us analyse the object \texttt{puech} created in the previous section. We noted that the object \texttt{puech} was not regular: <<>>= is.regular(puech) @ \noindent The function \texttt{is.regular} returns a Boolean... such a result can be obtained from a regular trajectory where just one relocation is missing, or from a completely irregular trajectory... we need more precision!\\ Have a look at the value of \texttt{dt} according to the date, using the function \texttt{plotltr}. Because \texttt{dt} is measured in seconds and that no more than one relocation is coellected every day, we convert this time lag into days by dividing it by 24 (hours/day) $\times$ 3600 (seconds / hour): <>= plotltr(puech, "dt/3600/24") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent The wild boar Chou was monitored during two successive summers (1992 and 1993). We need to ``cut'' this burst into two ``sub-burst''. We will use the function \texttt{cutltraj} to proceed. We first define a function \texttt{foo} that returns \texttt{TRUE} when the time lag between two successive relocations is greater than 100 days: <<>>= foo <- function(dt) { return(dt> (100*3600*24)) } @ \noindent Then, we use the function \texttt{cutltraj} to cut any burst relocations with a value of \texttt{dt} such that \texttt{foo(dt)} is true, into several bursts for which no value of \texttt{dt} fullfills this criterion: <<>>= puech2 <- cutltraj(puech, "foo(dt)", nextr = TRUE) puech2 @ \noindent Now, note that the burst of \texttt{Chou} has been splitted into two bursts: the first burst corresponds to the monitoring of Chou during 1992, and the second burst corresponds to the monitoring of Chou during 1993. We can give more explicit names to these bursts: <<>>= burst(puech2)[3:4] <- c("Chou.1992", "Chou.1993") puech2 @ \noindent Note that the function \texttt{id()} can be used similarly to replace the IDs of the animals.\\ \subsection{Playing with bursts} The bursts in an object \texttt{ltraj} can be easily managed. For example, consider te object \texttt{puech2} created previously: <<>>= puech2 @ \noindent Imagine that we want to work only on the males (Brock, Calou and Jean). We can subset this object using a classical extraction function: <<>>= puech2b <- puech2[c(1,2,5)] puech2b @ \noindent Or, if we want to study the animals monitored in 1993, we may combine this object with the monitoring of Chou in 1993: <<>>= puech2c <- c(puech2b, puech2[4]) puech2c @ \noindent It is also possible to select the bursts according to their id of their burst id (see the help page of \texttt{Extract.ltraj} for additional information, and in particular the example section).\\ The function \texttt{which.ltraj} can also be used to identify the bursts satisfying a condition. For example, imagine that we want to identify the bursts where the distance between successive relocations was greater than 2000 metres at least once: <<>>= bu <- which.ltraj(puech2, "dist>2000") bu @ \noindent This data frame contains the ID, burst ID and relocation numbers satisfying the specified criterion. We can then extract the bursts satisfying this criterion: <<>>= puech2[burst(puech2)%in%bu$burst] @ \subsection{Placing the missing values in the trajectory} \label{sec:NA} Now, look again at the time lag between successive relocations: <>= plotltr(puech2, "dt/3600/24") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent The relocations have been collected daily, but there are many days during which this relocation was not possible (storm, lack of field workers, etc.). We need to add missing values to define a regular trajectory. To proceed, we will use the function \texttt{setNA}. We have to define a reference date: <<>>= refda <- strptime("00:00", "%H:%M", tz="Europe/Paris") refda @ This reference date will be used to check that each date in the object of class \texttt{ltraj} is separated from this reference by an integer multiple of the theoretical \texttt{dt} (here, one day), and place the missing values at the times when relocations should theoretically have been collected. We use the function \texttt{setNA}: <<>>= puech3 <- setNA(puech2, refda, 1, units = "day") puech3 @ \noindent The trajectories are now regular, but there are now a lot of missing values! \subsection{Rounding the timing of the trajectories to define a regular trajectory} \label{sec:round} In some cases, despite the fact that the relocations were expected to be collected to return a regular trajectory, a minor delay is sometimes observed in this timing (e.g. the GPS collar needs some time to relocate). For example, consider the monitoring of four ibex in the Belledonne Mountains (French Alps): <<>>= data(ibexraw) ibexraw @ \noindent There is a variable time lag between successive relocations. Look at the time lag between successive relocations: <>= plotltr(ibexraw, "dt/3600") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent The relocations should have been collected every 4 hours, but there are some missing values. Use the function \texttt{setNA} to place the missing values, as in the section \ref{sec:NA}. We define a reference date and place the missing values: <<>>= refda <- strptime("2003-06-01 00:00", "%Y-%m-%d %H:%M", tz="Europe/Paris") ib2 <- setNA(ibexraw, refda, 4, units = "hour") ib2 @ \noindent Even when filling the gaps with NAs, the trajectory is still not regular. Now, look again at the time lag between successive relocations: <>= plotltr(ib2, "dt/3600") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent We can see that the time lag is only slightly different from 4 hour. The function \texttt{sett0} can be used to ``round'' the timing of the coordinates: <<>>= ib3 <- sett0(ib2, refda, 4, units = "hour") ib3 @ \noindent The trajectory is now regular.\\ \noindent \textbf{Important note}: The functions \texttt{setNA} and \texttt{sett0} are to be used to define a regular trajectory from a nearly regular trajectory. \textbf{It is NOT intended to transform an irregular trajectory into a regular one} (many users of \texttt{adehabitat} asked this question). \subsection{A special type of trajectories: same duration} \label{sec:sdtraj} In some cases, an object of class \texttt{ltraj} contains several regular bursts of the same duration characterized by relocations collected at the same time (same time lags between successive relocations, same number of relocations). We can check whether an object of class ``ltraj'' is of this type with the function \texttt{is.sd}. For example, consider again the movement of 4 ibexes monitored using GPS, stored in an object of class \texttt{ltraj} created in the previous section: <<>>= is.sd(ib3) @ \noindent This object is not of the type \texttt{sd} (same duration). However, theoretically, all the trajectories should have been sampled at the same time points. It is regular, but there are mismatches between the time of the relocations: <<>>= ib3 @ \noindent This is caused by the fact that there are missing relocations at the beginning and/or end of the monitoring for several animals (A160 and A286). We can use the function \texttt{set.limits} to define the time of beginning and ending of the trajectories. This function adds NAs to the beginning and ending of the monitoring when required: <<>>= ib4 <- set.limits(ib3, begin = "2003-06-01 00:00", dur = 14, units = "day", pattern = "%Y-%m-%d %H:%M", tz="Europe/Paris") ib4 @ \noindent All the trajectory are now covering the same time period: <<>>= is.sd(ib4) @ \noindent \textbf{Remark}: in our example, all the bursts are covering exactly the same time period (all begin at the same time and date and all stop at the same time and date). However, the function \texttt{set.limits} is much more flexible. Imagine for example that we are studying the movement of an animal during the night, from 00:00 to 06:00. If we have one burst per night, then it is possible to define an object of class \texttt{ltraj}, type \texttt{sd}, containing several nights of monitoring, even if the nights of monitoring do not correspond to the same date. If we consider that all the bursts cover the same period, then it is still possible to use the function \texttt{set.limits} to define an object of type \texttt{sd} (this is explained deeply on the help page of \texttt{set.limits}).\\ It is then possible to store some parameters of \texttt{sd} objects into a data frame (with one relocation per row and one burst per column), using the function \texttt{sd2df}. For example, considering the distance between successive relocations: <<>>= di <- sd2df(ib4, "dist") head(di) @ \noindent This data frame can then be used to study the interactions or similarities between the bursts. \subsection{Metadata on the trajectories (Precision of the relocations, etc.)} Sometimes, additional information is available for each relocation, and we may wish to store this information in the object of class \texttt{ltraj}, to allow the analysis of the relationships between these additional variables and the parameters of the trajectory.\\ This meta information can be stored in the attribute \texttt{infolocs} of each burst. This should be defined when creating the object \texttt{ltraj}, \textit{but can also be defined later} (see section \ref{sec:ras} for an example). For example, load the dataset \texttt{capreochiz}: <<>>= data(capreochiz) head(capreochiz) @ This dataset contains the relocations of one roe deer monitored using a GPS collar in the Chize forest (Deux-Sevres, France). This dataset contains the x and y coordinates (in kilometres), the date, and several variables characterizing the precision of the relocations. Note that the date is already of class \texttt{POSIXct}. We now define the object of class \texttt{ltraj}, storing the variables \texttt{Dop, Status, Temp, Act, Conv} in the attribute \texttt{infolocs} of the object: <<>>= capreo <- as.ltraj(xy = capreochiz[,c("x","y")], date = capreochiz$date, id = "Roe.Deer", infolocs = capreochiz[,4:8]) capreo @ The object \texttt{capreo} can be managed as usual. The function \texttt{infolocs()} can be used to retrieve the attributes \texttt{infolocs} of the bursts building up a trajectory: <<>>= inf <- infolocs(capreo) head(inf[[1]]) @ The function \texttt{removeinfo} can be used to set the attribute \texttt{infolocs} of all bursts to \texttt{NULL}.\\ Note that it is required that:\\ \begin{itemize} \item all the bursts are characterized by the same variables in the attribute \texttt{infolocs}. For example, it is not possible to store only the variable \texttt{Dop} for one burst and only the variable \texttt{Status} for another burst into the same object;\\ \item each row of the data frame stored as attributes \texttt{infolocs} correspond to one relocation (that is, the number of rows of the attribute should be the same as the number of relocations in the corresponding burst).\\ \end{itemize} Most functions of the package \texttt{adehabitatLT} do manage this attribute. For example, the functions \texttt{cutltraj} and \texttt{plotltr} can be used by calling variables stored in this attribute (as well as many other functions). For example: <>= plotltr(capreo, "log(Dop)") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Analyzing the trajectories} In this section, we will describe several tools available in \texttt{adehabitatLT} to analyse a trajectory. \subsection{Randomness of the missing values} A first important point is the examination of the distribution of the missing values in the trajectory. Missing values are frequent in the trajectories of animals collected using telemetry (e.g., GPS collar may not receive the signal of the satellite at the time of relocation, due for example to the habitat structure obscuring the signal, etc.). As noted by Graves and Waller (2006), the analysis of the patterns of missing values should be part of trajectory analysis.\\ The package \texttt{adehabitatLT} provides several tools for this analysis. For example, consider the object \texttt{ib4} created in section \ref{sec:sdtraj}, and containing 4 bursts describing the movements of 4 ibexes in the Belledonne moutain. We can first test whether the missing values occur at random in the monitoring using the function \texttt{runsNAltraj}: <>= runsNAltraj(ib4) @ <>= <> @ \begin{figure}[htbp] \begin{center} <>= <> <> <> @ \end{center} \end{figure} In this case, no difference appears between the number of runs actually observed in our trajectories and the distribution of the number of runs under the hypothesis of a random distribution of the NAs. The hypothesis of a random distribution of the NAs seems reasonable here. The statistics used here for the test is the number of runs in the sequence of relocations. For example, the sequence reloc-NA-NA-reloc-reloc-reloc-NA-NA-NA-reloc contains 5 runs, 3 runs of successful relocations and 2 runs of missing values. Under the hypothesis of random distribution of the missing values in the sequence, the theoretical expectation and standard deviation of the number of runs is known. The runs test is a randomization test that compares the standardized value of the number of runs (i.e. (value-expectation)/(standard deviation)) to the distribution of values obtained after randomizing the distribution of the NA in the sequence. Thus, a negative value of this standardized number of runs indicates that the missing values tend to be clustered together in the sequence.\\ But now, consider the distribution of the missing values in the case of the monitoring of one brown bear using a GPS collar: <<>>= data(bear) bear @ \noindent This trajectory is regular. The bear was monitored during one month, with one relocation every 30 minutes. We now test for a random distribution of the missing values for this trajectory: <>= runsNAltraj(bear) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent In this case, the missing values are not distributed at random. Have a look at the distribution of the missing values: <>= plotNAltraj(bear) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Because of the high number of relocations in this trajectory, this graph is not very clear. So a better way to study the distribution of the missing values is to work directly on the vector indicating whether the relocations are missing or not. That is: <<>>= missval <- as.numeric(is.na(bear[[1]]$x)) head(missval) @ \noindent This vector can then be analyzed using classical time series methods (e.g. Diggle 1990). We do not pursue on this aspect, as this is not the aim of this vignette to describe time series methods. \subsection{Should we consider the time?} \subsubsection{Type II or type I?} Until now, we have only considered trajectories of type II (time recorded). However, a common approach to the analysis of animal movements is to consider the movement as a discretized curve, and to study the geometrical properties of this curve (e.g., Turchin 1998; Benhamou 2004). That is, even if the data collection implied the recording of the time, it is often more convenient to consider the monitored movement as a trajectory of type I. There are two ways to define a type I trajectory with the functions of \texttt{adehabitatLT}. The first is to set the argument \texttt{typeII=FALSE} when calling the function \texttt{as.ltraj}. The second is to use the function \texttt{typeII2typeI}. For example, considering the trajectory of the bear loaded in the previous section, we can transform it into a type I object by: <<>>= bearI <- typeII2typeI(bear) bearI @ Nothing has changed, except that the time is replaced by an integer vector ordering the relocations in the trajectory. \subsubsection{Rediscretizing the trajectory in space} Several authors have advised to rediscretize type I trajectories so that they are built by steps with a constant length (e.g. Turchin 1998). This is a convenient approach to the analysis, as all the geometrical properties of the trajectory can be summarized by studying the variation of the relative angles.\\ The function \texttt{redisltraj} can be used for this rediscretization. For example, look at the trajectory of the brown bear stored in \texttt{bearI} (created in the previous section): <>= plot(bearI) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent Now, rediscretize this trajectory with constant step length of 500 metres: <<>>= bearIr <- redisltraj(bearI, 500) bearIr @ \noindent The number of relocations has increased. Have a look at the rediscretized trajectory: <>= plot(bearIr) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent Then, the geometrical properties can be studied by examining the distribution of the relative angles. For example, the function \texttt{sliwinltr} can be used to smooth the cosine of relative angle using a sliding window method: <>= sliwinltr(bearIr, function(x) mean(cos(x$rel.angle)), type="locs", step=30) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The beginning of the trajectory is characterized by mean cosine close to 0.5 (tortuous trajectory). Then the movements of the animal is more linear (i.e., less tortuous). A finer analysis should now be done on these data. So that we need to get the relative angles from this rediscretized trajectory: <<>>= cosrelangle <- cos(bearIr[[1]]$rel.angle) head(cosrelangle) @ This vector can now be analyzed using classical time series analysis methods. We do not pursue this analysis further, as this is beyond the scope of this vignette. \subsubsection{Rediscretizing the trajectory in time} \label{sec:redist} A common way to deal with missing values consist in linear interpolation. That is given relocations $r_1 = (x_1, y_1, t_1)$ and $r_3 = (x_3, y_3, t_2)$, it is possible to estimate the relocation $r_2$ collected at $t_2$ with: \begin{eqnarray} x_2 & = & \frac{(t_2 - t_1)}{(t_3 - t_1)} (x_3 - x_1) \nonumber \\ y_2 & = & \frac{(t_2 - t_1)}{(t_3 - t_1)} (y_3 - y_1) \nonumber \end{eqnarray} \textbf{Such an interpolation may be of limited value when many relocations are missing} (because it supposes implicitly that the animal is moving along a straight line). However, it can be useful to ``fill in'' a small number of relocations. Such an interpolation can be carried out with the function \texttt{redisltraj}, setting the argument \texttt{type = "time"}. In this case, this function rediscretizes the trajectory so that the time lag between successive relocations is constant. For example, consider the data set \texttt{porpoise}: <<>>= data(porpoise) porpoise @ We focus on the first three animals. David and Mitchell contain a small number of missing relocations, so that it can be a good idea to interpolate these relocations, to facilitate later calculations. We first remove the missing values with \texttt{na.omit} and then interpolate the trajectory so that the time lag between successive relocations is constant: <<>>= (porpoise2 <- redisltraj(na.omit(porpoise[1:3]), 86400, type="time")) @ The object \texttt{porpoise2} now does not contain any missing values. Note that the rediscretization step should be expressed in seconds. \subsection{Dynamic exploration of a trajectory} The package \texttt{adehabitatLT} provides a function very useful for the dynamic exploration of animal movement: the function \texttt{trajdyn}. This function allows to dynamically zoom/unzoom, measure distance between several bursts or several relocations, explore the trajectory in space and time, etc. For example, the ibex data set is explored by typing: <>= trajdyn(ib4) @ \begin{center} \includegraphics[height=10cm,keepaspectratio]{trajdyn} \end{center} \noindent Note that this function can draw a background defined by an object of class \texttt{SpatialPixelsDataFrame} or \texttt{SpatialPolygonsDataFrame}. \subsection{Analyzing autocorrelation} Dray et al. (2010) noted that the analysis of the sequential autocorrelation of the descriptive parameters presented in section \ref{sec:param} is essential to the understanding of the mechanisms driving these movements. The approach proposed by these authors is implemented in \texttt{adehabitatLT}. In this section, we describe the functions that can be used to carry out this kind of analysis.\\ A positive autocorrelation of a parameter means that values taken near to each other tend to be either more similar (positive autocorrelation) or less similar (negative autocorrelation) than would be expected from a random arrangement.\\ \subsubsection{Testing for autocorrelation of the linear parameters} The independence test of Wald and Wolfowitz (1944) is implemented in the generic function \texttt{wawotest}. Basically, this function can be used to test the sequential autocorrelation in a vector. However, the method \texttt{wawotest.ltraj} allows to test autocorrelation for the three \textit{linear} parameters \texttt{dx}, \texttt{dy} and \texttt{dist} for each burst in an object of class \texttt{ltraj}. For example, consider again the monitoring of movements of the bear: <<>>= wawotest(bear) @ Note that this function removes the missing values before the test. The row \texttt{p} indicates the P-value of the test. We can see that the three linear parameters are strongly positively autocorrelated. There are periods during which the animal is traveling at large speed and periods when the animals are walking at lower speed. Note that this was already apparent on the graph showing the changes with time of the distance between successive relocations: <>= plotltr(bear, "dist") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \noindent This was even clearer on the graph showing the moving average of the distance (with a window of 5 days): <>= sliwinltr(bear, function(x) mean(na.omit(x$dist)), 5*48, type="locs") @ \noindent (not executed in this report). \subsubsection{Analyzing the autocorrelation of the parameters} The autocorrelation function (ACF) $\rho(a)$ measures the correlation between a parameter measured at time $t$ and the same parameter measured at time $t-a$ in the same time series (Diggle, 1990). This allows to analyze the autocorrelation, and to identify the scales at which this autocorrelation occurs. Dray et al. (2010) noted that the autocorrelation function measured at lag 1 ($\rho(1)$) is mathematically equivalent to the independence test of Wald and Wolfowitz (1944).\\ Dray et al. (2010) extended the mathematical bases underlying the ACF to handle the missing data occurring frequently in the trajectories. Their approach is implemented in the function \texttt{acfdist.ltraj} (the management of NAs is described in detail on the help page of this function). For example, consider again the monitoring of a brown bear: <>= acfdist.ltraj(bear, lag=5, which="dist") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We have calculated here the ACF for the distance for a time lag up to 5 relocations. The interested reader can try to calculate the ACF for trajectories up to 100 relocations to see the cyclic patterns occurring in this trajectory. \subsubsection{Testing autocorrelation of the angles} The test of the autocorrelation of the angular parameters (relative or absolute angles, see section \ref{sec:param}) is based on the chord distance between successive angles (see Dray et al. 2010 for additional details): <>= opar <- par(mar = c(0,0,4,0)) plot(0,0, asp=1, xlim=c(-1, 1), ylim=c(-1, 1), ty="n", axes=FALSE, main="Criteria f for the measure of independence between successive angles at time i-1 and i") box() symbols(0,0,circle=1, inches=FALSE, lwd=2, add=TRUE) abline(h=0, v=0) x <- c( cos(pi/3), cos(pi/2 + pi/4)) y <- c( sin(pi/3), sin(pi/2 + pi/4)) arrows(c(0,0), c(0,0), x, y) lines(x,y, lwd=2, col="red") text(0, 0.9, expression(f^2 == 2*sum((1 - cos(alpha[i]-alpha[i-1])), i==1, n-1)), col="red") foo <- function(t, alpha) { xa <- sapply(seq(0, alpha, length=20), function(x) t*cos(x)) ya <- sapply(seq(0, alpha, length=20), function(x) t*sin(x)) lines(xa, ya) } foo(0.3, pi/3) foo(0.1, pi/2 + pi/4) foo(0.11, pi/2 + pi/4) text(0.34,0.18,expression(alpha[i]), cex=1.5) text(0.15,0.11,expression(alpha[i-1]), cex=1.5) par(opar) @ \begin{center} <>= <> <> <> @ \end{center} For example, the function \texttt{testang.ltraj} is a randomization test using the mean squared chord distance as a criteria. For example, considering again the trajectory of the bear, we can test the autocorrelation of the relative angles between successive moves: <<>>= testang.ltraj(bear, "relative") @ \subsubsection{Analyzing the autocorrelation of angular parameters} Dray et al. (2010) extended the ACF to angular parameters by considering the chord distance as a criteria to build the ACF. This approach is implemented in the function \texttt{acfang.ltraj}. Considering again the \texttt{bear} dataset: <>= acfang.ltraj(bear, lag=5) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We can see that the relative angle observed at time $i$ is significantly correlated with the angle observed at time $i-1$. \subsection{Segmenting a trajectory into segments characterized by a homogenous behaviour} \subsubsection{The method of Gueguen (2001)} We implemented a new approach to the segmentation of movement data, relying on a Bayesian partitioning of a sequence. This approach was originally developed in molecular biology, to partition DNA sequences (Gueguen 2001). We describe this approach in this section.\\ Biologically, a positive autocorrelation in any of the descriptive parameters may mean that the animal behaviour is changing with time (there are periods during which the animal is feeding, other during which the animal is resting, etc.). The idea is then to segment the trajectory of the animal into homogenous segments.\\ We will use the movements of a porpoise monitored using an Argos collar to illustrate this method. First we load the data: <<>>= data(porpoise) gus <- porpoise[1] gus @ \noindent The trajectory is regular and is built by relocations collected every 24 hours during two months. Plot the data: <>= plot(gus) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Visually, \textbf{the trajectory seems to be built by three segments}. At the very beginning of the trajectory, the animal is performing very short moves. Then, the animal is travelling faster toward the southwest, and finally, the animal is again performing very small moves.\\ We can draw the ACF for the distance between successive relocations to illustrate the autocorrelation pattern from another point of view: <>= acfdist.ltraj(gus, "dist", lag=20) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} There is a strong autocorrelation pattern present in the data, up to lag 8. We can plot the distances between successive relocations according to the date <>= plotltr(gus, "dist") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Now, let us suppose that the distances between successive relocations have been generated by a normal distribution, with different means corresponding to different behaviours. Let us built 10 models corresponding to 10 values of the mean distance ranging from 0 to 130 km/day: <<>>= (tested.means <- round(seq(0, 130000, length = 10), 0)) @ Based on the visual exploration of the distribution of distance, we set the standard deviation of the distribution to 5 km. We can now define 10 models characterized by 10 different values of means and with a standard deviation of 5 km: <<>>= (limod <- as.list(paste("dnorm(dist, mean =", tested.means, ", sd = 5000)"))) @ The approach of Gueguen (2001) allows, based on these a priori models, to find both the number and the limits of the segments building up the trajectory. Any model can be supposed for any parameter of the steps (the distance, relative angles, etc.), provided that the model is Markovian.\\ Given the set of \textit{a priori} models, for a given step of the trajectory, it is possible to compute the probability density that the step has been generated by each model of the set. The function \texttt{modpartltraj} computes the matrix containing the probability densities associated to each step (rows), under each model of the set (columns): <<>>= mod <- modpartltraj(gus, limod) mod @ Then, we can estimate the optimal number of segments in the trajectory, given the set of \textit{a priori} models, using the function \texttt{bestpartmod}, taking as argument the matrix \texttt{mod}: <>= bestpartmod(mod) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} This graph presents the value of the log-likelihood (y) that the trajectory is actually made of $K$ segments (x). Note that this log-likelihood is actually corrected using the method of Gueguen (2001) (which implies the Monte Carlo simulation of the independence of the steps in the trajectory -- explaining the boxplots --, see the help page of \texttt{bestpartmod} for further details on this procedure). In this case, the method indicates that 4 segments are a reasonable choice for the segmentation. \textbf{This is a surprise for us, as we rather expected 3 segments} (actually, the number of segments returned by the function depend on the models supposed \textit{a priori}).\\ Finally, the function \texttt{partmod.ltraj} can be used to compute the segmentation. The mathematical rationale underlying these two functions is the following: given an optimal $k$-partition of the trajectory, if the $i^{th}$ step of the trajectory belongs to the segment $k$ predicted by the model $d$, then either the relocation $i-1$ belongs to the same segment, in which case the segment containing $i-1$ is predicted by $d$, or the relocation $i-1$ belongs to another segment, and the other $(k-1)$ segments together constitute an optimal $(k-1)$ partition of the trajectory [1--$(i-1)$]. These two probabilities are computed recursively by the functions from the matrix \texttt{mod}, observing that the probability of a 1-partition (partition built by one segment) of the trajectory from 1 to $i$ described by the model $m$ is simply the product of the probability densities of the steps from 1 to i under the model m.\\ \noindent \textbf{Remark:} this approach relies on the hypothesis of the independence of the steps within each segment.\\ Now, use the function \texttt{partmod.ltraj} to partition the trajectory of the porpoise into 4 segments: <<>>= (pm <- partmod.ltraj(gus, 4, mod)) @ We can see that the models at the beginning of the trajectory and at the end of the trajectory are the same. Have a look at this segmentation: <>= plot(pm) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} This is very interesting: we already noted that the movements at the very beginning and the end of the trajectory were much slower that the rest of the trajectory, and this is confirmed by this partition. However, this segmentation illustrates a change of speed at the middle of the ``migration''. The end of the migration is much faster than the beginning. This is clearer on the graph showing the changes in distance between successive relocations with the date. Let us plot this graph together with the segmentation:\\ <>= ## Shows the segmentation on the distances: plotltr(gus, "dist") tmp <- lapply(1:length(pm$ltraj), function(i) { coul <- c("red","green","blue")[as.numeric(factor(pm$stats$mod))[i]] lines(pm$ltraj[[i]]$date, rep(tested.means[pm$stats$mod[i]], nrow(pm$ltraj[[i]])), col=coul, lwd=2) }) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The end of the migration is nearly two times faster than the beginning of the migration.\\ To conclude, have a look at the residuals of this segmentation: <>= ## Computes the residuals res <- unlist(lapply(1:length(pm$ltraj), function(i) { pm$ltraj[[i]]$dist - rep(tested.means[pm$stats$mod[i]], nrow(pm$ltraj[[i]])) })) plot(res, ty = "l") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} And a Wald and Wolfowitz test suggests that the residuals of this segmentation are independent, confirming the validity of the approach: <<>>= wawotest(res) @ \subsubsection{The method of Lavielle (1999, 2005)} Barraquand and Benhamou (2008) proposed an approach to partition the trajectory, to identify the places where the animal stays a long time (based on the calculation of the residence time: see the help page of the function \texttt{residenceTime}). Once the residence time has been calculated for each relocation, they propose to use the method of Lavielle (1999, 2005) to partition the trajectory. We describe this method in this section.\\ The method of Lavielle \textit{per se} finds the best segmentation of a time series, given that it is built by $K$ segments. It searches the segmentation for which a contrast function (measuring the contrast between the actual series and the model underlying the segmented series) is minimized. Let $Y_i$ be the value of the focus variable (e.g. the speed, the residence time, etc.) at time $i$, and $n$ the number of observations. We suppose that the data have been generated by the following model: $$ Y_i = \mu_i + \sigma_i \epsilon_i $$ \noindent where $\mu_i$ and $\sigma_i$ are respectively the mean and the standard deviation of $Y_i$, and $\epsilon_i$ is a sequence of zero mean random variables with unit variance. \textit{Note that it is not required that the $\epsilon_i$ are independent}. The model underlying the segmentation relies on the hypothesis that the parameters $\mu_i$ and $\sigma_i$ are constant within a segment, but also that they vary between successive segments. Several contrast functions can be used to measure the discrepancy between the observed series and a given segmentation model, depending on assumptions on this variation:\\ \begin{itemize} \item we can suppose \textit{a priori} that only the mean $\mu_i$ changes between segments, and that $\sigma_i = \sigma$ is the same for all segments. Then, for a given partition of the series built by $K$ segments with known limits, the following function can be used to measure the discrepancy between the observed series and the model supposing a mean constant within segments and changing from one segment to the next: $$ J_K(Y) = \sum_k G^1_k(Y_{i, i \in k}) $$ where $Y_{i, i \in k}$ is the set of observations in the segment $k$, and $$ G^1_k(Y_{i, i \in k}) = \sum_{i=t^k_1}^{t^k_{n(k)}} (Y_i - \bar Y_k)^2 $$ with $t^k_1$ and $t^k_{n(k)}$ the indices of the first and last observations of segment $k$ respectively, and $\bar Y_k$ the sample mean of the observations in the segment $k$;\\ \item we can suppose \textit{a priori} that only the standard deviation $\sigma_i$ changes between segments, and that $\mu_i = \mu$ is the same for all segments. Then, for a given partition of the series built by $K$ segments with known limits, the following function can be used to measure the discrepancy between the observed series and the model supposing a variance constant within segments and changing from one segment to the next: $$ J_K(Y) = \sum_k G^2_k(Y_{i, i \in k}) $$ where $$ G^2_k(Y_{i, i \in k}) = \frac{1}{n(k)} \log \left ( \frac{1}{n(k)} \sum_{i=t^k_1}^{t^k_{n(k)}} (Y_i - \bar Y)^2 \right ) $$ with $n(k)$ the number of observations in the segment $k$ and $\bar Y$ the mean of the whole series;\\ \item finally, we can suppose \textit{a priori} that both the standard deviation $\sigma_i$ and the mean $\mu_i$ change between segments. Then, for a given partition of the series built by $K$ segments with known limits, the following function can be used to measure the discrepancy between the observed series and the model supposing a mean and variance constant within segments but changing from one segment to the next: $$ J_k(Y) = \sum_k G^3_k(Y_{i, i \in k}) $$ where $$ G^3_k(Y_{i, i \in k}) = \frac{1}{n(k)} \log \left ( \frac{1}{n(k)} \sum_{i=t^k_1}^{t^k_{n(k)}} (Y_i - \bar Y_k)^2 \right ) $$ \end{itemize} The three contrast functions can be used to measure the contrast between a given series and the model underlying a given segmentation of this series. \textit{For a given number of segments} $K$, the method of Lavielle consists in the search of the best segmentation of the series, i.e. the segmentation for which the chosen contrast function is minimized.\\ For example, consider a series of 10 observations. First, we calculate an upper triangular $10\times 10$ matrix containing the value of the contrast function for all possible segments built by the successive observations of the series. This matrix is called the contrast matrix: <>= par(mar=c(0,0,0,0)) plot(c(0,1), xlim=c(-0.3,1.3), ylim=c(-0.3,1.3), axes=FALSE) ii <- seq(0.05, 0.95, by=0.1) tmp <- lapply(ii, function(i) { lapply(ii, function(j) { if (j-(1-i)> -0.01) { polygon(c(i-0.05, i+0.05, i+0.05, i-0.05), c(j-0.05, j-0.05, j+0.05, j+0.05), col="grey") } else { polygon(c(i-0.05, i+0.05, i+0.05, i-0.05), c(j-0.05, j-0.05, j+0.05, j+0.05)) } })}) polygon(c(0,1,1,0), c(0,0,1,1), lwd=2) text(rep(-0.05,10), ii, 10:1) text(ii, rep(1.05,10), 1:10) i <- 0.65 j <- 0.85 polygon(c(i-0.05, i+0.05, i+0.05, i-0.05), c(j-0.05, j-0.05, j+0.05, j+0.05), col="red") segments(0.4, 1.18, i, j) points(i,j, pch=16) text(0.4, 1.25, "Contrast function measured on the segment\n beginning at observation 2 and ending at observation 7") @ Given a number $K$ of segments, the Lavielle method consists in the search of the best path from the first to the last observation in this matrix. For example, consider the following path: <>= par(mar=c(0,0,0,0)) plot(c(0,1), xlim=c(-0.3,1.3), ylim=c(-0.3,1.3), axes=FALSE) ii <- seq(0.05, 0.95, by=0.1) tmp <- lapply(ii, function(i) { lapply(ii, function(j) { if (j-(1-i)> -0.01) { polygon(c(i-0.05, i+0.05, i+0.05, i-0.05), c(j-0.05, j-0.05, j+0.05, j+0.05), col="grey") } else { polygon(c(i-0.05, i+0.05, i+0.05, i-0.05), c(j-0.05, j-0.05, j+0.05, j+0.05)) } })}) polygon(c(0,1,1,0), c(0,0,1,1), lwd=2) text(rep(-0.05,10), ii, 10:1) text(ii, rep(1.05,10), 1:10) lines(c(0.00, 0.45, 0.50, 0.75, 0.80, 0.95), c(0.95, 0.95, 0.45, 0.45, 0.15, 0.15), col="red", lwd=3) points(c(0.45, 0.75, 0.95), c(0.95, 0.45, 0.15), pch=16, cex=2) @ This segmentation is built by 3 segments:\\ \begin{itemize} \item the first segment begins at the first observation and ends at the fifth observation. The value of the contrast function for this segment is at the intersection between the first row and the fifth column. \item the second segment begins at the sixth observation and ends at the eighth observation. The value of the contrast function for this segment is at the intersection between the sixth row and the eighth column. \item the last segment begins at the ninth observation and ends at the tenth observation. The value of the contrast function for this segment is at the intersection between the ninth row and the tenth column.\\ \end{itemize} The value of the contrast function for the whole segmentation is the sum of the values located at the dots. For a given value of the number of segments (in this example $K=3$), a dynamic programming algorithm allows to identify the best segmentation, i.e. the path through the matrix for which the contrast function is minimized.\\ \noindent \textit{Remark}: in practice, a minimum number of observations $L_{min}$ in the segments should be specified. For example, if $L_{min} = 3$, no segment will contain less than 3 observations.\\ Now, remains the question of the choice of the optimal number of segments $K_{opt}$. Several approaches have been proposed to guide the choice of $K_{opt}$. The approaches proposed by Lavielle (2005) rely on the examination of the decrease of the contrast function $J(K)$ with the number of segments $K$. Indeed, the more there are segments, and the smaller the contrast function will be (because a segmentation built by more segments will fit more closely to the actual series). However, there should be a clear "break" in the decrease of this function after the optimal value $K_{opt}$. Therefore, the number of segments can be chosen graphically, based on the examination of the decrease of this contrast function J(K) with $K$. Lavielle (2005) also suggested an alternative way to estimate automatically the optimal number of segments, also relying on the presence of a "break" in the decrease of the contrast function. He proposed to choose the last value of $K$ for which the second derivative $D(K)$ of a \textit{standardized} constrast function is greater than a threshold $S$ (see Lavielle, 2005 for details). Based on numerical experiments, he proposed to choose the value $S = 0.75$.\\ Note that the standardization of the contrast function is based on the specification of a value of $K_{max}$, the maximum number of segments expected by the user (see Lavielle 2005 for details). Note that the value of $K_{max}$ chosen by the user may have a strong effect on the resulting value of $K_{opt}$, in particular for small time series (i.e. less than 500 observations). Barraquand (com. pers) noted from simulations that values of $K_{max}$ set to about 3 or 4 times the expected value of $K_{opt}$ seem to give the best results. \\ The Lavielle method is implemented in the function \texttt{lavielle}. For example, consider the following time series, built by 6 segments of 100 observations each, drawn from a normal distribution with a mean varying from one segment to the next (0, 2, 0, -3, 0, 2): <>= suppressWarnings(RNGversion("3.5.0")) set.seed(129) seri <- c(rnorm(100), rnorm(100, mean=2), rnorm(100), rnorm(100, mean=-3), rnorm(100), rnorm(100, mean=2)) plot(seri, ty="l", xlab="time", ylab="Series") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We will use the Lavielle method to see whether we can identify the limits used to generate this series. We can use the function \texttt{lavielle} to calculate the contrast matrix: <<>>= (l <- lavielle(seri, Lmin=10, Kmax=20, type="mean")) @ Note that we indicate that each segment should contain at least 10 observations (\texttt{Lmin=10}) and that the whole series is built by at most 20 segments (\texttt{Kmax = 20}). The result of the function is a list of class \texttt{"lavielle"} containing various results, but the most interesting one here is the contrast matrix. We can have a look at the decrease of the contrast function when the number $K$ of segments increase, with the function \texttt{chooseseg}: <>= chooseseg(l) @ This function returns: (i) the value of the contrast function \texttt{Jk} for various values of $K$, (ii) the value of the second derivative of the contrast function \texttt{D} for various values of $K$, and (iii) a graph showing the decrease of \texttt{Jk} with $K$. The slope of the contrast function is strongly negative when $K<6$, but there is a sharp break at $K = 6$. Visually, the value of $K_{opt} = 6$ seems therefore reasonable. Moreover, the last value of $K$ for which $D>0.75$ is $K=6$, which confirms our first visual choice (the second derivative is actually very strong for this value of $K$). We can have a look at the best segmentation with 6 segments with the function \texttt{findpath}: <>= fp <- findpath(l, 6) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The Lavielle method does a pretty good job in this ideal situation. Note that \texttt{fp} is a list containing the indices of the first and last observations of the series building the segment: <<>>= fp @ \noindent The limits found by the method are very close to the actual limits here.\\ The function \texttt{lavielle} is a generic function that has a method for objects of class \texttt{"ltraj"}. For example, consider again the porpoise described in the previous section. We will perform a segmentation of this trajectory based on the distances travelled by the animal from one relocation to the next. <>= plotltr(gus, "dist") @ We can use the function \texttt{lavielle} to partition this trajectory into segments with homogeneous speed. A priori, based on a visual examination, there is no reason to expect that this series is built by more than 20 segments. We will use the Lavielle method, with $K_{max} = 20$ and $L_{min} = 2$ (we suppose that there are at least 2 observations in a segment): <>= lav <- lavielle(gus, Lmin=2, Kmax=20) chooseseg(lav) @ There is a clear break in the decrease of the contrast function after $K=4$. In addition, this is also the last value of $K$ for which $D>0.75$. Look at the segmentation of the series of distances: <>= kk <- findpath(lav, 4) kk @ Note that the red lines indicate the beginning of new segments. The function \texttt{findpath} here returns an object of class \texttt{"ltraj"} with one burst per segment: <>= plot(kk) @ Note that this segmentation is nearly identical to the one found with the method of Gueguen (2001).\\ \noindent \textit{Important Remark}: The Lavielle method implies the calculation of a contrast matrix with dimensions equal to $n \times n$, where $n$ is the number of observations in the series. However, as $n$ increases, the size of the contrast matrix may exceed the memory limit of the computer. In this case, a possibility consists in the calculation of the contrast function along a grid of size $l_d$. Therefore the contrast function will not be calculated for all possible segments in the series, but for all possible segments containing multiple of $l_d$ observations. For example, if $l_d = 3$, possible segments beginning by observation 1 are (1,2,3), (1, 2, ..., 6), (1, 2, ..., 9), (1, 2, ..., 12), etc. Therefore, $l_d$ defines the resolution of the segmentation: the segment limits can only occur for observations located at indices multiple of $l_d$. This sub-sampling approach is possible with the function \texttt{lavielle}, by setting the parameter \texttt{ld} to a value greater than 1. \textit{Note that in this case, \texttt{Lmin} should necessarily be a multiple of \texttt{ld}} (the function fails otherwise).\\ \noindent \textit{Remark}: The Lavielle method was originally propose to partition a trajectory based on the residence time of the animal. The residence time method is implemented in the function \texttt{residenceTime}. The use of this function will add a new variable in the component \texttt{"infolocs"} of the object of class \texttt{"ltraj"} that can be used with the Lavielle method (see the help page of this function). \subsection{Rasterizing a trajectory} \label{sec:ras} In some cases, it may be useful to rasterize a trajectory. In particular, when the aim of the study is to examine the habitat traversed by the animal, this approach may be useful. For example, consider the dataset \texttt{puechcirc}, containing 3 trajectories of 2 wild boars. It may be useful to identify the habitat traversed by the animal during each step. A habitat map is available in the dataset \texttt{puechabonsp}: <>= data(puechcirc) data(puechabonsp) mimage(puechabonsp$map) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We can rasterize the trajectories of the wild boars: <<>>= ii <- rasterize.ltraj(puechcirc, puechabonsp$map) @ The result is a list containing 3 objects of class "SpatialPointsDataFrame" (one per animal). Let us examine the first one: <<>>= tr1 <- ii[[1]] head(tr1) @ This data frame contains the coordinates of the pixels traversed by each step. For example, the rasterized trajectory for the first animal is: <>= plot(tr1) points(tr1[tr1[[1]]==3,], col="red") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The red points indicate the pixels traversed by the third step. These results can be used to identify the habitat characteristics of each step. For example, we may calculate the mean elevation for each step. To proceed, we use the function \texttt{over} of the package sp. <>= ov <- over(tr1, puechabonsp$map) mo <- tapply(ov[[1]], tr1[[1]], mean) plot(mo, ty="l") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Here, we can see that the first animal stays on the plateau at the beginning at the monitoring, then goes down to the crops, and goes back to the plateau. It is easy to repeat the operation for all the animals. We will make use of the \texttt{infolocs} attribute. We first build a list containing data frames, each data frame containing on variable describing the mean elevation traversed by the animal between relocation $i-1$ and relocation $i$: <<>>= val <- lapply(1:length(ii), function(i) { ## get the rasterized trajectory tr <- ii[[i]] ## over with the map ov <- over(tr, puechabonsp$map) ## calculate the mean elevation mo <- tapply(ov[[1]], tr[[1]], mean) ## prepare the output elev <- rep(NA, nrow(puechcirc[[i]])) ## place the average values at the right place ## names(mo) contains the step number (i.e. relocation ## number +1) elev[as.numeric(names(mo))+1] <- mo df <- data.frame(elevation = elev) ## same row.names as in the ltraj row.names(df) <- row.names(puechcirc[[i]]) return(df) }) @ \noindent Then, we define the \texttt{infolocs} attribute: <<>>= infolocs(puechcirc) <- val @ \noindent and finally, we can plot the mean elevation as a function of date: <>= plotltr(puechcirc, "elevation") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \subsection{Models of animal movements} Several movement models have been proposed in the litteratured to describe animal movements. The package \texttt{adehabitatLT} contains several functions allowing to simulate these models. Such simulations can be very useful to test hypotheses concerning a trajectory, because all the descriptive parameters of the steps are also generated by the functions. Actually, the package proposes 6 functions to simulate such models:\\ \begin{itemize} \item \texttt{simm.brown} can be used to simulate a Brownian motion; \item \texttt{simm.crw} can be used to simulate a correlated random walk. This model has been often used to describe animal movements (Kareiva and Shigesada 1983); \item \texttt{simm.mba} can be used to simulate an arithmetic Brownian motion (with a drift parameter and a covariance between the coordinates, see Brillinger et al. 2002); \item \texttt{simm.bb} can be used to simulate a Brownian bridge motion (i.e. a Brownian motion constrained by a fixed start and end point); \item \texttt{simm.mou} can be used to simulate a bivariate Ornstein-Uhlenbeck motion (often used to describe the sedentarity of an animal, e.g. Dunn and Gipson 1977); \item \texttt{simm.levy} can be used to simulate a Levy walk, as described (Bartumeus et al. 2005).\\ \end{itemize} All these functions return an object of class \texttt{ltraj}. For example, simulate a correlated random walk built by 1000 steps characterized by a mean cosine of the relative angles equal to 0.95 and a scale parameter for the step length equal to 1 (see the help page of \texttt{simm.crw} for additional details on the meaning of these parameters): <<>>= sim <- simm.crw(1:1000, r=0.95) sim @ Note that the vector \texttt{1:1000} passes as an argument is considered here as a vector of dates (it is converted to the class \texttt{POSIXct} by the function, see section \ref{sec:ltraj} for more details on this class). Other dates can be passed to the functions. Have a look at the simulated trajectory: <>= plot(sim, addp=FALSE) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Null models of animal movements} \subsection{What is a null model?} The package \texttt{adehabitatLT} provides several functions to perform a null model analysis of trajectories. Null models are frequently used in community ecology to test hypotheses about the processes that generated the data. Gotelli and Graves (1996) defined the null model as ``\textit{a pattern generating model that is based on randomization of ecological data or random sampling from a known or imagined distribution. The null model is designed with respect to some ecological or evolutionary process of interest. Certain elements of the data are held constant, and others are allowed to vary stochastically to create new assemblage patterns. The randomization is designed to produce a pattern that would be expected in the absence of a particular ecological mechanism}''.\\ Gotelli (2001) gave a more detailed description of the null model approach: ``\textit{to build a null model, an index of community structure, such as the amount of niche overlap (...), is first measured for the real data. Next, a null community is generated according to an algorithm or set of rules for randomization, and this same index is measured for the null community. A large number of null communities are used to generate a frequency histogram of index values expected if the null hypothesis is true. The position of the observed index in the tails of this null distribution is then used to assign a probability value to the pattern, just as in a conventional statistical analysis}.\\ Although mostly used in community ecology, this approach was also advocated for trajectory data, e.g. in the analysis of habitat selection (Martin et al. 2008) and the study of the interactions between animals (Richard et al. 2012). This approach can indeed be very useful to test hypotheses related to animal movements. The package adehabitatLT propose several general null models that can be used to test biological hypotheses. \\ \subsection{The problem} For example, consider the first animal in the object \texttt{puechcirc} loaded previously. We plot this trajectory on an elevation map: <>= data(puechcirc) data(puechabonsp) boar1 <- puechcirc[1] xo <- coordinates(puechabonsp$map) ## Note that xo is a matrix containing the coordinates of the ## pixels of the elevation map (we will use it later to define ## the limits of the study area). plot(boar1, spixdf=puechabonsp$map, xlim=range(xo[,1]), ylim=range(xo[,2])) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} At first sight, it seems that the animal occupies a large range of values of elevation. We may want to test whether the variance of elevation values at animal relocations could have been observed under the hypothesis of random habitat use. The question is therefore: what do we mean by ``random habitat use''? We can precise our aim here. We may consider the animal as a random sample of all animals living on the study area. If we focus on the second scale of habitat selection as defined by Johnson (1980), i.e. on the selection of the home range in the study area, under the hypothesis of random habitat use, the observed trajectory could have been located anywhere on the study area.\\ Therefore, if we want to know if the actual elevation variance could have been observed under the hypothesis of random habitat use, one possibility is to simulate this kind of random habitat use a large number of times, and to calculate the elevation variance for each simulated data set. If the observed value is far from the distribution of simulated values, this would allow to rule out the null model.\\ One possibility to simulate this kind of ``random habitat use'' for this animal could be to randomly rotate and shift the observed trajectory over the study area. Rotating and shifting the trajectory as a whole allows to keep the trajectory structure unchanged (therefore taking into account the internal constraints of the animal, see Martin et al. 2008). The function \texttt{NMs.randomShiftRotation} from the package \texttt{adehabitatLT} allows to define this type of null model (but other null models are also available in \texttt{adehabitatLT}, see section \ref{sec:imprem}), and the function \texttt{testNM} allows to simulate it.\\ \subsection{The principle of the approach} The basic principle of the null model approach implemented in the package is the following:\\ \begin{itemize} \item first write a treatment function that will be used to calculate a criterion characterizing the process under study from the simulated data (see below for how to write such a function);\\ \item then choose a null model and define the constraints to be satisfied by simulated datasets. Define this null model with one function \texttt{NMs.*} (see section \ref{sec:imprem} or the help page of \texttt{testNM} for a list of the available null models), and simulate $N$ random data sets generated by this model using the function \texttt{testNM};\\ \item finally compare the observed criterion with the distribution of simulated values to determine whether the observed data set could have been generated by the null model.\\ \end{itemize} The function \texttt{testNM}, used to simulate the model, generates at each step of the simulation process a data frame \texttt{x} with three columns: the coordinates x, y, and the date (whatever the chosen type of null model). The treatment function defines how to calculate the criterion used in the test from this data frame. The treatment function can be any function defined by the user, but \textbf{must} take arguments \texttt{x} and \texttt{par}, where:\\ \begin{itemize} \item \texttt{x} is the data frame generated by the null model;\\ \item \texttt{par} can be any R object (e.g. a list) containing the parameters required by the treatment function. Note that this argument is needed even if no parameters are required by the treatment function. In this case, the argument \texttt{par} will not be used by the function, but must be present in the function definition).\\ \end{itemize} We will define a treatment function that calculate the elevation variance from a simulated trajectory later. For the moment, we will define a treatment function that just plot the simulated trajectories on an elevation map. Consider the following function: <<>>= plotfun <- function(x, par) { image(par) points(x[,1:2], pch=16) lines(x[,1:2]) return(x) } @ This function will be used to plot the simulations of the null model. In this case, the argument \texttt{par} correspond to a map of the study area. Note that this function also returns the data frame \texttt{x}. Now, define the following null model: <>= nmo <- NMs.randomShiftRotation(na.omit(boar1), rshift = TRUE, rrot = TRUE, rx = range(xo[,1]), ry = range(xo[,2]), treatment.func = plotfun, treatment.par = puechabonsp$map[,1], nrep=9) nmo @ We have removed the missing values from the object of class \texttt{ltraj} (this function does not accept missing values in the trajectory). The arguments \texttt{rshift} and \texttt{rrot} indicate that we want to randomly rotate the trajectory and shift it over the study area. The study area is necessarily a rectangle defined by its x and y limits \texttt{rx} and \texttt{ry}. We indicate that the treatment function is the function \texttt{plotfun} that we just wrote, and that the argument \texttt{par} that will be passed to the treatment function (i.e. \texttt{treatment.par}) is the elevation map. We only define 9 repetitions of the null model. Now, we simulate the model using the function \texttt{testNM}: <>= suppressWarnings(RNGversion("3.5.0")) set.seed(90909) par(mfrow=c(3,3), mar=c(0,0,0,0)) resu <- testNM(nmo, count = FALSE) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} This figure illustrates the datasets generated by the null model. We can see that each data set corresponds to the original trajectory after a random rotation and shift over the study area. First, we can see a problem: some of the generated trajectories are located outside the study area (because the study area is necessarily defined as a rectangle in this type of null model). It here necessary to define a \textit{constraint function} when defining the null model. Like the treatment function, the constraint function takes two arguments \texttt{x} and \texttt{par}, and \textit{must return a logical value} indicating whether the constraint(s) are satisfied or not. In our example, we would like that all the relocations building the trajectory fall within the study area. In other words, if we overlap spatially the relocations and the elevation map, there should be no missing value. Define the following constraint function: <>= confun <- function(x, par) { ## Define a SpatialPointsDataFrame from the trajectory coordinates(x) <- x[,1:2] ## overlap the relocations x to the elevation map par jo <- join(x, par) ## checks that there are no missing value res <- all(!is.na(jo)) ## return this check return(res) } @ Now, define again the null model, but also define the constraint function: <<>>= nmo2 <- NMs.randomShiftRotation(na.omit(boar1), rshift = TRUE, rrot = TRUE, rx = range(xo[,1]), ry = range(xo[,2]), treatment.func = plotfun, treatment.par = puechabonsp$map[,1], constraint.func = confun, constraint.par = puechabonsp$map[,1], nrep=9) @ Now, if we simulate the null model, only those data sets satisfying the constraint will be returned by the function: <>= suppressWarnings(RNGversion("3.5.0")) set.seed(90909) par(mfrow=c(3,3), mar=c(0,0,0,0)) resu <- testNM(nmo2, count = FALSE) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \subsection{Back to the problem} Now consider again our problem: is the variance of the elevation values at the relocations greater than expected under the hypothesis of random use of space? We can write a treatment function that calculates the variance of elevation values from a data frame \texttt{x} containing the relocations and a SpatialPixelsDataFrame \texttt{par} containing the elevation map: <<>>= varel <- function(x, par) { coordinates(x) <- x[,1:2] jo <- join(x, par) return(var(jo)) } @ We can define again a null model to calculate a distribution of values of variance expected under the null model. We use the function \texttt{varel} as treatment function in the null model. To save some time, we calculate only 99 values under this null model, but the user is encouraged to try the function with a larger value: <<>>= nmo3 <- NMs.randomShiftRotation(na.omit(boar1), rshift = TRUE, rrot = TRUE, rx = range(xo[,1]), ry = range(xo[,2]), treatment.func = varel, treatment.par = puechabonsp$map[,1], constraint.func = confun, constraint.par = puechabonsp$map[,1], nrep=99) @ Note that we define the same constraint function as before (all relocations should be located inside the study area. We now simulate the null model: <>= sim <- testNM(nmo3, count=FALSE) @ <>= load("sim.Rdata") @ Now, calculate the observed value of variance: <<>>= (obs <- varel(na.omit(boar1)[[1]], puechabonsp$map[,1])) @ And compare this observation with the distribution obtained under the null model, using the function \texttt{as.randtest} from the package ade4: <>= (ran <- as.randtest(unlist(sim[[1]]), obs)) plot(ran) @ The P-value is rather low, which seems to indicate that the range of elevations used by the boar is important in comparison to what would be expected under the hypothesis of random habitat use. \subsection{Null model with several animals} Now, consider again the data set \texttt{puechcirc}. This data set contains three trajectories of two wild boars: <>= puechcirc plot(puechcirc) @ Note that the two trajectories of CH93 are roughly located at the same place on the study area. We decide to bind these two trajectories into one: <<>>= (boar <- bindltraj(puechcirc)) @ Now, we can reproduce the null model analysis separately for each animal. When the object of class \texttt{ltraj} passed as argument contains several trajectories, the simulations are performed separately for each one. Therefore, to define the null model for all animals in \texttt{boar}, we can use the same command line as before, just replacing \texttt{boar1} by \texttt{boar}: <<>>= nmo4 <- NMs.randomShiftRotation(na.omit(boar), rshift = TRUE, rrot = TRUE, rx = range(xo[,1]), ry = range(xo[,2]), treatment.func = varel, treatment.par = puechabonsp$map[,1], constraint.func = confun, constraint.par = puechabonsp$map[,1], nrep=99) @ We now simulate this null model with the function \texttt{testNM}: <>= sim2 <- testNM(nmo4, count=FALSE) @ <>= load("sim2.Rdata") @ \texttt{sim2} is a list with two components (one per trajectory), each component being itself a list with \texttt{nrep=99} elements (the elevation variance calculated for each of the 99 datasets generated by the null model). We can calculate the observed elevation variance for each observed trajectory: <<>>= (obs <- lapply(na.omit(boar), function(x) { varel(x, puechabonsp$map[,1]) })) @ And calculate a P-value for each animal separately: <<>>= lapply(1:2, function(i) { as.randtest(unlist(sim2[[i]]), obs[[i]]) }) @ In this case, none of the two tests are significant at the conventional $\alpha = $ 5\% level. \subsection{Single null models and multiple null models} In the previous section, we showed how to build a null model and simulate this null model for several trajectories separately. Now, we may find more convenient to design a global criterion (measured on all animals) to test whether the elevation variance is greater than expected under the null model. This is the principle of ``multiple null models''. This approach is the following:\\ \begin{itemize} \item First define a global criterion that will be used to test the hypothesis under study; \item Define ``single null models'' using any of the functions \texttt{NMs.*} (see below), i.e. the randomization approach and the constraints that will be used to simulate a trajectory for each animal; \item Define a ``multiple null model'' from the ``single null models'', by defining the constraints that should be satisfied by each set of simulated trajectories and the treatment function that will be applied to each set; \item Simulate the null model \texttt{nrep} times to obtain \texttt{nrep} values of the criterion under the defined null model, using the function \texttt{testNM} \item Calculate the observed value of the criterion and calculate a P-value by comparing the observed value to the distribution of values expected under the null model. \end{itemize} We illustrate this approach on our example. We first define, as a global criterion, the mean elevation variance, i.e. the elevation variance averaged over all animals. We need to define a treatment function allowing to calculate this global criterion. The treatment function should take two arguments named \texttt{x} and \texttt{par}. The argument \texttt{x} should be an object of class \texttt{"ltraj"} (i.e. the set of trajectories simulated by the null model at each step of the randomization process) and the argument \texttt{par} can be any R object (e.g. a list) containing the parameters needed by the treatment function. In our example, the treatment function will be the following: <<>>= meanvarel <- function(x, par) { livar <- lapply(x, function(y) { coordinates(y) <- y[,1:2] jo <- join(y, par) return(var(jo)) }) mean(unlist(livar)) } @ We have defined the single null models in the previous section. We now define a multiple null model from this object using the function \texttt{NMs2NMm}: <<>>= nmo5 <- NMs2NMm(nmo4, treatment.func = meanvarel, treatment.par = puechabonsp$map, nrep = 99) @ Note that both the treatment function and the number of repetitions that we have defined in the function \texttt{NMs.randomShiftRotation} will be ignored when the multiple null model will be simulated. We can now simulate the model: <>= sim3 <- testNM(nmo5, count=FALSE) @ <>= load("sim3.Rdata") @ At each step of the randomization process, two trajectories are simulated (under the null model and satisfying the constraints) and the treatment function is applied to the simulated trajectories. The result is therefore a list with \texttt{nrep} component, each component containing the result of the treatment function. We now calculate the observed criterion: <<>>= (obs <- meanvarel(na.omit(boar), puechabonsp$map)) @ And we define a randomization test from these results, using the function \texttt{as.randtest} from the package ade4: <>= (ran <- as.randtest(unlist(sim3), obs)) plot(ran) @ As a conclusion, we are unable to highlight any difference between the observed mean elevation variance and the distribution of values expected under the null model. \subsection{Several important remarks} \label{sec:imprem} We give in this section several important remarks:\\ \begin{itemize} \item We have illustrated null models in this vignette with the help of the function \texttt{NMs.randomShiftRotation}. However, several other null models are available in the package. The function \texttt{NMs.randomCRW} can be used to simulate a correlated random walk, with the distribution of turning angles and distances between successive relocations calculated from the data. The function \texttt{NMs.CRW} implements a purely parametric correlated random walk (see \texttt{?simm.crw}), where parameters are specified by the used. Finally, the function \texttt{NMs.randomCs} is similar to the function \texttt{NMs.randomShiftRotation}, but give a finer control on the distances between trajectory barycenter and a set of capture sites. The help page of \texttt{testNM} given further details and examples of use.\\ \item We have illustrated the use of null models to test for habitat selection, but it can be used for a wide variety of hypotheses (interactions between animals, etc.).\\ \item The main argument of the functions \texttt{NMs.*} is an object of class \texttt{ltraj}. However, in some cases, the null models can be useful to test hypotheses on other types of data. For example, imagine that we want to measure the overlap between the home-range of several animals monitored using radio-tracking, and that, for each animal, only the X and Y coordinates are available. It is therefore impossible, in theory, to define an object of class \texttt{ltraj} (since a date is needed for such objects), though it might be biologically sensible to compare the observed overlap with the overlap expected under a random distribution of the animals. It is still possible to use the functions \texttt{NMs.randomShiftRotation} and \texttt{NMs.randomCs} for this kind of analysis. In this case, the user can create a ``fake'' date variable (e.g. if there are \texttt{nr} rows in the data frame containing the X and Y coordinates, a ``fake'' date variable can be created by \texttt{fa <- 1:nr}, and can be used to create an object of class \texttt{ltraj} that can be used in null model analysis.\\ \item Note that for multiple null models, two kinds of constraints are possible: it is possible to define constraints for the individual trajectories in the function \texttt{NMs.*} (e.g. the simulated trajectory should fall within the study area), and for the whole set of trajectories in the function \texttt{NMs2NMm} (e.g. at least 80\% of the trajectories should be located in a given habitat type, see Richard et al., 2012 for an example). The two types of constraint function are taken into account when simulating the null model.\\ \item Note that missing relocations are not allowed in null model analysis. Therefore, if there are missing relocations in the trajectory, the use of \texttt{na.omit.ltraj} is required prior to the analysis. \end{itemize} \section{Conclusion and perspectives} Several other methods can be used to analyze a trajectory. Thus, the first passage time method, developed by Fauchald and Tveraa (2003) to identify the areas where \textit{area restricted search} occur is implemented in the function \texttt{fpt}. Several methods are available in the package \texttt{adehabitatHR} to estimate a home range based on objects of class \texttt{ltraj}. Thus, the Brownian bridge kernel method (Bullard 1999, Horne et al. 2007), the biased random bridge kernel method (Benhamou and Cornelis 2010, Benhamou 2011), and the product kernel algorithm (Keating and Cherry 2009) are implemented in the functions \texttt{kernelbb} and \texttt{kernelkc} respectively.\\ But one thing is important: at many places in this vignette, we have noted that the descriptive parameters of the steps can be analysed as a (possibly multiple) time series. The R environment provides many functions to perform such analyses, and \textbf{we stress that the package \texttt{adehabitatLT} should be considered as a springboard toward such functions}.\\ We included in the package \texttt{adehabitatLT} several functions allowing the analysis of animal movements. All the brother packages \texttt{adehabitat*} contain a vignette similar to this one, which explains not only the functions, but also in some cases the philosophy underlying the analysis of animal space use. \section*{References} \begin{description} \item Bartumeus, F., da Luz, M.G.E., Viswanathan, G.M. and Catalan, J. 2005. Animal search strategies: a quantitative random-walk analysis. Ecology, 86: 3078--3087. \item Barraquand, F. and Benhamou, S. 2008. Animal movements in heterogeneous landscapes: identifying profitable places and homogeneous movement bouts. Ecology, 89, 3336--3348. \item Benhamou, S. and Cornelis, D. 2010. Incorporating movement behavior and barriers to improve kernel home range space use estimates. Journal of Wildlife Management, 74, 1353--1360. \item Benhamou, S. 2011. Dynamic approach to space and habitat use based on biased random bridges. PLOS ONE, 6, 1--8. \item Benhamou, S. 2004. How to reliably estimate the tortuosity of an animal's path: straightness, sinuosity, or fractal dimension? Journal of Theoretical Biology, 229, 209--220. \item Brillinger, D., Preisler, H., Ager, A., Kie, J. and Stewart, B. 2002. Employing stochastic differential equations to model wildlife motion. Bulletin of the Brazilian Mathematical Society, 2002, 33, 385-408. \item Brillinger, D., Preisler, H., Ager, A. and Kie, J. 2004. An exploratory data analysis (EdA) of the paths of moving animals. Journal of Statistical Planning and Inference, 122, 43-63. \item Bullard, F. 1999. Estimating the home range of an animal: a Brownian bridge approach Johns Hopkins University. \item Calenge, C. 2005. Des outils statistiques pour l'analyse des semis de points dans l'espace ecologique. Universite Claude Bernard Lyon 1. \item Calenge, C. 2006. The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals. Ecological modelling, 197, 516--519. \item Calenge, C., Dray, S. and Royer-Carenzi, M. 2009. The concept of animals' trajectories from a data analysis perspective. Ecological Informatics, 4, 34--41. \item Diggle, P. 1990. Time series. A biostatistical introduction Oxford University Press. \item Dray, S., Royer-Carenzi, M. and Calenge, C. 2010. The exploratory analysis of autocorrelation in animal-movement studies. Ecological Research, 4, 34-41. \item Dunn, J. and Gipson, P. 1977. Analysis of radio telemetry data in studies of home range. Biometrics, 33, 85-101 \item Fauchald, P. and Tveraa, T. 2003. Using first-passage time in the analysis of area-restricted search and habitat selection/ Ecology, 84, 282-288. \item Gotelli, N. 2001. Research frontiers in null model analysis. Global ecology and biogeography, 10, 337--343. \item Gotelli, N. and Graves, G. 1996. Null models in Ecology. Smithsonian Institution Press. \item Graves, T. and Waller, J. 2006. Understanding the causes of missed global positioning system telemetry fixes. Journal of Wildlife Management, 70, 844--851. \item Gueguen, L. 2001. Segmentation by maximal predictive partitioning according to composition biases. Pp 32-44 in: Gascuel, O. and Sagot, M.F. (Eds.), Computational Biology, LNCS, 2066. \item Horne, J., Garton, E., Krone, S. and Lewis, J. 2007. Analyzing animal movements using Brownian bridges. Ecology, 88, 2354--2363. \item Kareiva, P. and Shigesada, N. 1983. Analysing insect movement as a correlated random walk. Oecologia, 56, 234--238. \item Keating, K. and Cherry, S. 2009. Modeling utilization distributions in space and time. Ecology, 90, 1971--1980. \item Lavielle, M. 2005. Using penalized contrasts for the change-point problem. Signal Processing, 85, 1501--1510. \item Lavielle, M. 1999. Detection of multiple changes in a sequence of dependent variables. Stochastic Processes and their Applications. 83, 79--102. \item Marsh, L. and Jones, R. 1988. The form and consequences of random walk movement models. Journal of Theoretical Biology, 133, 113--131. \item Martin, J., Calenge, C., Quenette, P.Y. and Allain{\'e}, D. 2008. Importance of movement constraints in habitat selection studies. Ecological Modelling, 213, 257--262. \item Pebesma, E. and Bivand, R.S. 2005. Classes and Methods for Spatial data in R. R News, 5, 9--13. \item Richard, E., Calenge, C., Said, S., Hamann, J.L. and Gaillard, J.M. (2012) A null model approach to study the spatial interactions between roe deer (\textit{Capreolus capreolus}) and red deer (\textit{Cervus elaphus}). Ecography, in press. \item Root, R. and Kareiva, P. 1984. The search for resources by cabbage butterflies (Pieris Rapae): Ecological consequences and adaptive significance of markovian movements in a patchy environment. Ecology, 65, 147--165. \item Turchin, P. 1998. Quantitative Analysis of Movement: measuring and modeling population redistribution in plants and animals. Sinauer Associates, Sunderland, MA. \item Wald, A. and Wolfowitz, J. 1944. Statistical Tests Based on Permutations of the Observations The Annals of Mathematical Statistics, 15, 358--372. \item Wiktorsson, M., Ryden, T., Nilsson, E. and Bengtsson, G. 2004. Modeling the movement of a soil insect. Journal of Theoretical Biology, 231, 497--513. \end{description} \end{document} % vim:syntax=tex