Behavioral Change Point Analysis

Description

A collection of functions that allows one to perform the behavioral change point analysis (BCPA) as described by Gurarie et al. (2009, Ecology Letters, 12: 395-408). The key features are estimation of discrete changes in time-series data, notable linear and turning components of gappy velocity times series extracted from movement data.

Details

For more movement-appropriate change point analysis, users are encouraged to apply correlated velocity change point analysis as implemented in the smoove package (as of this writing on GitHub at https://github.com/EliGurarie/smoove) which implements methods described in Gurarie et al. 2017.

There is a fairly detailed vignette - type in vignette("bcpa"), and an updated vignette vignette("UnivariateBCPA") with a univariate example. The key analysis function is WindowSweep, and reading its documentation is a good way to start using this package.

WindowSweep uses a suite of functions that might be useful for more narrow analysis. These functions are all available and documented, and are listed here hierarchically, from the most fundamental to the most complex:

• GetRho maximizes the likelihood to estimate autocorrelation rho or characteristic time-scale tau.

• GetDoubleL estimates the parameters and returns the log-likelihood at either side of a given break

• GetBestBreak finds the single best change point according to the likelihood returned by GetDoubleL'](GetDoubleL)

• GetModels uses a (modified) BIC model selection for all combinations from M0 (\mu_1 = \mu_2, \sigma_1 = \sigma_2, \rho_1 = \rho_2) to M7 (\mu_1 \neq \mu_2, \sigma_1 \neq \sigma_2, \rho_1 \neq \rho_2) to characterize the "Best Break"

• Finally, WindowSweep sweeps a longer time series with the Best Break / Model Selection analysis, identifying most likely break points and BIC selected models across the time series.

Summary, diagnostic, and plotting functions are:

• PartitionParameters - estimated parameters of a BCPA.

• ChangePointSummary - summary table of the change points.

• plot.bcpa - a plotting method for visualizing the time series with vertical lines as change points.

• PathPlot - a method for drawing a color-coded track of the analysis.

• DiagPlot - diagnostic plots for BCPA.

A few preprocessing functions are also available:

• plot.track - method for plotting a generic "track" object.

• GetVT - returns step-lengths, absolute and turning angles from track data.

Author(s)

Eliezer Gurarie

References

• Gurarie, E., R. Andrews and K. Laidre. 2009. A novel method for identifying behavioural changes in animal movement data. Ecology Letters. 12: 395-408.

• Gurarie, E., C. Fleming, W.F. Fagan, K. Laidre, J. Hernandez-Pliego, O. Ovaskainen. 2017. Correlated velocity models as a fundamental unit of animal movement: synthesis and applications. Movement Ecology, 5:13.

Examples

# Running through a complete analysis:
## loading data
  data(Simp)
## plotting the track (using the plot.track method)
  plot(Simp)
## Obtaining the movement summary table (with turning angles and step lengths)
  Simp.VT <- GetVT(Simp)
## Applying the analysis
 if(interactive()){
  Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, 
                         windowstep = 1, progress=TRUE)
                         } else
 Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, 
                         windowstep = 1, progress=FALSE)
                                            
## plotting outpots
  plot(Simp.ws, threshold=7)
  plot(Simp.ws, type="flat", clusterwidth=3)
  PathPlot(Simp, Simp.ws)
  PathPlot(Simp, Simp.ws, type="flat")
## Diagnostic of assumptions
  DiagPlot(Simp.ws)

Obtain summary of BCPA analysis

Description

Produces a summary of change points for a "flat" analysis, identifying phases (periods between change points) with estimated parameters, clustering neighboring ones according to a kernel density of the windowsweep breaks.

Usage

ChangePointSummary(windowsweep, clusterwidth = 1, tau = TRUE)

Arguments

Value

a list containing two elements:

Author(s)

Eliezer Gurarie

See Also

WindowSweep

Examples

if(!exists("Simp.VT")){
 data(Simp)
 Simp.VT <- GetVT(Simp)}
if(!exists("Simp.ws"))
 Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, windowstep = 1, progress=TRUE)
# too many change points:
ChangePointSummary(Simp.ws)
# about the right number of change points:
ChangePointSummary(Simp.ws, clusterwidth=3)

Diagnostic plot for BCPA

Description

Draws diagnostic plots for BCPA analysis. Specifically: a qqplot, a histogram (with a N(0,1) density curve), and an acf of the standardized residuals of an analysis.

Usage

DiagPlot(
  windowsweep,
  type = c("smooth", "flat")[1],
  plotme = TRUE,
  values = FALSE,
  ...
)

Arguments

Value

If values is TRUE, returns the values of the standardized residuals.

Author(s)

Eliezer Gurarie

See Also

PartitionParameters

Examples

data(Simp)
if(!exists("Simp.VT"))
 Simp.VT <- GetVT(Simp)
if(!exists("Simp.ws"))
 Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, windowstep = 1, progress=TRUE)
DiagPlot(Simp.ws)
DiagPlot(Simp.ws, type="flat")
# The Simp's diagnostic plots are excellent.

Find most likely change point in irregular time series

Description

Finds the single best change point according to the likelihood function. Used internally within WindowSweep.

Usage

GetBestBreak(x, t, range = 0.6, ...)

Arguments

Value

returns a single row (vector) with elements: breaks,tbreaks,mu1,sigma1,rho1,LL1,mu2,sigma2,rho2,LL2,LL. The breakpoint is calculated for a range of possible values of width range*l (where l is the length of the time series). The output of this function feeds WindowSweep.

Author(s)

Eliezer Gurarie

See Also

WindowSweep which uses it, and GetDoubleL for the likelihood estimation.

Examples

# An example with a single break:
x <- c(arima.sim(list(ar = 0.9), 20) + 10, arima.sim(list(ar = 0.1), 20)) 
t <- 1:length(x)
plot(t,x, type="l")
(bb <- GetBestBreak(x,t, tau=FALSE))
abline(v = bb[2], col=2)

Obtain log-likelihood and parameter estimates for a given break point.

Description

Takes a time series with values x obtained at time t and a time break tbreak, and returns the estimates of \mu, \sigma and \tau (or \rho) as well as the negative log-likelihood of those estimates before and after the break. Mostly for use internally within GetBestBreak.

Usage

GetDoubleL(x, t, tbreak, ...)

Arguments

Value

a vector containing the parameters and the negative log-likelihoods in order: mu1, sigma1, tau1, LL1, mu2, sigma2, tau2, LL2

Author(s)

Eliezer Gurarie

See Also

GetBestBreak uses this function, while this function uses GetRho

Obtain likelihood estimates of gappy Gaussian time series

Description

Obtain likelihood of gappy standardized Gaussian time series "x" sampled at times "t" given parameter "rho" (autocorrelation). Alternatively computes the characteristic time scale "tau".

Arguments

Value

Returns the log-likelihood of the data.

Author(s)

Eliezer Gurarie

See Also

Core function of BCPA, used directly in GetRho

Examples

# simulate autocorrelated time series
  rho.true <- 0.8
  x.full <- arima.sim(1000, model=list(ar = rho.true))
  t.full <- 1:1000
  
# subsample time series
  keep <- sort(sample(1:1000, 200))
  x <- x.full[keep]
  t <- t.full[keep]
  plot(t,x, type="l")
  
# Obtain MLE of rho
  rhos <- seq(0,.99,.01)
  L <- sapply(rhos, function(r) GetL(x, t, r))
  rho.hat <- rhos[which.max(L)]
  plot(rhos, L, type = "l")
  abline(v = c(rho.true, rho.hat), lty=3:2, lwd=2)
  legend("bottomleft", legend=c("true value","MLE"), lty=3:2, lwd=2, 
         title = expression(rho))
         
# Why tau is better
  tau.true <- -1/log(rho.true)
  taus <- seq(1,10,.1)
  L <- sapply(taus, function(r) GetL(x, t, r, tau = TRUE))
  tau.hat <- taus[which.max(L)]

  plot(taus, L, type = "l")
  abline(v = c(tau.true, tau.hat), lty=3:2, lwd=2)
  legend("bottomleft", legend=c("true value","MLE"), lty=3:2, lwd=2, 
         title = expression(tau))

Model selection at a known breakpoint

Description

Returns all parameter estimates and log-likelihoods for all possible models at a selected breakpoint. These are:

• M0 - all parameters equal

• M1 - \mu_1 != \mu_2

• M2 - \sigma_1 != \sigma_2

• M3 - \tau_1 != \tau_2

• M4 - \mu_1 != \mu_2 AND \sigma_1 != \sigma_2

• M5 - \mu_1 != \mu_2 AND \tau_1 != \tau_2

• M6 - \sigma_1 != \sigma_2 AND \tau_1 != \tau_2

• M7 - all parameters unequal

Usage

GetModels(x, t, breakpoint, K = 2, tau = TRUE)

Arguments

Value

Returns a names matrix with 8 rows (one for each model) and columns: Model, LL, bic, mu1, s1, rho1, mu2, s2, rho2. Fairly self-explanatory. Note that the rho columns include the tau values, if tau is TRUE (as it usually should be).

Author(s)

Eliezer Gurarie

See Also

Used directly within WindowSweep. Relies heavily on GetRho.

Characteristic time / auto-correlation for irregular time series

Description

Estimates characteristic time \tau or auto-correlation \rho from a gappy time series dataset. It works first by estimating the mean and standard deviation directly from the time series X, using these to standardize the time series, and then optimizes for the likelihood of the value of \tau or \rho.

Usage

GetRho(x, t, tau = TRUE)

Arguments

Value

Returns a vector of length two: the estimate and the negative log-likelihood

Author(s)

Eliezer Gurarie

Obtain VT table from Track

Description

The VT table computes speeds, step lengths, orientations, turning angles and ancillary variables from a track. The output of this function is (typically) meant to feed the WindowSweep function.

Usage

GetVT(Data, units = "hour", skiplast = TRUE)

Arguments

Value

a data frame with the following columns:

Author(s)

Eliezer Gurarie

See Also

WindowSweep

Examples

data(Simp)
plot(Simp)
Simp.VT <- GetVT(Simp)
head(Simp.VT)
# Distribution of estimated speeds
hist(Simp.VT$V, col="grey", breaks=20)
# Distribution of turning angles

if (requireNamespace("circular", quietly = TRUE)) 
   circular::rose.diag(Simp.VT$Theta, bins=24) else
     hist(Simp.VT$Theta)

Make Track

Description

Simple convenience function for creating a track class object from X-Y-Time movement data. A track class object can be conveniently plotted and analyzed within bcpa.

Usage

MakeTrack(X, Y, Time)

Arguments

Value

a track class data frame, with three columns: X, Y and Time.

See Also

plot.track

Examples

X <- cumsum(arima.sim(n=100, model=list(ar=0.8)))
Y <- cumsum(arima.sim(n=100, model=list(ar=0.8)))
Time <- 1:100
mytrack <- MakeTrack(X,Y,Time)
plot(mytrack)

Partition parameters

Description

Partitions - and, ultimately, estimates - all parameters of a BCPA, either as a rolling average (smooth BCPA), or as constant values within fixed change points (flat BCPA).

Usage

PartitionParameters(
  windowsweep,
  type = c("smooth", "flat")[1],
  clusterwidth = 1
)

Arguments

Value

a data frame containing the three estimates: mu.hat, s.hat, rho.hat.

Author(s)

Eliezer Gurarie

See Also

used in ChangePointSummary and PhasePlot

Path plot of BCPA output

Description

Plots the animal's trajectory, with segments color-coded according to the time scale / auto-correlation of the BCPA output, and width of segments proportional to the estimated mean of the BCPA.

Usage

PathPlot(
  Data,
  windowsweep,
  type = c("smooth", "flat")[1],
  clusterwidth = 1,
  plotlegend = FALSE,
  tauwhere = "topleft",
  n.legend = 5,
  ncol.legend = 1,
  bty.legend = "n",
  ...
)

Arguments

Author(s)

Eliezer Gurarie

Examples

if(!exists("Simp.ws"))
{
 data(Simp)
 Simp.ws <- WindowSweep(GetVT(Simp), "V*cos(Theta)", windowsize = 50, 
 windowstep = 1, progress=TRUE)
}

PathPlot(Simp, Simp.ws, plotlegend=TRUE, n.legend=3)
PathPlot(Simp, Simp.ws, type="flat", clusterwidth=3, plotlegend=TRUE)

Phase plot of BCPA output

Description

Behavioral phase plot of BCPA output. Mean and standard deviation are on the x and y axis. Size and color of points represent the time scale of autocorrelation

Usage

PhasePlot(
  windowsweep,
  type = c("smooth", "flat")[1],
  clusterwidth = 1,
  ...,
  legend.where = "bottomright"
)

Arguments

Author(s)

Eliezer Gurarie

See Also

WindowSweep, PartitionParameters

Examples

if(!exists("Simp.ws"))
{
 data(Simp)
 Simp.ws <- WindowSweep(GetVT(Simp), "V*cos(Theta)", windowsize = 50, windowstep = 1, progress=TRUE)
}

PhasePlot(Simp.ws)

Simulated Chim (Simp) Data

Description

This is a simulated movement track from a continuous auto-correlated process of total duration 60 time units with four behavioral phases that switch at times T=10, 20 and 50 from, variously, higher or lower velocity and longer or shorter characteristic time scale of autocorrelation. The original data was simulated with time intervals of 0.01 (leading to a complete track with 6000 data points). 200 points were randomly sampled from the "true" movement track, such that the intervals between locations are random with mean interval 0.3 units.

Usage

data(Simp)

Format

Time

times of observation

X, Y

coordinates

Source

simulated data

Examples

data(Simp)
plot(Simp)

Perform window sweep for BCPA

Description

This is the workhorse function of the BCPA. It performs a sweep of the time series, searching for most significant change points and identifying the parsimonious model according to an adjusted BIC.

Usage

WindowSweep(
  data,
  variable,
  time.var = "T.mid",
  windowsize = 50,
  windowstep = 1,
  units = "hours",
  K = 2,
  tau = TRUE,
  range = 0.6,
  progress = TRUE,
  plotme = FALSE,
  ...
)

Arguments

Value

an object of class windowsweep, which is a list containing:

Author(s)

Eliezer Gurarie

See Also

for internal functions: GetModels, GetBestBreak, GetDoubleL; for summarizing output: ChangePointSummary; for plotting output: plot.bcpa

Examples

# Using the included simulated  movement data
require(bcpa)
data(Simp)
plot(Simp)
Simp.VT <- GetVT(Simp)
## note: column names of Simp.VT include:
###  "S" - step length
### "Theta" - turning angle
###  T.start" "T.end"   "T.mid"  "T.POSIX" - time variables

if(interactive())
  Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, 
                         windowstep = 1, progress=TRUE) else
  Simp.ws <- WindowSweep(Simp.VT, "V*cos(Theta)", windowsize = 50, 
                         windowstep = 1, progress=FALSE)

plot(Simp.ws, threshold = 7)
plot(Simp.ws, type = "flat", clusterwidth = 3)
PathPlot(Simp, Simp.ws)
PathPlot(Simp, Simp.ws, type="flat")
DiagPlot(Simp.ws)

# Example of application on one dimensional data (e.g. depth)

# simulate some data with three change points: surface, medium, deep, surface

## random times within 100 hours of now
n.obs <- 100
time = sort(Sys.time() - runif(n.obs, 0, n.obs) * 60 * 60)

d1 <- 50; d2 <- 100
t1 <- 25; t2 <- 65; t3 <- 85
sd1 <- 1; sd2 <- 5; sd3 <- 10

dtime <- as.numeric(difftime(time, min(time), units="hours"))
phases <- cut(dtime, c(-1, t1, t2, t3, 200), labels = c("P1","P2","P3","P4")) 
means <- c(0,d1,d2,0)[phases]
sds <- c(sd1,sd2,sd3,sd1)[phases]
depth <- rnorm(n.obs,means, sds)
# make all depths positive!
depth <- abs(depth)
mydata <- data.frame(time, depth)

# perform windowsweep
ws <- WindowSweep(mydata, "depth", time.var = "time", windowsize = 30, 
                  windowstep = 1, progress=interactive())
plot(ws)
plot(ws, type = "flat", cluster = 6)
ChangePointSummary(ws, cluster = 6)

Plotting method for BCPA output

Description

Plotting method for the output of a BCPA analysis with vertical break points, superimposed estimates of the partitioned mean and variance estimates and color-coded autocorrelation estimates.

Usage

## S3 method for class 'bcpa'
plot(
  x,
  type = c("smooth", "flat")[1],
  threshold = 3,
  clusterwidth = 1,
  col.cp = rgb(0.5, 0, 0.5, 0.5),
  pt.cex = 0.5,
  legend = TRUE,
  rho.where = "topleft",
  mu.where = "nowhere",
  col.sd = "red",
  col.mean = "black",
  ...
)

Arguments

Author(s)

Eliezer Gurarie

See Also

Plots output of the WindowSweep function.

Examples

if(!exists("Simp.ws"))
{
 data(Simp)
 Simp.ws <- WindowSweep(GetVT(Simp), "V*cos(Theta)", windowsize = 50, 
 windowstep = 1, progress=TRUE)
}

plot(Simp.ws)
# this actually provides basically the exact original changepoints
plot(Simp.ws, threshold=7)
# here's a flat analysis
plot(Simp.ws, type="flat", clusterwidth=3, legend=FALSE)

Plot Track

Description

Default method for plotting tracks

Usage

## S3 method for class 'track'
plot(x, pch = 19, col = rgb(0, 0, 0, 0.2), ...)

Arguments

Examples

data(Simp)
is(Simp)
plot(Simp)