distantia: an R package to compute the dissimilarity between

multivariate time-series

Summary

The package distantia allows to measure the dissimilarity betweenmultivariate ecological time-series (METS hereafter). The packageassumes that the target sequences are ordered along a given dimension,being depth and time the most common ones, but others such as latitudeor elevation are also possible. Furthermore, the target METS can beregular or irregular, and have their samples aligned (sameage/time/depth) or unaligned (different age/time/depth). The onlyrequirement is that the sequences must have at least two (but ideallymore) columns with the same name and units representing differentvariables relevant to the dynamics of an ecological system.

In this document I explain the logics behind the method, show how to useit, and demonstrate how the distantia package introduces usefultools to compare multivariate time-series. The topics covered in thisdocument are:

• Installation of the package.
• Comparing two irregular METS.
• Comparing multiple irregular METS.
• Comparing regular and aligned METS.
• Restricted permutation test to assess the significance ofdissimilarity scores
• Assessing the contribution of every variable to the dissimilaritybetween two METS.
• Partial matching: finding the section in a long METS more similar toa given shorter one.
• Sequence slotting: combining samples of two METS into a singlecomposite sequence.
• Tranferring an attribute from one METS to another: direct andinterpolated modes.

Installation

You can install the released version of distantia (currently v1.0.0)from CRAN with:

rinstall.packages("distantia")

And the development version (currently v1.0.1) fromGitHub with:

rlibrary(devtools)devtools::install_github("BlasBenito/distantia")

Loading the library, plus other helper libraries:

rlibrary(distantia)library(ggplot2)library(viridis)library(kableExtra)library(qgraph)library(tidyr)

Comparing two irregular METS

In this section I will use two example datasets based on the Abernethypollen core (Birks and Mathewes, 1978) to fully explain the logicalbackbone of the dissimilarity analyses implemented in distantia.

``` r

loading sequences

data(sequenceA)data(sequenceB)

showing first rows

kable(sequenceA[1:15, ], caption = "Sequence A")```

rkable(sequenceB[1:15, ], caption = "Sequence B")

Data preparation

Notice that sequenceB has a few NA values (that were introduced toserve as an example). The function prepareSequences gets them readyfor analysis by matching colum names and handling empty data. It allowsto merge two or more METS into a single dataframe ready for furtheranalyses. Note that, since the data represents pollen abundances, aHellinger transformation (square root of the relative proportions ofeach taxa, it balances the relative abundances of rare and dominanttaxa) is applied. This transformation balances the relative importanceof very abundant versus rare taxa. The function prepareSequenceswill generally be the starting point of any analysis performed with thedistantia package.

``` r

checking the function help-file.

help(prepareSequences)

preparing sequences

AB.sequences <- prepareSequences( sequence.A = sequenceA, sequence.A.name = "A", sequence.B = sequenceB, sequence.B.name = "B", merge.mode = "complete", if.empty.cases = "zero", transformation = "hellinger")

showing first rows of the transformed data

kable(AB.sequences[1:15, ], digits = 4, caption = "Sequences A and B ready for analysis.")```

The computation of dissimilarity between the datasets A and B requiresseveralsteps.

1. Computation of a distance matrix among the samples of both sequences.

It is computed by the distanceMatrix function, which allows the userto select a distance metric (so far the ones implemented aremanhattan, euclidean, chi, and hellinger). The functionplotMatrix allows an easy visualization of the resulting distancematrix.

``` r

computing distance matrix

AB.distance.matrix <- distanceMatrix( sequences = AB.sequences, method = "euclidean")

plotting distance matrix

plotMatrix( distance.matrix = AB.distance.matrix, color.palette = "viridis", margins = rep(4,4))```

2. Computation of the least-cost path within the distance matrix.

This step uses a dynamic programming algorithm to find the least-costpath that connnects the cell 1,1 of the matrix (lower left in the imageabove) and the last cell of the matrix (opposite corner). This can bedone via in two different ways.

• an orthogonal search by moving either one step on the x axisor one step on the y axis at a time (see Equation 1).

Equation 1

• an orthogonal and diagonal search (a.k.a diagonal) whichincludes the above, plus a diagonal search.

Equation 2
![AB_{between} = 2 \times (D(A_{1}, B_{1}) +\sum_{i=1}^{m}\sum_{j=1}^{n} min\left(\begin{array}{c}D(A_{i},B_{j+1}), \\ D(A_{i+1}, B_{j} \\ D(A_{i+1}, B_{j+1})\end{array}\right))](https://latex.codecogs.com/png.latex?AB%7Bbetween%7D%20%3D%202%20%5Ctimes%20%28D%28A%7B1%7D%2C%20B%7B1%7D%29%20%2B%20%5Csum%7Bi%3D1%7D%5E%7Bm%7D%5Csum%7Bj%3D1%7D%5E%7Bn%7D%20min%5Cleft%28%5Cbegin%7Barray%7D%7Bc%7DD%28A%7Bi%7D%2C%20B%7Bj%2B1%7D%29%2C%20%5C%5C%20D%28A%7Bi%2B1%7D%2C%20B%7Bj%7D%20%5C%5C%20D%28A%7Bi%2B1%7D%2C%20B%7Bj%2B1%7D%29%20%5Cend%7Barray%7D%5Cright%29%29"AB{between} = 2 \times (D(A{1}, B{1}) + \sum{i=1}^{m}\sum{j=1}^{n} min\left(\begin{array}{c}D(A{i}, B{j+1}), \\ D(A{i+1}, B{j} \\ D(A{i+1}, B{j+1}) \end{array}\right))")

Where:

• and are the number ofsamples of the multivariate time-series and.
• and are the indices ofthe samples in and being considered oneach step of the recursive algorithm.
• is a function thatreturns the multivariate distance (i.e. Manhattan) between any givenpair of samples of and .
• is a functionreturning the minimum distance of a subset of distances defined inthe neigborhood of the given indices and

The equation returns, which is the double of the sum of distances that liewithin the least-cost path, and represent the distance between thesamples of A and B. The value of is computed by using the functions leastCostMatrix,which computes the partial solutions to the least-cost problem,leastCostPath, which returns the best global solution, andleastCost function, which sums the distances of the least-cost pathand multiplies them by 2.

The code below performs these steps according to both equations

``` r

ORTHOGONAL SEARCH

computing least-cost matrix

AB.least.cost.matrix <- leastCostMatrix( distance.matrix = AB.distance.matrix, diagonal = FALSE)

extracting least-cost path

AB.least.cost.path <- leastCostPath( distance.matrix = AB.distance.matrix, least.cost.matrix = AB.least.cost.matrix, diagonal = FALSE )

DIAGONAL SEARCH

computing least-cost matrix

AB.least.cost.matrix.diag <- leastCostMatrix( distance.matrix = AB.distance.matrix, diagonal = TRUE)

extracting least-cost path

AB.least.cost.path.diag <- leastCostPath( distance.matrix = AB.distance.matrix, least.cost.matrix = AB.least.cost.matrix.diag, diagonal = TRUE )

plotting solutions

plotMatrix( distance.matrix = list( 'A|B' = AB.least.cost.matrix[[1]], 'A|B' = AB.least.cost.matrix.diag[[1]] ), least.cost.path = list( 'A|B' = AB.least.cost.path[[1]], 'A|B' = AB.least.cost.path.diag[[1]] ), color.palette = "viridis", margin = rep(4,4), plot.rows = 1, plot.columns = 2 )```

Computing from these solutions is straightforward with thefunction leastCost

``` r

orthogonal solution

AB.between <- leastCost( least.cost.path = AB.least.cost.path )

diagonal solution

AB.between.diag <- leastCost( least.cost.path = AB.least.cost.path.diag )```

Which returns a value for of 33.7206 for the orthogonal solution, and 22.7596 forthe diagonal one. Diagonal solutions always yield lower values for than orthogonal ones.

Notice the straight vertical and horizontal lines that show up in someregions of the least cost paths shown in the figure above. These areblocks, and happen in dissimilar sections of the compared sequences.Also, an unbalanced number of rows in the compared sequences cangenerate long blocks. Blocks inflate the value of because the distance to a given sample is countedseveral times per block. This problem often leads to false negatives,that is, to the conclusion that two sequences are statisticallydifferent when actually they are not.

This package includes an algorithm to remove blocks from the least costpath, which offers more realistic values for. The function leastCostPathNoBlocks reads a leastcost path, and removes all blocks as follows.

``` r

ORTHOGONAL SOLUTION

removing blocks from least cost path

AB.least.cost.path.nb <- leastCostPathNoBlocks( least.cost.path = AB.least.cost.path )

computing AB.between again

AB.between.nb <- leastCost( least.cost.path = AB.least.cost.path.nb )

DIAGONAL SOLUTION

removing blocks

AB.least.cost.path.diag.nb <- leastCostPathNoBlocks( least.cost.path = AB.least.cost.path.diag )

diagonal solution without blocks

AB.between.diag.nb <- leastCost( least.cost.path = AB.least.cost.path.diag.nb )```

Which now yields 11.2975 for the orthogonal solution, and 16.8667 forthe diagonal one. Notice how now the diagonal solution has a highervalue, because by default, the diagonal method generates less blocks.That is why each measure of dissimilarity (orthogonal, diagonal,orthogonal no-blocks, and diagonal no-blocks) lies within a differentcomparative framework, and therefore, outputs from different methodsshould not be compared.

Hereafter only the diagonal no-blocks option will be considered in theexample cases, since it is the most general and safe solution of thefour mentioned above.

``` r

changing names of the selected solutions

AB.least.cost.path <- AB.least.cost.path.diag.nbAB.between <- AB.between.diag.nb

removing unneeded objects

rm(AB.between.diag, AB.between.diag.nb, AB.between.nb, AB.distance.matrix, AB.least.cost.matrix, AB.least.cost.matrix.diag, AB.least.cost.path.diag, AB.least.cost.path.diag.nb, AB.least.cost.path.nb, sequenceA, sequenceB)```

3. Autosum, or sum of the distances among adjacent samples on each sequence.

This step requires to compute the distances between adjacent samples ineach sequence and sum them, as shown in Equation 3.

Equation 3

This operation is performed by the autoSum function shown below.

``` rAB.within <- autoSum( sequences = AB.sequences, least.cost.path = AB.least.cost.path, method = "euclidean" )AB.within

> $A|B

> [1] 19.69168

```

4. Compute dissimilarity score .

The dissimilarity measure was firstdescribed in the book “Numerical methods in Quaternary pollenanalysis”(Birks and Gordon, 1985). Psi is computed as shown in Equation 4a:

Equation 4a

This equation has a particularity. Imagine two identical sequences A andB, with three samples each. In this case, is computed as

Since the samples of each sequence with the same index are identical,this can be reduced to

which in turn equals as shown in Equation 4, yielding a value of0.

This equality does not work in the same way when the least-cost pathsearch-method includes diagonals. When the sequenes are identical,diagonal methods yield an of 0, leading to a equal to-1. To fix this shift, this package uses Equation 4b instead when is selected, which adds 1 to the final solution.

Equation 4b

In any case, the psi function only requires the least-cost, and theautosum of both sequences to compute. Since weare working with a diagonal search, 1 has to be added to the finalsolution.

rAB.psi <- psi( least.cost = AB.between, autosum = AB.within )AB.psi[[1]] <- AB.psi[[1]] + 1

Which yields a psi equal to 1.7131. The output of psi is a list,that can be transformed to a dataframe or a matrix by using theformatPsi function.

``` r

to dataframe

AB.psi.dataframe <- formatPsi( psi.values = AB.psi, to = "dataframe")kable(AB.psi.dataframe, digits = 4)```

workflowPsi: doing it all at once

All the steps required to compute psi, including the format optionsprovided by formatPsi are wrapped together in the functionworkflowPsi. It includes options to switch to a diagonal method, andto ignore blocks, as shown below.

``` r

checking the help file

help(workflowPsi)

computing psi for A and B

AB.psi <- workflowPsi( sequences = AB.sequences, grouping.column = "id", method = "euclidean", format = "list", diagonal = TRUE, ignore.blocks = TRUE)AB.psi

> $A|B

> [1] 1.713075

```

The function allows to exclude particular columns from the analysis(argument exclude.columns), select different distance metrics(argument method), use diagonals to find the least-cost path (argumentdiagonal), or measure psi by ignoring blocks in the least-cost path(argument ignore.blocks). Since we have observed several blocks in theleast-cost path, below we compute psi by ignoring them.

``` r

cleaning workspace

rm(list = ls())```

Comparing multiple irregular METS

The package can work seamlessly with any given number of sequences, aslong as there is memory enough available (but check the new functionworkflowPsiHP, it can work with up to 40k sequences, if you have acluster at hand, and a few years to waste). To do so, almost everyfunction uses the packages“doParallel” and“foreach”, that togetherallow to parallelize the execution of the distantia functions by usingall the processors in your machine but one.

The example dataset sequencesMIS contains 12 sections of the samesequence belonging to different marine isotopic stages identified by acolumn named “MIS”. MIS stages with odd numbers are generallyinterpreted as warm periods (interglacials), while the odd ones areinterpreted as cold periods (glacials). In any case, this interpretationis not important to illustrate this capability of the library.

rdata(sequencesMIS)kable(head(sequencesMIS, n=15), digits = 4, caption = "Header of the sequencesMIS dataset.")

``` runique(sequencesMIS$MIS)

> [1] "MIS-1" "MIS-2" "MIS-3" "MIS-4" "MIS-5" "MIS-6" "MIS-7"

> [8] "MIS-8" "MIS-9" "MIS-10" "MIS-11" "MIS-12"

```

The dataset is checked and prepared with prepareSequences.

rMIS.sequences <- prepareSequences( sequences = sequencesMIS, grouping.column = "MIS", if.empty.cases = "zero", transformation = "hellinger")

The dissimilarity measure psi can be computed for every combinationof sequences through the function workflowPsi shown below.

``` rMIS.psi <- workflowPsi( sequences = MIS.sequences, grouping.column = "MIS", method = "euclidean", diagonal = TRUE, ignore.blocks = TRUE)

there is also a "high-performance" (HP) version of this function with a much lower memory footprint. It uses the options method = "euclidean", diagonal = TRUE, and ignore.blocks = TRUE by default

MIS.psi <- workflowPsiHP( sequences = MIS.sequences, grouping.column = "MIS")

ordered with lower psi on top

kable(MIS.psi[order(MIS.psi$psi), ], digits = 4, caption = "Psi values between pairs of MIS periods.")```

A dataframe like this can be transformed into a matrix to be plotted asan adjacency network with the qgraph package.

``` r

psi values to matrix

MIS.psi.matrix <- formatPsi( psi.values = MIS.psi, to = "matrix")

dissimilariy to distance

MIS.distance <- 1/MIS.psi.matrix**4

plotting network

qgraph::qgraph( MIS.distance, layout='spring', vsize=5, labels = colnames(MIS.distance), colors = viridis::viridis(2, begin = 0.3, end = 0.8, alpha = 0.5, direction = -1) )```

Or as a matrix with ggplot2.

``` r

ordering factors to get a triangular matrix

MIS.psi$A <- factor(MIS.psi$A, levels=unique(sequencesMIS$MIS))MIS.psi$B <- factor(MIS.psi$B, levels=unique(sequencesMIS$MIS))

plotting matrix

ggplot(data=na.omit(MIS.psi), aes(x=A, y=B, size=psi, color=psi)) + geompoint() + viridis::scalecolor_viridis(direction = -1) + guides(size = FALSE)```

The dataframe of dissimilarities between pairs of sequences can be alsoused to analyze the drivers of dissimilarity. To do so, attributes suchas differences in time (when sequences represent different times) ordistance (when sequences represent different sites) between sequences,or differences between physical/climatic attributes between sequencessuch as topography or climate can be added to the table, so models suchas (were A, B, and C are these attributes) can befitted.

``` r

cleaning workspace

rm(list = ls())```

Comparing regular aligned METS

The package distantia is also useful to compare synchronic sequencesthat have the same number of samples. In this particular case, distancesto obtain are computed only between samples with the sametime/depth/order, and no distance matrix (nor least-cost analysis) isrequired. When the argument paired.samples in prepareSequences isset to TRUE, the function checks if the sequences have the same numberof rows, and, if time.column is provided, it selects the samples thathave valid time/depth columns for every sequence in the dataset.

Here we test these ideas with the climate dataset included in thelibrary. It represents simulated palaeoclimate over 200 ky. at foursites identified by the column sequenceId. Note that this time thetransformation applied is “scaled”, which uses the scale function ofR base to center and scale the data.

``` r

loading sample data

data(climate)

preparing sequences

climate <- prepareSequences( sequences = climate, grouping.column = "sequenceId", time.column = "time", paired.samples = TRUE, transformation = "scale" )```

In this case, the argument paired.samples of workflowPsi must beset to TRUE. Additionally, if the argument same.time is set to TRUE,the time/age of the samples is checked, and samples without the sametime/age are removed from the analysis.

``` r

computing psi

climate.psi <- workflowPsi( sequences = climate, grouping.column = "sequenceId", time.column = "time", method = "euclidean", paired.samples = TRUE, #this bit is important same.time = TRUE, #removes samples with unequal time format = "dataframe")

ordered with lower psi on top

kable(climate.psi[order(climate.psi$psi), ], digits = 4, row.names = FALSE, caption = "Psi values between pairs of sequences in the 'climate' dataset.")```

``` r

cleaning workspace

rm(list = ls())```

Restricted permutation test to assess the significance of dissimilarity values

One question that may arise when comparing time series is “to whatextent are dissimilarity values a result of chance?”. Answering thisquestion requires to compare a given dissimilarity value with adistribution of dissimilarity values resulting from chance. However… howdo we simulate chance in a multivariate time-series? The natural answeris “permutation”. Since samples in a multivariate time-series areordered, randomly re-shuffling samples is out of the question, becausethat would destroy the structure of the data. A more gentler alternativeis to randomly switch single data-points (a case of a variable)independently by variable. This kind of permutation is named “restrictedpermutation”, and preserves global trends within the data, but changeslocal structure.

A restricted permutation test on psi values requires the followingsteps:

• Compute the real psi on two given sequences A and B.
• Repeat the following steps several times (99 to 999):
  • For each case of each column of A and B, randomly apply one ofthese actions:
    • Leave it as is.
    • Replace it with the previous case.
    • Replace it with the next case.
  • Compute randomized psi between A and B and store the value.
• Add real psi to the pool of randomized psi.
• Compute the proportion of randomized psi that is equal or lowerthan real psi.

Such a proportion represents the probability of obtaining a value lowerthan real psi by chance.

Since the restricted permutation only happens at a local scale withineach column of each sequence, the probability values returned are veryconservative and shouldn’t be interpreted in the same way p-values areinterpreted.

The process described above has been implemented in theworkflowNullPsi function. We will apply it to three groups of thesequencesMIS dataset.

``` r

getting example data

data(sequencesMIS)

working with 3 groups (to make this fast)

sequencesMIS <- sequencesMIS[sequencesMIS$MIS %in% c("MIS-4", "MIS-5", "MIS-6"),]

preparing sequences

sequencesMIS <- prepareSequences( sequences = sequencesMIS, grouping.column = "MIS", transformation = "hellinger")```

The computation of the null psi values goes as follows:

``` rrandom.psi <- workflowNullPsi( sequences = sequencesMIS, grouping.column = "MIS", method = "euclidean", diagonal = TRUE, ignore.blocks = TRUE, repetitions = 99, #recommended value: 999 parallel.execution = TRUE)

there is also a high-performance version of this function with fewer options (diagonal = TRUE, ignore.blocks = TRUE, and method = "euclidean" are used by default)

random.psi <- workflowNullPsiHP( sequences = sequencesMIS, grouping.column = "MIS", repetitions = 99, #recommended value: 999 parallel.execution = TRUE)```

Note that the number of repetitions has been set to 9 in order tospeed-up execution. The actual number would ideally be 999.

The output is a list with two dataframes, psi and p.

The dataframe psi contains the real and random psi values. Thecolumn psi contains the dissimilarity between the sequences in thecolumns A and B. The columns r1 to r9 contain the psi valuesobtained from permutations of thesequences.

rkable(random.psi$psi, digits = 4, caption = "True and null psi values generated by workflowNullPsi.")

The dataframe p contains the probability of obtaining the real psivalue by chance for each combination ofsequences.

rkable(random.psi$p, caption = "Probability of obtaining a given set of psi values by chance.")

``` r

cleaning workspace

rm(list = ls())```

Assessing the contribution of a variable to the dissimilarity between two sequences

What variables are more important in explaining the dissimilaritybetween two sequences?, or in other words, what variables contributethe most to the dissimilarity between two sequences? One reasonableanswer is: the one that reduces dissimilarity the most when removed fromthe data. This section explains how to use the functionworkflowImportance follows such a principle to evaluate theimportance of given variables in explaining differences betweensequences.

First, we prepare the data. It is again sequencesMIS, but with onlythree groups selected (MIS 4 to 6) to simplify the analysis.

``` r

getting example data

data(sequencesMIS)

getting three groups only to simplify

sequencesMIS <- sequencesMIS[sequencesMIS$MIS %in% c("MIS-4", "MIS-5", "MIS-6"),]

preparing sequences

sequences <- prepareSequences( sequences = sequencesMIS, grouping.column = "MIS", merge.mode = "complete")```

The workflow function is pretty similar to the ones explained above.However, unlike the other functions in the package, that parallelizeacross the comparison of pairs of sequences, this one parallelizes thecomputation of psi on combinations of columns, removing one columneach time.

WARNING: the argument ‘exclude.columns’ of ‘workflowImportance’ doesnot work in version 1.0.0 (available in CRAN), but the bug is fixed inversion 1.0.1 (available in GitHub). If you are using 1.0.0, I recommendyou to subset ‘sequences’ so only the grouping column and the numericcolumns to be compared are available for the function.

``` rpsi.importance <- workflowImportance( sequences = sequencesMIS, grouping.column = "MIS", method = "euclidean", diagonal = TRUE, ignore.blocks = TRUE )

there is also a high performance version of this function, but with fewer options (it uses euclidean, diagonal, and ignores blocks by default)

psi.importance <- workflowImportanceHP( sequences = sequencesMIS, grouping.column = "MIS" )```

The output is a list with two slots named psi and psi.drop.

The dataframe psi contains psi values for each combination ofvariables (named in the coluns A and B) computed for all columns inthe column All variables, and one column per variable named Withoutvariable_name containing the psi value when that variable is removedfrom the comparedsequences.

rkable(psi.importance$psi, digits = 4, caption = "Psi values with all variables (column 'All variables'), and without one variable at a time.")

This table can be plotted as a bar plot as follows:

``` r

extracting object

psi.df <- psi.importance$psi

to long format

psi.df.long <- tidyr::gather(psi.df, variable, psi, 3:ncol(psi.df))

creating column with names of the sequences

psi.df.long$name <- paste(psi.df.long$A, psi.df.long$B, sep=" - ")

plot

ggplot(data=psi.df.long, aes(x=variable, y=psi, fill=psi)) + geombar(stat = "identity") + coordflip() + facetwrap("name") + scalefill_viridis(direction = -1) + ggtitle("Contribution of separated variables to dissimilarity.") + labs(fill = "Psi")```

The second table, named psi.drop describes the drop in psi values,in percentage, when the given variable is removed from the analysis.Large positive numbers indicate that dissimilarity drops (increase insimilarity) when the given variable is removed, confirming that thevariable is important to explain the dissimilarity between bothsequences. Negative values indicate an increase in dissimilarity betweenthe sequences when the variable is dropped.

In summary:

• High psi-drop value: variable contributes to dissimilarity.
• Low or negative psi-drop value: variable contributes tosimilarity.

rkable(psi.importance$psi.drop, caption = "Drop in psi, as percentage of the psi value obtained when using all variables. Positive values indicate that the sequences become more similar when the given variable is removed (contribution to dissimilarity), while negative values indicate that the sequences become more dissimilar when the variable is removed (contribution to similarity).")

``` r

extracting object

psi.drop.df <- psi.importance$psi.drop

to long format

psi.drop.df.long <- tidyr::gather(psi.drop.df, variable, psi, 3:ncol(psi.drop.df))

creating column with names of the sequences

psi.drop.df.long$name <- paste(psi.drop.df.long$A, psi.drop.df.long$B, sep=" - ")

plot

ggplot(data=psi.drop.df.long, aes(x=variable, y=psi, fill=psi)) + geombar(stat = "identity") + coordflip() + geomhline(yintercept = 0, size = 0.3) + facetwrap("name") + scalefillviridis(direction = -1) + ggtitle("Drop in dissimilarity when variables are removed.") + ylab("Drop in dissimilarity (%)") + labs(fill = "Psi drop (%)")```

``` r

cleaning workspace

rm(list = ls())```

Partial matching: finding the section in a long sequence more similar to a given short sequence

In this scenario the user has one short and one long sequence, and thegoal is to find the section in the long sequence that better matches theshort one. To recreate this scenario we use the dataset sequencesMIS.The first 10 samples will serve as short sequence, and the first 40samples as long sequence. These small subsets are selected to speed-upthe execution time of this example.

``` r

loading the data

data(sequencesMIS)

removing grouping column

sequencesMIS$MIS <- NULL

subsetting to get the short sequence

MIS.short <- sequencesMIS[1:10, ]

subsetting to get the long sequence

MIS.long <- sequencesMIS[1:40, ]```

The sequences have to be prepared and transformed. For simplicity, thesequences are named short and long, and the grouping column is namedid, but the user can name them at will. Since the data representscommunity composition, a Hellinger transformation is applied.

rMIS.short.long <- prepareSequences( sequence.A = MIS.short, sequence.A.name = "short", sequence.B = MIS.long, sequence.B.name = "long", grouping.column = "id", transformation = "hellinger")

The function workflowPartialMatch shown below is going to subset thelong sequence in sizes between min.length and max.length. In theexample below this search space has the same size as MIS.short tospeed-up the execution of this example, but wider windows are possible.If left empty, the length of the segment in the long sequence to bematched will have the same number of samples as the short sequence. Inthe example below we look for segments of the same length, two samplesshorter, and two samples longer than the shorter sequence.

rMIS.psi <- workflowPartialMatch( sequences = MIS.short.long, grouping.column = "id", method = "euclidean", diagonal = TRUE, ignore.blocks = TRUE, min.length = nrow(MIS.short), max.length = nrow(MIS.short))

The function returns a dataframe with three columns: first.row (firstrow of the matched segment of the long sequence), last.row (last rowof the matched segment of the long sequence), and psi (ordered fromlower to higher). In this case, since the long sequence contains theshort sequence, the first row shows a perfectmatch.

rkable(MIS.psi[1:15, ], digits = 4, caption = "First and last row of a section of the long sequence along with the psi value obtained during the partial matching.")

Subsetting the long sequence to obtain the segment best matching withthe short sequence goes as follows.

``` r

indices of the best matching segment

best.match.indices <- MIS.psi[1, "first.row"]:MIS.psi[1, "last.row"]

subsetting by these indices

best.match <- MIS.long[best.match.indices, ]```

``` r

cleaning workspace

rm(list = ls())```

Sequence slotting: combining samples of two sequences into a single composite sequence

Under this scenario, the objective is to combine two sequences into asingle composite sequence. The basic assumption followed by thealgorithm building the composite sequence is most similar samplesshould go together, but respecting the original ordering of thesequences. Therefore, the output will contain the samples in bothsequences ordered in a way that minimizes the multivariate distancebetween consecutive samples. This scenario assumes that at least one ofthe sequences do not have a time/age/depth column, or that the values insuch a column are uncertain. In any case, time/age/depth is notconsidered as a factor in the generation of the composite sequence.

The example below uses the pollenGP dataset, which contains 200samples, with 40 pollen types each. To create a smalle case study, thecode below separates the first 20 samples of the sequence into twodifferent sequences with 10 randomly selected samples each. Even thoughthis scenario assumes that these sequences do not have depth or age,these columns will be kept so the result can be assessed. That is whythese columns are added to the exclude.columns argument. Also, notethat the argument transformation is set to “none”, so the output isnot transformed, and the outcome can be easily interpreted. This willgive more weight to the most abundant taxa, which will in fact guide theslotting.

``` r

loading the data

data(pollenGP)

getting first 20 samples

pollenGP <- pollenGP[1:20, ]

sampling indices

set.seed(10) #to get same result every timesampling.indices <- sort(sample(1:20, 10))

subsetting the sequence

A <- pollenGP[sampling.indices, ]B <- pollenGP[-sampling.indices, ]

preparing the sequences

AB <- prepareSequences( sequence.A = A, sequence.A.name = "A", sequence.B = B, sequence.B.name = "B", grouping.column = "id", exclude.columns = c("depth", "age"), transformation = "none")```

Once the sequences are prepared, the function workflowSlotting willallow to combine (slot) them. The function computes a distance matrixbetween the samples in both sequences according to the methodargument, computes the least-cost matrix, and generates the least-costpath. Note that it only uses an orthogonal method considering blocks,since this is the only option really suitable for this task.

rAB.combined <- workflowSlotting( sequences = AB, grouping.column = "id", time.column = "age", exclude.columns = "depth", method = "euclidean", plot = TRUE)

The function reads the least-cost path in order to find the combinationof samples of both sequences that minimizes dissimilarity, constrainedby the order of the samples on each sequence. The output dataframe has acolumn named original.index, which has the index of each sample in theoriginaldatasets.

rkable(AB.combined[1:15,1:10], digits = 4, caption = "Combination of sequences A and B.")

Note that several samples show inverted ages with respect to theprevious samples. This is to be expected, since the slotting algorithmonly takes into account distance/dissimilarity between adjacent samplesto generate theordering.

``` r

cleaning workspace

rm(list = ls())```

Transferring an attribute from one sequence to another

This scenario assumes that the user has two METS, one of them with agiven attribute (age/time) that needs to be transferred to the othersequence by using similarity/dissimilarity (constrained by sample order)as a transfer criterion. This case is relatively common inpalaeoecology, when a given dataset is dated, and another taken at aclose location is not.

The code below prepares the data for the example. The sequencepollenGP is the reference sequence, and contains the column age.The sequence pollenX is the target sequence, without an agecolumn. We generate it by taking 40 random samples between the samples50 and 100 of pollenGP. The sequences are prepared withprepareSequences, as usual, with the identificators “GP” and “X”

``` r

loading sample dataset

data(pollenGP)

subset pollenGP to make a shorter dataset

pollenGP <- pollenGP[1:50, ]

generating a subset of pollenGP

set.seed(10)pollenX <- pollenGP[sort(sample(1:50, 40)), ]

we separate the age column

pollenX.age <- pollenX$age

and remove the age values from pollenX

pollenX$age <- NULLpollenX$depth <- NULL

removing some samples from pollenGP

so pollenX is not a perfect subset of pollenGP

pollenGP <- pollenGP[-sample(1:50, 10), ]

prepare sequences

GP.X <- prepareSequences( sequence.A = pollenGP, sequence.A.name = "GP", sequence.B = pollenX, sequence.B.name = "X", grouping.column = "id", time.column = "age", exclude.columns = "depth", transformation = "none")```

The transfer of “age” values from GP to X can be done in two ways,both constrained by sample order:

• Direct: each sample in X gets the age of its most similarsample in GP.
• Interpolated: each sample in X gets an age interpolated fromthe ages of the two most similar samples in GP. The interpolationis weighted by the similarity between the samples.

Direct transfer

A direct transfer of an attribute from the samples of one sequence tothe samples of another requires to compute a distance matrix betweensamples, the least-cost matrix and its least-cost path (both with theoption diagonal activated), and to parse the least-cost path file toassign attribute values. This is done by the functionworkflowTransfer with the option .

``` r

parameters

X.new <- workflowTransfer( sequences = GP.X, grouping.column = "id", time.column = "age", method = "euclidean", transfer.what = "age", transfer.from = "GP", transfer.to = "X", mode = "direct" )

kable(X.new[1:15, ], digits = 4)```

The algorithm finds the most similar samples, and transfers attributevalues directly between them. This can result in duplicated attributevalues, as highlighted in the table above. The Pearson correlationbetween the original ages (stored in pollenX.age) and the assignedones is 0.9996, so it can be concluded that in spite of its simplicity,this algorithm yields accurate results.

Interpolated transfer

If we consider:

• Two samples and ofthe sequence .
• Each one with a value the attribute,and .
• One sample of the sequence .
• With an unknown attribute.
• The multivariate distance between the samplesand .
• The multivariate distance between the samplesand .
• The weight , computed as
• The weight , computed as

The unknwon value iscomputed as:

The code below exemplifies the operation, using the samples 1 and 4 ofthe dataset pollenGP as and , and thesample 3 as .

``` r

loading data

data(pollenGP)

samples in A

Ai <- pollenGP[1, 3:ncol(pollenGP)]Aj <- pollenGP[4, 3:ncol(pollenGP)]

ages of the samples in A

Ati <- pollenGP[1, "age"]Atj <- pollenGP[4, "age"]

sample in B

Bk <- pollenGP[2, 3:ncol(pollenGP)]

computing distances between Bk, Ai, and Aj

DBkAi <- distance(Bk, Ai)DBkAj <- distance(Bk, Aj)

normalizing the distances to 1

wi <- DBkAi / (DBkAi + DBkAj)wj <- DBkAj / (DBkAi + DBkAj)

computing Btk

Btk <- wi * Ati + wj * Atj```

The table below shows the observed versus the predicted values for.

rtemp.df <- data.frame(Observed = pollenGP[3, "age"], Predicted = Btk)kable(t(temp.df), digits = 4, caption = "Observed versus predicted value in the interpolation of an age based on similarity between samples.")

Below we create some example data, where a subset of pollenGP will bethe donor of age values, and another subset of it, named pollenX willbe the receiver of the age values.

``` r

loading sample dataset

data(pollenGP)

subset pollenGP to make a shorter dataset

pollenGP <- pollenGP[1:50, ]

generating a subset of pollenGP

set.seed(10)pollenX <- pollenGP[sort(sample(1:50, 40)), ]

we separate the age column

pollenX.age <- pollenX$age

and remove the age values from pollenX

pollenX$age <- NULLpollenX$depth <- NULL

removing some samples from pollenGP

so pollenX is not a perfect subset of pollenGP

pollenGP <- pollenGP[-sample(1:50, 10), ]

prepare sequences

GP.X <- prepareSequences( sequence.A = pollenGP, sequence.A.name = "GP", sequence.B = pollenX, sequence.B.name = "X", grouping.column = "id", time.column = "age", exclude.columns = "depth", transformation = "none")```

To transfer attributes from GP to X we use the workflowTransferfunction with the option mode = “interpolate”.

``` r

parameters

X.new <- workflowTransfer( sequences = GP.X, grouping.column = "id", time.column = "age", method = "euclidean", transfer.what = "age", transfer.from = "GP", transfer.to = "X", mode = "interpolated" )

kable(X.new[1:15, ], digits = 4, caption = "Result of the transference of an age attribute from one sequence to another. NA values are expected when predicted ages for a given sample yield a higher number than the age of the previous sample.") %>% row_spec(c(8, 13), bold = T)```

When interpolated values of the age column (transferred attribute viainterpolation) show the value NA, it means that the interpolationyielded an age lower than the previous one. This happens when the same and are used to evaluatetwo or more different samples, and the second is more similar to than the first one.These NA values can be removed with na.omit(), or interpolated withthe functionsimputeTS::na.interpolationorzoo::na.approx.

Without taking into account these NA values, the Pearson correlationof the interpolated ages with the real ones is 0.9985.

IMPORTANT: the interpretation of the interpolated ages requires acareful consideration. Please, don’t do it blindly, because thisalgorithm has its limitations. For example, significant hiatuses in thedata can introduce wild variations in interpolated ages.

``` r

cleaning workspace

rm(list = ls())```

wget https://archive.org/download/github.com-BlasBenito-distantia_-_2021-03-03_23-13-27/BlasBenito-distantia_-_2021-03-03_23-13-27.bundle

 git clone BlasBenito-distantia_-_2021-03-03_23-13-27.bundle 

plus-circle Add Review

comment

Reviews