PhysMove includes 5 metrics for quantifying movement patterns that are based on ten functions, including:
rms()calcDisp() and
plotDispPDF()fitDist(), compDist(), and
plotDist()randomise() and
plotRandomTracks()turningAngles() and
plotAngles()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:
timeUnit: time unit used to calculate the time between
locations (timeUnit= "days", by default),wBins: width of the time bins used to calculate how
frequently displacements occurred (wBins=1.1, by
default),plot: create a scatter plot (plot=TRUE, by
default), andlm: fit a linear model to examine the relationship
between root-mean-square displacements and time (lm=TRUE,
by default).rms() results are output as a list. The first list
element is a data frame of results with three columns:
timeUnitThe 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.4973449The 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:
min_hr and max_hr: set the minimum and
maximum times between location estimates in hours, respectively
(min_hr=24 and max_hr=240, by default),interval_hr: set the time interval in hours. This
parameter creates a sequence of time windows between the minimum and
maximum times over the set time interval (interval_hr=24,
by default), andrange_hr: set the range in hours. This parameter allows
the code to identify location estimates that are close to, but not
exactly separated by the interval_hr input value
(range_hr=6, by default).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.8543Probability 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:
normalised: normalise the data before plotting, which
divides all displacements in a time window by the mean displacement for
that time window (normalised=TRUE, by default).colours: change the point colours
(colours=rainbow, by default) andlegend: add or remove a legend
(legend=TRUE, by default).plotDispPDF() outputs a data frame of all data used to
create the plot, including:
normalised=TRUE the displacements are normalised values),
and# 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.
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:
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:
dist: distributions you want to fit to the displacement
data. fitDist() can fit continuous power-law (“pl”),
exponential (“exp”), and lognormal (“lnorm”) distributions
(dist=c("pl","exp","lnorm"), by default). To fit only one
or two distributions simply remove the distribution(s) you are not
interested in running, e.g., dist=c("exp","lnorm").set_dmin: To limit the fitted distribution to values
above a specified value. If your data are going to be normalised this
value will have to be a normalised value as well. Default = NULL.full: determines if distributions are fit over the full
range of displacement data (full=TRUE), or to displacements
truncated by a minimum value (full=FALSE, by default)normalise: normalise displacements before fitting
distributions (normalise=TRUE, by default). Displacements
should be normalised if they were calculated over multiple temporal
periods.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().
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:
fitLines: add fit lines for each distribution
(fitLines=TRUE, by default),setDist: plot only specific distributions
(setDist=NULL, by default, which will plot all
distributions),colours: change the colours of the fit lines
(colours=c("red","gold2","blue"), by default),legend: add a legend (legend=TRUE, by
default), andlabel: X axis label. Note that “Normalised” will
automatically be added if distributions were fit to normalised data.
Default = NULL and will result in x-axis label of “input data”.plotDist() outputs a plot and a data frame of the
displacements (x values) and ccdf values (y values) used in the
plot.
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:
force_AICc: force compDist() to calculate
an AICc (force_AICc=FALSE, by default). By default,
compDist() compares distribution fits using weighted AICc
scores (AIC scores corrected for small sample sizes) when the ratio of
sample size (nTail) to the number of parameters (K) is ≤ 40,
based on the model with the largest K; else, weighted AIC scores are
calculated, following Burnham and Anderson (2004).
force_AICc is used to calculate an AICc instead of an AIC
score (if force_AICc = TRUE). The highest wAIC or wAICc
score from each comparison indicates the best-fit distribution.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.
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.
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-156We 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 distributionOnce 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-31Overall 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.
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:
randTrack: change the number of randomised tracks that
are created (randTrack=100, by default),gridCell: change the grid cell size in degrees
(gridCell=0.25, by default),plot: create a scatter plot of the results
(plot=TRUE, by default), andlm: fit a linear model to the average number of grid
cells visited by the randomised tracks and the number of grid cells
visited by the original tracks (lm=TRUE, by default). The
slope of this model can be used to discuss how correlations may have
influenced interpretations of movement.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.8709073To visualise the tracks created with randomise() you can
use plotRandomTracks() (Figure V10).
plotRandomTracks() requires three parameters:
randomise().plotRandomTracks() also includes 6 optional
parameters:
numPlot: the number of randomised tracks to plot
(numPlot=1:5, by default, which will plot the first 5
randomised versions of each track),colours: the colours of the original and randomised
location estimates, respectively
(colours=c("black","grey70"), by default),tracks: connect points with lines
(tracks=TRUE, by default),startCol and endCol: change the colours of
the starting and ending points of each track, respectively
(startCol="red" and endCol = "blue", by
default), andlegend: add a legend (legend=TRUE, by
default).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.
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:
min_hr and max_hr: set the minimum and
maximum times between location estimates in hours
(min_hr=24 and max_hr=240, by default),interval_hr: set the time interval in hours, which
creates a sequence of time windows between the minimum and maximum times
over a set time interval (interval_hr=24, by default),range_hr: set the range in hours, which allows the code
to identify location estimates that are close to, but not exactly
separated by the interval_hr input value
(range_hr=6, by default), andhistPlot: output a histogram and control if “all” time
windows are plotted or if only a specific time window is plotted, in
which case “all” is replaced with a number corresponding to the desired
time window, e.g., histPlot=c(TRUE,1) will plot the first
time window (histPlot=c(TRUE,"all"), by default).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.99861Results 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:
timePlot: control if “all” time windows or only
specific windows are plotted (timePlot="all", by
default),colours: change line colours
(colours=rainbow, by default), andlegend: add a legend (legend=TRUE, by
default).plotAngles() outputs a data frame of all data used to
create the circle plot, including:
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.
Burnham, K.P. & Anderson, D.R. (2004) Multimodel Inference: Understanding AIC and BIC in Model Selection. Sociological Methods & Research, 33, 261-304.
Calich, H.J. et al. (2021) Comprehensive analytical approaches reveal species-specific search strategies in sympatric apex predatory sharks. Ecography, 44, 1544-1556.
Farage, C. et al. (2021) Identifying flow modules in ecological networks using Infomap. Methods in Ecology and Evolution, 12, 778–786.
Méndez, V., et al. (2013). Stochastic Foundations in Movement Ecology: Anomalous Diffusion, Front Propagation and Random Searches. Berlin, Heidelberg, Germany, Springer Berlin / Heidelberg.
Rodríguez, J.P. et al. (2017) Big data analyses reveal patterns and drivers of the movements of southern elephant seals. Scientific Reports, 7, 1-10.
Viswanathan, G. M., et al. (2011). The Physics of Foraging: An Introduction to Biological Encounters and Random Searches. Cambridge, Cambridge University Press.
Wickham, H. (2016) ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag, New York.