PhysMove: Movement Patterns

Hannah J. Calich, Jorge Rodríguez, Víctor Eguíluz & Ana M. M. Sequeira

Last updated: 2026-07-25

Index

  1. Introduction and data preparation
  2. Movement patterns
  3. Space-use patterns
  4. Intraspecific movements

Movement patterns

PhysMove includes 5 metrics for quantifying movement patterns that are based on ten functions, including:

Scale of movement

The rms() function provides insights into how movement scales with time by calculating mean and root-mean-square (RMS) displacements and plotting them over time (Figure V2).

rms() requires a data frame with telemetry data (see data formatting) and includes four optional parameters:

rms() results are output as a list. The first list element is a data frame of results with three columns:

The second list element is a data frame of results from the linear model and is only exported when lm=TRUE. The slope of the linear model describes how displacement scales with time.

Note that because rms() calculates all displacements in each track, this function can take time. Progress updates will appear when calculations are 25%, 50%, 75%, and 100% complete.

# Calculate RMS values with default parameters
rms.result <- rms(tracks)
#> 25% complete
#> 50% complete
#> 75% complete
#> Calculations complete

#> Scaling exponent = 0.4973

Figure V2 Scatter plot of mean (grey points; q=1) and root-mean-square (RMS; black points; q=2) displacements (d) in kilometers (km) from 'tracks' dataset over time (T) in days, fit to a linear model (red line with standard error shaded in grey). Plot created with rms() default parameters.

# Summarise RMS results
summary(rms.result[["rmsResults"]])
#>    timeWindow        meanDisplacements rmsDisplacements
#>  Min.   :   0.9755   Min.   :  8.23    Min.   : 11.41  
#>  1st Qu.:  16.2477   1st Qu.: 38.37    1st Qu.: 44.46  
#>  Median :  64.6389   Median : 78.68    Median : 90.40  
#>  Mean   : 190.3053   Mean   :111.46    Mean   :122.62  
#>  3rd Qu.: 257.7394   3rd Qu.:154.68    3rd Qu.:174.23  
#>  Max.   :1025.3728   Max.   :583.63    Max.   :584.91

# Summarise linear model results and identify the scaling exponent 
RMSlinearModel <- rms.result[["lm"]]
print(RMSlinearModel)
#>                          term  estimate   std.error statistic      p.value
#> 1                 (Intercept) 2.4177080 0.029749856  81.26789 1.243332e-60
#> 2 log(RMS_Result$timeBin_log) 0.4973449 0.006669422  74.57092 1.599496e-58

# Determine the scaling exponent 
RMSlinearModel$estimate[2]
#> [1] 0.4973449

Back to top

Movement patterns across temporal scales

The calcDisp() function calculates displacements travelled in kilometres over set time windows.

calcDisp() requires a data frame with telemetry data (see data formatting) and has four optional parameters that allow you to change different aspects of the time windows:

For example, by default, calcDisp() calculates displacements between location estimates separated by 10 time windows: 24 ± 6 hours, 48 ± 6 hours, 72 ± 6 hours, etc., until 240 ± 6 hours. Note that displacements smaller than 0.001 km (1 m) are excluded to avoid zero or near-zero values arising from duplicated or stationary locations.

calcDisp() outputs a list where each list element contains the displacements calculated over a time window, such that the first list element contains data from the first time window and so on. For example, by default, the first list element includes displacements calculated over 24 ± 6 hours, and the tenth list element includes displacements calculated over 240 ± 6 hours.

# Calculate displacements with default parameters
dispAll <- calcDisp(tracks)

# [1] "15598 displacements in 24 +/- 6 hour(s)"
# [1] "15573 displacements in 48 +/- 6 hour(s)"
# [1] "15548 displacements in 72 +/- 6 hour(s)"
# [1] "15523 displacements in 96 +/- 6 hour(s)"
# [1] "15498 displacements in 120 +/- 6 hour(s)"
# [1] "15473 displacements in 144 +/- 6 hour(s)"
# [1] "15448 displacements in 168 +/- 6 hour(s)"
# [1] "15423 displacements in 192 +/- 6 hour(s)"
# [1] "15398 displacements in 216 +/- 6 hour(s)"
# [1] "15373 displacements in 240 +/- 6 hour(s)"
# Summarise displacements calculated over the first time window (24 ± 6 hrs)
summary(unlist(dispAll[[1]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.2605  2.7049  5.5149  8.2299 11.1227 76.8543

Probability density function (pdf) of displacements

Probability density functions (pdfs) of displacements describe the probability distribution of the displacements and can be used to calculate the probability of displacements occurring (Figures V3 to Figure V4).

The plotDispPDF() function requires a list of displacements calculated using calcDisp() and includes 3 optional parameters:

plotDispPDF() outputs a data frame of all data used to create the plot, including:

# Create a probability density function (pdf) plot of normalised 
# displacements
plot.data <- plotDispPDF(dispAll)

Figure V3 Probability density function (pdf) plot of normalised displacements from the 'tracks' dataset calculated over 10 time windows, 24 to 240 hours at 24 ± 6-hour time intervals with calcDisp(). Plot created with plotDispPDF() default parameters (i.e., where normalised=TRUE).

# Create a probability density function (pdf) plot of raw (i.e., not 
# normalised) displacements
plot.data.norm <- plotDispPDF(dispAll, normalised=FALSE)

Figure V4 Probability density function (pdf) plot of displacements from the 'tracks' dataset calculated over 10 time windows, 24 to 240 hours at 24 ± 6-hour time intervals with calcDisp(). Plot created with plotDispPDF() where normalised=FALSE.

Back to top

Search patterns

PhysMove can be used to identify the best-fit distribution of displacements, which can provide insights into the search pattern(s) a species may use to locate resources. Determining the best-fit distribution for the displacements involves 3 functions:

  1. fitDist():

Fits cdfs of continuous power-law, exponential, and lognormal distributions over the full range of displacements (i.e., full distributions) or to displacements truncated by a minimum value (i.e., truncated distributions). The fitDist() function requires a list of values (e.g., displacements calculated using calcDisp()) and includes four optional parameters:

fitDist() outputs a list including two list elements. The first list element is a data frame that includes:

The second list element is a logical argument that records if the displacements were normalised (TRUE) or not (FALSE). This information is required for plotDist() and compDist().

  1. plotDist():

Uses the results from fitDist() to plot ccdfs of the displacements with fit lines for each distribution. The plotDist() function requires a list of values (e.g., calculated using calcDisp()) and results from fitDist(), and includes four optional parameters:

plotDist() outputs a plot and a data frame of the displacements (x values) and ccdf values (y values) used in the plot.

  1. compDist():

Compares distribution fits from fitDist() and identifies the best-fit distribution for the displacements. Note that compDist() can only be used when all distributions are fit to the same range of data (e.g., when full=TRUE or if set_dmin≠NULL). See Figure V5 for a methods overview. The compDist() function requires a list of displacements calculated using calcDisp(), the results from fitDist(), and includes one optional parameter:

compDist() outputs a data frame that contains the summary statistics for each distribution fit (as described for fitDist()) with the corresponding AICc/AIC scores and weighted AICc/AIC scores (wAICc/wAIC). The distribution with the highest wAIC or wAICc score from each comparison is the best-fit distribution.

Figure V5 Diagram outlining the procedure for identifying the best-fit distribution of displacements.

Search patterns example

In the example below we calculate displacements over 24 ± 6 hours, plot a pdf of the displacements, and identify the best-fit distributions for both full and truncated datasets. Using a single time window helps keep run times manageable, while still capturing the overall displacement patterns. In general, we recommend fitting distributions to both full and truncated datasets to gain a comprehensive understanding of displacement patterns.

We begin by calculating displacements over 24 ± 6 hours with calcDisp() and plotting a pdf of the displacements with plotDispPDF() (Figure V6).

# Calculate displacements over 24 ± 6 hours
disp <- calcDisp(tracks, max_hr=24)

# Summarise displacements
summary(unlist(disp))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.2605  2.7049  5.5149  8.2299 11.1227 76.8543

# Plot displacements (as displacements were only calculated over one time window they do not need to be normalised)
plot.data.pdf <- plotDispPDF(disp, normalised=FALSE)

Figure V6 Probability density function (pdf) plot of displacements calculated using calcDisp() with max_hr=24. Plot created with plotDispPDF() and normalised=FALSE.

Fitting full distributions

We use fitDist() to fit the full range of distributions calculated over 24 ± 6 hours to power-law, exponential, and lognormal distributions, and plotDist() to visualise the results (Figure V7).

# Fit all distributions to the full range of displacement data 
distResults <- fitDist(disp, full=TRUE, normalise=FALSE) 
#> Fitting a power law distribution
#> Fitting an exponential distribution
#> Fitting a lognormal distribution

distResults[["distResults"]]
#>   distribution      dmin parameter1 parameter2 nTail
#> 1           pl 0.2605242  1.3341746         NA 15598
#> 2          exp 0.2605242  0.1254743         NA 15598
#> 3        lnorm 0.2605242  1.6405463   1.044668 15598
# Create a ccdf plot of displacements with fit lines illustrating 
# distributions fit to the full range of displacements
plot.data.all.pdf <- plotDist(disp, distResults, label="Displacements (km)")

Figure V7 Complementary cumulative distribution function (ccdf) of displacements (calculated using calcDisp() with max_hr=24). Plot includes fit lines for power-law (pl), exponential (exp), and lognormal (lnorm) distributions based on results from fitDist() with full=TRUE. Plot created using plotDist() default parameters.

Next, we use compDist() to compare distribution fits over the full range of displacement data.

# Identify the best-fit distribution for the full range of displacement data
compResults <- compDist(disp, distResults)
compResults
#>   distribution      dmin parameter1 parameter2 nTail       AIC          wAIC
#> 1           pl 0.2605242  1.3341746         NA 15598 116783.60  0.000000e+00
#> 2          exp 0.2605242  0.1254743         NA 15598  95948.77  1.000000e+00
#> 3        lnorm 0.2605242  1.6405463   1.044668 15598  96664.80 3.274395e-156

Fitting truncated distributions

We use fitDist() to identify the best-fit dmin for each distribution and fit truncated power-law, exponential, and lognormal distributions to displacements calculated over 24 ± 6 hours. We visualise results with plotDist() (Figure V8).

# Fit all distributions and identify the best-fit dmin for each distribution
distResultsTrunc <- fitDist(disp, full=FALSE, normalise=FALSE)
print(distResultsTrunc[["distResults"]])
#>   distribution      dmin parameter1 parameter2 nTail
#> 1           pl 20.874018  4.8875853         NA  1387
#> 2          exp  5.771445  0.1207768         NA  7484
#> 3        lnorm  1.968525  1.9017072   0.829428 12657
# Create a ccdf plot of displacements with fit lines illustrating 
# distributions fit to the best-fit dmin for each distribution
plot.data.all.trunc <- plotDist(disp, distResultsTrunc, label="Displacements (km)")

Figure V8 Complementary cumulative distribution function (ccdf) of displacements calculated using calcDisp() with max_hr=24 including fit lines for truncated power-law (pl), exponential (exp), and lognormal (lnorm) distributions based on the best-fit dmin results from fitDist(). Plot created using plotDist() default parameters.

Since each distribution was fit to a different range of data (i.e., nTail values are different for each distribution), we cannot run compDist() directly. Instead, we must make pairwise comparisons where fitDist() is re-run three times (once for each of the three distributions) where the set_dmin parameter is set to each of the best-fit dmin values in turn.

# Fit all distributions using the dmin value for the 
# power-law distribution
dmin <- distResultsTrunc[["distResults"]][1,2]
distResultsPl <- fitDist(disp, set_dmin=dmin, normalise=FALSE)
#> Fitting a power law distribution
#> Fitting an exponential distribution
#> Fitting a lognormal distribution
# Fit all distributions using the dmin value for the 
# exponential distribution
dmin <- distResultsTrunc[["distResults"]][2,2]
distResultsExp <- fitDist(disp, set_dmin=dmin, normalise=FALSE)
#> Fitting a power law distribution
#> Fitting an exponential distribution
#> Fitting a lognormal distribution
# Fit all distributions using the dmin value for the 
# lognormal distribution
dmin <- distResultsTrunc[["distResults"]][3,2]
distResultsLnorm <- fitDist(disp, set_dmin=dmin, normalise=FALSE)
#> Fitting a power law distribution
#> Fitting an exponential distribution
#> Fitting a lognormal distribution

Once all distributions are fit using each of the best-fit dmin values, the distribution fits can be compared using compDist(). An important consideration for interpreting the compDist() results from pairwise comparisons is that if a dmin was set to favour a specific distribution, but the wAIC (or wAICc) scores do not identify that distribution as the best fit, the distribution corresponding to the dmin value is not the best-fit distribution for the displacements. For example, in the first pairwise comparison below, you’ll see the dmin was set to the best-fit dmin for a power-law (20.87); however, the wAIC scores identified an exponential distribution as the best fit. Therefore, we conclude that a power-law is not the best-fit distribution for the data.

# Compare distribution fits based on the best-fit dmin value for the power-law distribution
compResultsPl <- compDist(disp, distResultsPl)
compResultsPl
#>   distribution     dmin parameter1 parameter2 nTail      AIC         wAIC
#> 1           pl 20.87402  4.8875857         NA  1387 8016.717 1.164441e-08
#> 2          exp 20.87402  0.1457233         NA  1387 7980.295 9.440862e-01
#> 3        lnorm 20.87402  2.4000922  0.5294741  1387 7985.948 5.591380e-02
# Compare distribution fits based on the best-fit dmin value for the exponential distribution
compResultsExp <- compDist(disp, distResultsExp)
compResultsExp
#>   distribution     dmin parameter1 parameter2 nTail      AIC         wAIC
#> 1           pl 5.771445  2.3168676         NA  7484 48305.18 0.000000e+00
#> 2          exp 5.771445  0.1207767         NA  7484 46464.64 1.000000e+00
#> 3        lnorm 5.771445  2.2645346  0.6574655  7484 46545.40 2.906385e-18
# Compare distribution fits based on the best-fit dmin value for the lognormal distribution
compResultsLnorm <- compDist(disp, distResultsLnorm)
compResultsLnorm
#>   distribution     dmin parameter1 parameter2 nTail      AIC         wAIC
#> 1           pl 1.968525  1.7440593         NA 12657 83839.95 0.000000e+00
#> 2          exp 1.968525  0.1263384         NA 12657 77567.29 1.000000e+00
#> 3        lnorm 1.968525  1.9017083  0.8294269 12657 77708.50 2.172308e-31

Overall conclusion: The data were best-fit to both full and truncated exponential distributions, when compared with power-law and lognormal distributions. Because the full distribution includes all of the data we will refer to the full distribution when presenting and discussing our final results.

Back to top

Influence of correlations on movement decisions

The randomise() function can be used to gain insights into how correlations influenced a species’ movements and space-use (Figure V9).

randomise() requires a data frame with telemetry data (see data formatting) and includes four optional parameters:

randomise() outputs a list with three list elements. The first list element is a data frame with three columns:

The second and third list elements contain the randomised longitude and latitude values, respectively, which are needed for the plotRandomTracks() function

# randomise() involves random number selection, so setting a seed enables the replication of results
set.seed(1)

# Randomise tracks from the 'tracks' dataset with default parameters
randomise.result <- randomise(tracks)
#> Randomizing tracks, step 1/3
#> Calculating average number of cells visited by randomised tracks, step 2/3
#> Calculating number of cells visited by original tracks, step 3/3

#> Slope = 0.8709

Figure V9 Scatter plot illustrating the relationship between the number of grid cells visited by the original tracks from the 'tracks' dataset and the average number of grid cells visited by the randomised tracks. The solid black line represents the linear model fit to this data, the grey shaded area reflects the standard error of the fit, and the dashed black line represents a 1:1 relationship. Plot created with randomise() default parameters.

# Summarise RMS results
summary(randomise.result[["resultsDF"]])
#>       ref     CellsInOriginalTracks AvgCellsInRandomisedTracks
#>  Min.   : 1   Min.   :27.00         Min.   :25.29             
#>  1st Qu.: 7   1st Qu.:39.00         1st Qu.:40.50             
#>  Median :13   Median :59.00         Median :54.93             
#>  Mean   :13   Mean   :57.32         Mean   :56.50             
#>  3rd Qu.:19   3rd Qu.:71.00         3rd Qu.:69.90             
#>  Max.   :25   Max.   :94.00         Max.   :99.58

# Determine the slope of the linear model
RandomiselinearModel <- randomise.result[["lm"]]
print(RandomiselinearModel)
#>                            term  estimate std.error statistic      p.value
#> 1                   (Intercept) 6.5763920 6.7960283 0.9676817 3.432696e-01
#> 2 plot.df$CellsInOriginalTracks 0.8709073 0.1122742 7.7569703 7.251785e-08

# Determine the slope without displaying the full linear model summary
RandomiselinearModel$estimate[2]
#> [1] 0.8709073

Plot randomised tracks

To visualise the tracks created with randomise() you can use plotRandomTracks() (Figure V10).

plotRandomTracks() requires three parameters:

  1. data frame with telemetry data (see data formatting),
  2. reference ID of the track you want to map (ref must be included in the telemetry data frame), and
  3. results from randomise().

plotRandomTracks() also includes 6 optional parameters:

plotRandomTracks() outputs the data used to create the map in three columns:

# Plot random tracks for 'tracks' dataset reference ID 1
plot.data.random.tracks <- plotRandomTracks(tracks, ref=1, randomise.result)

Figure V10 Map illustrating the original track for reference ID 1 from the 'tracks' dataset (black points and line) and the first 5 randomised tracks for track reference ID 1 calculated using randomise() (grey points and lines). The starting and ending locations are in red and blue, respectively. Plot created with plotRandomTracks() default parameters and ref=1.

Back to top

Turning angles

The turningAngles() function calculates turning angles between sets of three consecutive location estimates separated by set time windows to describe how species explore their habitats (Figure V11).

turningAngles() requires a data frame with telemetry data (see data formatting) and includes 5 optional parameters:

Results are output in a list where each list element contains the angles calculated over a time window, such that the first list element contains data from the first time window and so on.

# Calculate turning angles in the 'tracks' dataset using default parameters
angleListAll <- turningAngles(tracks)

# [1] "15573 angles in 24 +/- 6 hour(s)"
# [1] "15523 angles in 48 +/- 6 hour(s)"
# [1] "15473 angles in 72 +/- 6 hour(s)"
# [1] "15423 angles in 96 +/- 6 hour(s)"
# [1] "15373 angles in 120 +/- 6 hour(s)"
# [1] "15323 angles in 144 +/- 6 hour(s)"
# [1] "15273 angles in 168 +/- 6 hour(s)"
# [1] "15223 angles in 192 +/- 6 hour(s)"
# [1] "15173 angles in 216 +/- 6 hour(s)"
# [1] "15123 angles in 240 +/- 6 hour(s)"

Figure V11 Histogram of turning angles from the 'tracks' dataset over ten time windows (24 to 240 hours at 24 ± 6 hour intervals). Plot created with turningAngles() default parameters.

# Summarise turning angles calculated over the first time window 
summary(angleListAll[[1]])
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
#> -179.97980  -93.41046    0.79527   -0.01537   92.75782  179.99861

Create a circle plot

Results from turningAngles() can be visualised with the plotAngles() function, which creates a circle plot (also known as a spider or radar plot) showing the frequency of turning angles over each time window (Figure V12).

plotAngles() requires the list of angles output from turningAngles() and includes 3 optional parameters:

plotAngles() outputs a data frame of all data used to create the circle plot, including:

# Plot angles with a circle plot
plot.data.angles <- plotAngles(angleListAll)

Figure V12 Circle plot of turning angles recorded from the 'tracks' dataset during ten time windows (24 to 240 hours at 24 ± 6 hour intervals). Plot created with plotAngles() default parameters.

Proceed to Space-Use Patterns

Back to top