Vector-valued functional data: gait cycles and hurricane tracks

A tf vector represents a sample of functions \(f: \mathbb{R} \to \mathbb{R}\) (see the tf vectors vignette for an introduction to these classes and their operations). Many real measurement processes, though, produce several coupled signals that share a common argument: the (hip, knee) joint-angle pair sampled across a gait cycle, the (longitude, latitude) position of a hurricane sampled in time, the (x, y, z) body coordinates of a moving animal. The natural object is a vector-valued function \(f: \mathbb{R} \to \mathbb{R}^d\), and tf represents a sample of those with the tf_mv family: tfd_mv for raw evaluations and tfb_mv for basis representations. Internally a tf_mv is a bundle of \(d\) ordinary tf vectors, so the multivariate methods shown below reuse the existing univariate numerical kernels by mapping over components.

This article uses two real datasets to put the API through its paces:

Code 
library(tf)
library(tidyfun)
library(dplyr)

1. A quick tour of vector-valued functional data

Before the case studies, this part introduces the objects: what a vector-valued function is, the two classes that represent samples of them, how to construct them, and the handful of geometric operations the case studies lean on.

What is a vector-valued function?

A tf vector stores a sample of scalar functions \(f: [a, b] \to \mathbb{R}\). A vector-valued function bundles \(d\) of them on a shared argument,

\[ f(t) = \big(f_1(t), \dots, f_d(t)\big), \qquad t \in [a, b], \]

so that as \(t\) runs over the domain the value \(f(t)\) traces a path through \(\mathbb{R}^d\). tf represents a sample of such paths with the tf_mv family. Internally a tf_mv is a bundle of d ordinary tf vectors sharing one argument, so every univariate verb – evaluation, arithmetic, derivatives, integration, smoothing, basis fitting – extends component-wise for free.

Two kinds of variation run through everything below: amplitude (how large the values get) and phase (where in \(t\) the features happen). Much of the work in functional data analysis (Ramsay and Silverman 2005) is telling the two apart (Marron et al. 2015), and Section 2 is largely about exactly that.

The tf_mv classes: tfd_mv and tfb_mv

Like their univariate cousins tfd/tfb, vector-valued functions come in two flavours:

Components are named, and a few accessors reach into the bundle: tf_ncomp() (number of components), tf_components() (the list of univariate tf vectors), tf_component(x, k) (one component, by name or index), and names().

Constructing tf_mv objects

tfd_mv() accepts three interchangeable input layouts. We use the built-in gait data (Olshen et al. 1989; Ramsay and Silverman 2005) throughout the gait case study – hip and knee angles, in degrees, for 39 boys, each sampled at 20 points across one gait cycle. The most direct construction is a named list of univariate tf vectors:

Code 
data(gait)
g <- tfd_mv(list(hip = gait$hip_angle, knee = gait$knee_angle))
g
#> tfd_mv<d=2>[39] (hip, knee): [0.025, 0.975] -> [-12, 64] x [0, 82]
#> components based on 20 evaluations each, interpolation by tf_approx_linear
#> [1]: ▆▆▅▅▄▃▃▃▂▂▂▂▃▄▅▆▇▇▆▆ | ▂▂▂▂▂▂▂▂▂▃▃▄▆▇█▇▆▅▃▂
#> [2]: ▇▇▇▅▅▄▃▃▂▁▁▁▃▄▅▇▇▇▇▇ | ▂▃▃▃▂▂▂▁▁▁▂▄▆▇█▇▆▄▃▂
#> [3]: ▇▇▆▅▅▅▄▃▂▁▁▁▂▄▅▇███▇ | ▂▃▄▄▃▃▂▂▁▂▂▄▆███▇▅▃▂
#> [4]: ▆▆▅▄▃▃▂▂▁▁▁▂▂▃▄▅▅▆▅▅ | ▁▂▂▂▂▁▁▁▁▁▃▅▆▇▇▇▅▂▁▁
#> [5]: ▄▃▂▂▂▂▁▁▁▁▁▁▂▃▅▆▅▅▅▅ | ▁▁▁▁▁▁▁▁▁▂▃▅▆▇█▇▅▃▁▁
#> [6]: █▇▇▆▅▅▄▃▃▂▂▂▂▄▅▆▇███ | ▂▂▃▃▃▂▂▂▁▂▂▃▅▆▇▇▇▅▃▁
#> 
#>     [....]   (33 not shown)

The print-out shows d = 2 components on a shared grid, with a sparkline pair per subject. It behaves like a vector of 39 bivariate curves; component accessors and bracket extraction expose the two levels separately – curves on rows, argument values on columns, components on the third array dimension:

Code 
tf_ncomp(g)
#> [1] 2
names(tf_components(g))
#> [1] "hip"  "knee"
dim(g[1:3, c(0.1, 0.5, 0.9)])
#> [1] 3 3 2
class(g[, , component = "hip"])
#> [1] "tfd_reg"    "tfd"        "tf"         "vctrs_vctr" "list"

The same object can be built from a 3-d array [curve, arg, component] – the layout as.matrix() returns –

Code 
arr <- as.matrix(g)            # [curve, arg, component]
dim(arr)
#> [1] 39 20  2
g_arr <- tfd_mv(arr, arg = tf_arg(g))
identical(as.matrix(g_arr), arr)
#> [1] TRUE

or from a wide data frame with one column per component (the long (id, arg, component, value) schema is what as.data.frame(unnest = TRUE) emits; long = FALSE gives the wide schema the constructor consumes):

Code 
df <- as.data.frame(g, unnest = TRUE, long = FALSE)   # id, arg, hip, knee
head(df, 2)
#>     id   arg hip knee
#> 1 boy1 0.025  37   10
#> 2 boy1 0.075  36   15
g_df <- tfd_mv(df, id = "id", arg = "arg", value = c("hip", "knee"))

Each route also takes an optional domain and a per-component evaluator (interpolation rule, default tf_approx_linear). A tfd_mv becomes a basis object with tfb_mv(), choosing a "spline" or "fpc" basis; basis arguments apply globally or per component when passed as a named list:

Code 
g_spline <- tfb_mv(g, basis = "spline", k = list(hip = 5, knee = 12),
                   verbose = FALSE)
g_spline
#> tfb_mv<d=2>[39] (hip, knee): [0.025, 0.975] -> [-12.68226, 64.90069] x [-0.3005807, 81.6327]
#>   hip: in basis representation: s(arg, bs = "cr", k = 5, sp = -1)
#>   knee: in basis representation: s(arg, bs = "cr", k = 12, sp = -1)
#> [1]: ▆▅▅▅▄▄▃▂▂▂▂▂▃▄▅▆▆▆▆▆ | ▂▂▂▂▂▂▂▂▂▃▃▄▆▇█▇▆▅▃▂
#> [2]: ▇▆▆▆▅▄▃▂▁▁▁▂▃▄▅▆▆▇▇▇ | ▂▃▃▃▂▂▂▁▁▁▂▄▆▇█▇▆▄▃▂
#> [3]: ▆▆▆▆▅▅▃▂▁▁▁▂▃▄▅▆▇▇▇█ | ▂▃▄▄▃▃▂▂▁▂▂▄▆███▇▅▃▂
#> [4]: ▆▅▅▄▃▃▂▂▁▁▁▂▃▃▄▅▅▅▅▅ | ▁▂▂▂▂▁▁▁▁▁▃▄▆▇▇▇▅▂▁▁
#> [5]: ▃▃▃▂▂▂▁▁▁▁▁▂▃▄▄▅▅▅▅▅ | ▁▁▁▁▁▁▁▁▁▂▃▅▆██▇▅▃▁▁
#> [6]: █▇▆▆▅▄▄▃▂▂▂▃▃▄▅▆▇▇██ | ▂▂▃▃▃▂▂▂▁▂▂▃▅▆▇▇▆▅▃▁
#> 
#>     [....]   (33 not shown)

Geometry on the bundle: speed, arc length, reparametrization

Because a tf_mv traces a path, the usual curve-geometry quantities are available, and they return univariate summaries of the multivariate object:

Code 
head(tf_arclength(g))
#>     boy1     boy2     boy3     boy4     boy5     boy6 
#> 166.6539 200.2389 223.4584 193.6448 177.4843 191.0575
plot(tf_speed(g), alpha = 0.3, main = expression("pointwise speed " * "||" * f*"'"*(t) * "||"))

Aligning vector-valued curves: a ladder

When curves differ mainly in timing, we usually want to factor that out – but “factor out” can mean several different things. There is a ladder of increasingly aggressive operations, distinguished by what each treats as a nuisance (Marron et al. 2015):

rung operation removes keeps tf entry point
1 arc-length reparametrization the clock (parametrization only) shape and size tf_reparam_arclength()
2 warp from a 1-d reference signal phase, onto a chosen reference amplitude, size, orientation tf_register(., method = "cc")
3 joint multivariate warp phase, from all components at once amplitude, size, orientation tf_register(., method = "srvf_mv")
4 full elastic shape registration phase + rotation + scale shape only tf_register_shape()

Rung 1 re-labels time within each curve without picking a template; rungs 2 and 3 warp every curve onto a common template (from one reference signal, or from all components jointly) but leave the values untouched; rung 4 additionally rotates and rescales each curve, landing in a shape space where only the geometric form remains. Section 2 walks all four rungs on the gait data – and shows where the top rung is, and is not, the right tool.

2. How variable is a “normal” gait cycle?

We built the two-component object g above. The two natural views of these 39 curves are the time-series (type = "facet", one panel per component) and the trajectory in phase space (type = "trajectory", the default for two-component objects).

Code 
plot(g, type = "facet", alpha = 0.4)

Code 
plot(g, alpha = 0.4)   # type = "trajectory" by default for d = 2

The facet view shows that both joints flex twice per cycle; the trajectory view collapses time and reveals the characteristic “butterfly” loop traced out as the leg cycles through stance and swing. Both views suggest the between-subject spread is far from uniform along the cycle.

Pointwise mean and standard deviation

mean() and sd() are vctrs group generics and dispatch component-wise: they each return a length-1 tfd_mv. Plot the raw curves with plot.tf, then overlay the mean +/- 2 sd envelope with lines.tf – all arithmetic on the components is component-wise on tf_mv, so mu + 2 * s is itself a tfd_mv and tf_component(...) gives the per-axis envelope:

Code 
mu <- mean(g)
s  <- sd(g)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
for (k in seq_len(tf_ncomp(g))) {
  nm <- names(tf_components(g))[k]
  plot(tf_component(g, k), alpha = 0.25,
       ylab = nm, main = nm)
  lines(tf_component(mu, k),              col = "red", lwd = 2)
  lines(tf_component(mu + 2 * s, k),      col = "red", lwd = 1.2, lty = 2)
  lines(tf_component(mu - 2 * s, k),      col = "red", lwd = 1.2, lty = 2)
}

Code 
par(op)

The pointwise sd peaks roughly where the angle itself is changing fastest – around the stance/swing transitions – so most of the between-subject variability lives in timing of those transitions, not in the angles attained at rest.

Most- and least-varied subjects

tf_arclength() measures the path length traced out in (hip, knee)-space. A short orbit = a subject whose (hip, knee) excursions are small or who moves the two joints in lock-step; a long orbit = a subject with large-amplitude, less synchronised excursions.

Code 
arc <- tf_arclength(g)
extreme <- c(which.min(arc), which.max(arc))
extreme
#>  boy1 boy39 
#>     1    39

plot(g, alpha = 0.15)
lines(g[extreme[1]], col = "steelblue", lwd = 2.2)
lines(g[extreme[2]], col = "firebrick", lwd = 2.2)
legend("bottomright", bty = "n",
       legend = c(sprintf("min arc length (%.0f deg)", arc[extreme[1]]),
                  sprintf("max arc length (%.0f deg)", arc[extreme[2]])),
       col = c("steelblue", "firebrick"), lwd = 2)

A ladder of registrations

The pointwise sd above conflates two kinds of between-subject variability: amplitude (how far the knee swings) and phase (at what fraction of the cycle heel-strike happens). Removing phase – “registration” – can mean progressively more, depending on what we are willing to treat as a nuisance. Below are four rungs, applied in turn to the gait sample, each followed by a facet plot of the result; we quantify them together at the end.

Rungs 3 and 4 (and the rung-comparison at the end) use the elastic registration routines from the suggested fdasrvf package – those code chunks are only evaluated if fdasrvf is installed.

Rung 1 – arc-length reparametrization

The mildest operation re-labels time within each curve and uses no template at all. tf_reparam_arclength() traverses every curve at approximately constant speed in its value space: the (hip, knee) path is unchanged – the same set of points – but equal time intervals now cover equal arc length.

Code 
g_unit <- tf_reparam_arclength(g)
plot(g_unit, type = "facet", alpha = 0.3)

tf_speed() makes the effect explicit. The raw speeds swing wildly (fast through swing, near-still in stance); after reparametrization they are nearly flat – each curve traversed at constant speed:

Code 
sp_raw  <- tf_speed(g)
sp_unit <- tf_speed(g_unit)
yl <- c(0, max(c(unlist(tf_evaluations(sp_raw)),
                 unlist(tf_evaluations(sp_unit))), na.rm = TRUE))

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(sp_raw,  alpha = 0.4, ylim = yl, main = "raw parameterization",
     ylab = "speed (deg / cycle)")
lines(mean(sp_raw),  col = "firebrick", lwd = 2)
plot(sp_unit, alpha = 0.4, ylim = yl, main = "arc-length parameterization",
     ylab = "speed (deg / cycle)")
lines(mean(sp_unit), col = "firebrick", lwd = 2)

Code 
par(op)

No curve is moved toward any other: this removes the parametrization only, leaving shape and size untouched and aligning nothing to a common reference.

Rung 2 – alignment to a 1-d reference signal

The next rung aligns curves to each other by estimating one shared time-warp per curve and applying it to both components. The warp is driven by a single univariate reference signal – a component, or any derived quantity. The choice matters: the knee carries the cycle’s sharpest event (heel-strike), so it locks the alignment best.

Code 
r_knee <- tf_register(g, method = "cc", ref_component = "knee")
plot(tf_aligned(r_knee), type = "facet", alpha = 0.3)

The knee component visibly tightens; the hip, phase-coupled to it, tightens only incidentally. Choosing ref_component = "hip" or = tf_speed would emphasise different events – a 1-d reference forces that modelling choice. The warps themselves are univariate tfds mapping raw to aligned phase, fanning out around the swing phase where heel-strike timing varies between subjects:

Code 
plot(tf_invert(tf_inv_warps(r_knee)), alpha = 0.5,
     main = "reference (knee) warps", ylab = "aligned phase")
abline(0, 1, lty = 3)

Rung 3 – joint multivariate reparametrization (srvf_mv)

method = "srvf_mv" removes the choice of reference: it estimates the single time-warp that best aligns the joint (hip, knee) trajectories in the elastic (square-root-velocity) sense, using all components at once (Srivastava, Wu, et al. 2011; Tucker, Wu, and Srivastava 2013). It is still only a re-timing – values are never rotated or rescaled – and the template is the multivariate Karcher mean in tf_template(reg_mv).

Code 
reg_mv <- tf_register(g, method = "srvf_mv", max_iter = 2)
plot(tf_aligned(reg_mv), type = "facet", alpha = 0.3)

Code 
plot(tf_invert(tf_inv_warps(reg_mv)), alpha = 0.5,
     main = "srvf_mv warps", ylab = "aligned phase")
abline(0, 1, lty = 3)

Because the warp must compromise between both components rather than chase one, it need not shrink any single channel as hard as a reference-targeted warp – but it needs no reference choice and respects both axes symmetrically.

Rung 4 – full elastic shape registration

tf_register_shape() goes furthest: on top of the time-warp it also fits, per curve, a rotation and a scale, and reports the aligned curves in a centered, normalised shape space (Srivastava, Klassen, et al. 2011; Srivastava and Klassen 2016). Translation, rotation and size are all quotiented out.

Code 
reg_shape <- tf_register_shape(g, max_iter = 2)
plot(tf_aligned(reg_shape), type = "facet", alpha = 0.3)   # note the y-axis: shape space

The (hip, knee) trajectory view makes the collapse vivid – the original butterfly loops on the left, the shape-registered curves on the right reduced to essentially one normalised shape (note the very different axis ranges):

Code 
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(g,                     type = "trajectory", alpha = 0.3, main = "original")
plot(tf_aligned(reg_shape), type = "trajectory", alpha = 0.3,
     main = "shape-registered (shape space)")

Code 
par(op)

The estimated rotations sit near the identity and the scales near 1 (little genuine rotation or scaling is present), yet the normalisation alone flattens the real between-subject differences:

Code 
round(tf_rotations(reg_shape)[, , 1], 3)   # rotation for subject 1
#>         hip  knee
#> hip   1.000 0.023
#> knee -0.023 1.000
summary(as.numeric(tf_scales(reg_shape)))  # scales, relative to the template
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   178.0   199.9   213.2   214.7   229.5   256.8

Gait’s two axes are not interchangeable: hip angle and knee angle are distinct physical quantities in fixed units. Quotienting out rotation mixes them into meaningless combinations, and quotienting out scale discards amplitude – which here is the signal of interest. Shape registration is the wrong tool for a bundle of fixed-unit channels; we return to where it is right below.

Quantifying the alignment

Two views across the rungs. First, a visual overview: the two components (rows hip, knee) and the estimated time-warp (row warp) for each registration, side by side – as-observed, arc-length reparametrized, registered to a 1-d reference (knee), registered jointly (srvf_mv), and fully shape-registered. Raw data has no warp; arc-length has no template warp but does reparametrize, so its warp panel shows the normalised cumulative arc length.

Code 
# warp implied by arc-length reparametrization: normalised cumulative arc length
# (raw cycle phase -> arc-length phase, same direction as the registration warps)
arclen_warp <- function(mv) {
  sp  <- tf_speed(mv)
  arg <- as.numeric(tf_arg(sp))
  dom <- tf_domain(sp)
  d   <- diff(arg)
  w <- sapply(tf_evaluations(sp), function(s) {
    cum <- c(0, cumsum((head(s, -1) + tail(s, -1)) / 2 * d))
    dom[1] + diff(dom) * cum / cum[length(cum)]
  })
  tfd(t(w), arg = arg, domain = dom)
}
Code 
cols <- list(
  list(nm = "as-observed", data = g,                     warp = NULL),
  list(nm = "arc-length",  data = g_unit,                warp = arclen_warp(g)),
  list(nm = "ref = knee",  data = tf_aligned(r_knee),    warp = tf_invert(tf_inv_warps(r_knee))),
  list(nm = "srvf_mv",     data = tf_aligned(reg_mv),    warp = tf_invert(tf_inv_warps(reg_mv))),
  list(nm = "shape",       data = tf_aligned(reg_shape), warp = tf_invert(tf_inv_warps(reg_shape)))
)

op <- par(mfrow = c(3, 5), mar = c(2.6, 3.2, 2, 0.8), mgp = c(1.9, 0.6, 0))
for (j in seq_along(cols))
  plot(tf_component(cols[[j]]$data, "hip"),  alpha = 0.3,
       main = cols[[j]]$nm, ylab = if (j == 1) "hip" else "")
for (j in seq_along(cols))
  plot(tf_component(cols[[j]]$data, "knee"), alpha = 0.3,
       main = "", ylab = if (j == 1) "knee" else "")
for (j in seq_along(cols)) {
  w <- cols[[j]]$warp
  if (is.null(w)) {
    plot.new(); text(0.5, 0.5, "(identity)", col = "grey55")
  } else {
    plot(w, alpha = 0.4, main = "", ylab = if (j == 1) "warp" else "")
    abline(0, 1, lty = 3)
  }
}

Code 
par(op)

Across the warp row, the registration warps fan out around the swing phase where heel-strike timing varies between subjects; the shape warp is comparable, since it shares the same time-warp step. Down the shape column the components have been rescaled into normalised shape space – note the compressed vertical axis – which is why those curves look so much tighter than the others.

Second, the numbers: the maximum pointwise sd of each component, before and after. Rungs 1–3 keep the data in its original degree units and are directly comparable; rung 4 (marked *) lives in normalised shape space, so its row is on a different scale and only illustrates the collapse.

Code 
peak_sd <- function(f, k) {
  arg <- tf_arg(f)
  max(unlist(tf_evaluate(tf_component(sd(f), k), arg)))
}
Code 
rungs <- list(
  "raw"                = g,
  "1 arc-length"       = g_unit,
  "2 reference (knee)" = tf_aligned(r_knee),
  "3 srvf_mv"          = tf_aligned(reg_mv),
  "4 shape (*)"        = tf_aligned(reg_shape)
)
data.frame(
  registration = names(rungs),
  sd_hip       = sapply(rungs, peak_sd, "hip"),
  sd_knee      = sapply(rungs, peak_sd, "knee"),
  row.names    = NULL
)
#>         registration     sd_hip    sd_knee
#> 1                raw 8.28498797 9.54533186
#> 2       1 arc-length 8.04009724 7.39283787
#> 3 2 reference (knee) 7.91270379 6.33047672
#> 4          3 srvf_mv 8.20232891 8.22321203
#> 5        4 shape (*) 0.01935199 0.02937269

In degree units the phase rungs chip away at the spread only modestly – gait’s phase variation is real but subtle. The knee-reference warp shrinks the knee component most (it targets that channel); srvf_mv spreads a gentler correction across both; and arc-length reparametrization, which aligns nothing to a common template, barely moves the pointwise sd at all. Rung 4’s near-zero entries are not a better alignment – they are the signal being quotiented away.

Shape registration and its quotient spaces

Shape registration earns its keep when the components really are interchangeable spatial coordinates and position, orientation and size are all nuisances: handwriting, gesture or movement paths, object outlines. Here is a synthetic example – a single base curve \((t,\, t^2)\), copied three times with each copy rotated, rescaled and shifted in the plane:

Code 
t      <- seq(0, 1, length.out = 51)
base   <- rbind(t, t^2)
scales <- c(1, 0.7, 1.3)
angles <- c(0, 0.4, -0.25)
offset <- rbind(c(0.2, -0.1), c(0.4, -0.2), c(0.6, -0.3))
beta   <- array(NA_real_, dim = c(3, length(t), 2),
                dimnames = list(c("a", "b", "c"), NULL, c("x", "y")))
for (i in 1:3) {
  rot   <- matrix(c(cos(angles[i]),  sin(angles[i]),
                    -sin(angles[i]), cos(angles[i])), nrow = 2)
  curve <- scales[i] * (rot %*% base) + matrix(offset[i, ], 2, length(t))
  beta[i, , 1] <- curve[1, ]
  beta[i, , 2] <- curve[2, ]
}
shapes <- tfd_mv(beta, arg = t)

What counts as “the same shape” depends on which transformations we quotient out, and tf_register_shape() lets us choose via the rotation and scale flags (translation and the time-warp are always removed). The three quotients differ in what survives (again, these chunks need the suggested fdasrvf package):

Code 
reg_full  <- tf_register_shape(shapes, max_iter = 2, rotation = TRUE,  scale = TRUE)
reg_rot   <- tf_register_shape(shapes, max_iter = 2, rotation = TRUE,  scale = FALSE)
reg_scale <- tf_register_shape(shapes, max_iter = 2, rotation = FALSE, scale = TRUE)

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(shapes,                asp = 1, col = 1:3, lwd = 2, main = "input: 3 placements")
plot(tf_aligned(reg_rot),   asp = 1, col = 1:3, lwd = 2, main = "rotation-only (keeps size)")
plot(tf_aligned(reg_scale), asp = 1, col = 1:3, lwd = 2, main = "scale-only (keeps orientation)")
plot(tf_aligned(reg_full),  asp = 1, col = 1:3, lwd = 2, main = "rotation + scale (full)")

Code 
par(op)

tf_scales() reports the per-curve size factors that were removed (relative to the template; 1 whenever scaling is off):

Code 
data.frame(curve      = c("a", "b", "c"),
           injected   = scales,
           full       = round(as.numeric(tf_scales(reg_full)),  3),
           rot_only   = round(as.numeric(tf_scales(reg_rot)),   3),
           scale_only = round(as.numeric(tf_scales(reg_scale)), 3))
#>   curve injected  full rot_only scale_only
#> 1     a      1.0 1.463        1      1.466
#> 2     b      0.7 1.024        1      1.026
#> 3     c      1.3 1.902        1      1.905

Modes of variation via FPC

Back on the gait sample, with the most visible phase variation reduced (we carry the knee-aligned curves forward), fit a per-component FPC basis (Ramsay and Silverman 2005). The leading PCs now mostly describe within-component amplitude modes:

Code 
fpc_var <- function(comp) attr(comp, "score_variance")
Code 
g_aligned <- tf_aligned(r_knee)
g_b <- tfb_mv(g_aligned, basis = "fpc", verbose = FALSE)

v_hip  <- fpc_var(tf_component(g_b, "hip"));  v_hip  <- v_hip  / sum(v_hip)
v_knee <- fpc_var(tf_component(g_b, "knee")); v_knee <- v_knee / sum(v_knee)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
barplot(head(v_hip,  6), names.arg = 1:6, main = "hip: variance share",
        ylab = "proportion", col = "grey70")
barplot(head(v_knee, 6), names.arg = 1:6, main = "knee: variance share",
        ylab = "proportion", col = "grey70")

Code 
par(op)

# residual: how well does the low-rank FPC basis approximate the curves?
g_round <- vctrs::vec_cast(g_b, g_aligned)
resid   <- g_aligned - g_round
rmse_per_subject <- sqrt(unlist(lapply(tf_evaluations(resid), function(m) {
  mean(as.matrix(m[, -1L, drop = FALSE])^2)
})))
summary(rmse_per_subject)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1685  0.2931  0.3377  0.3485  0.4048  0.5658

The first 2-3 FPCs explain most of the remaining within-component variance, and the RMSE between the registered and FPC-reconstructed (hip, knee) curves is well below a degree for most subjects.

Joint modes of variation via MFPCA

The per-component FPC above gives each component its own scores: a subject’s hip-mode-1 score and knee-mode-1 score are separate numbers, so the basis cannot, on its own, express that the two joints co-vary. Multivariate FPCA (Happ and Greven 2018), via tfb_mfpc(), instead produces a single score per subject per mode, paired with a vector-valued eigenfunction \(\Psi_m = (\Psi_m^{\text{hip}}, \Psi_m^{\text{knee}})\). One shared score \(s_{im}\) scales both halves at once, so each mode encodes coupled hip-knee co-variation in a single coordinate system: \[ f_i(t) \approx \mu(t) + \sum_m s_{im}\,\Psi_m(t). \]

The default weights = "inverse_variance" normalises each component to equal total variance before combining, so the (larger-range) knee does not simply dominate the hip; "snr", "equal", or a numeric vector are also available.

Code 
g_m <- tfb_mfpc(g_aligned, pve = 0.95)
g_m
#> tfb_mv<d=2>[39] (hip, knee): [0.025, 0.975] -> [-10.93848, 65.14214] x [0.405835, 80.19255]
#> components in basis representation: 7 MFPCs
#> [1]: ▆▆▅▅▄▄▃▃▂▂▂▂▃▄▅▆▇▇▇▆ | ▁▂▂▂▂▂▂▂▂▂▃▄▆▇▇▇▆▄▃▂
#> [2]: ▇▇▆▅▄▄▃▂▁▁▁▂▃▄▆▇▇▇▇▇ | ▂▃▃▃▂▂▁▁▁▂▃▄▆▇█▇▆▅▃▂
#> [3]: ▇▇▆▅▄▃▂▂▁▁▁▂▃▄▆▇███▇ | ▂▃▄▃▃▂▂▁▁▂▃▅▇███▇▅▃▂
#> [4]: ▆▆▅▄▄▃▃▂▁▁▁▁▂▃▄▅▅▆▅▅ | ▁▂▂▂▂▁▁▁▁▁▂▄▆▇▇▇▆▄▁▁
#> [5]: ▄▃▂▂▂▁▁▁▁▁▁▁▂▃▄▅▆▆▅▅ | ▁▁▁▁▁▁▁▁▁▂▃▄▆▇█▇▆▄▁▁
#> [6]: █▇▆▅▄▄▃▃▂▂▂▂▃▄▆▇████ | ▂▃▃▂▂▂▂▁▂▂▃▄▆▇▇▇▆▄▂▂
#> 
#>     [....]   (33 not shown)

scores <- tf_mfpc_scores(g_m)          # 39 x M, shared across both components
nu     <- attr(g_m, "mfpc")$evalues    # joint (multivariate) eigenvalues
ve     <- nu / sum(nu)                  # variance share per shared mode
dim(scores)
#> [1] 39  7
round(head(ve, 4), 3)
#> [1] 0.465 0.224 0.138 0.084

The multivariate eigenfunctions come back as a tfd_mv – one bivariate “curve” per mode – so each MFPC can be inspected component by component. The leading mode’s hip and knee parts move together, which is exactly the coupling the independent per-component analysis could not represent:

Code 
psi    <- tf_mfpc_efunctions(g_m)
k_show <- min(3L, length(psi))
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(tf_component(psi, "hip")[seq_len(k_show)], col = seq_len(k_show), lwd = 2,
     main = "MFPC eigenfunctions: hip", ylab = "loading")
plot(tf_component(psi, "knee")[seq_len(k_show)], col = seq_len(k_show), lwd = 2,
     main = "MFPC eigenfunctions: knee", ylab = "loading")
legend("topright", bty = "n", lwd = 2, col = seq_len(k_show),
       legend = paste("MFPC", seq_len(k_show)), cex = 0.85)

Code 
par(op)

The shared scores live in one coordinate system, so MFPCA summarises each subject with a single set of M numbers rather than one set per component. That is a genuine trade-off: the independent per-component basis is free to optimise each component separately, so it reconstructs more accurately for a given per-component truncation, while MFPCA spends fewer, coupled scores and accepts a little more error in exchange for the shared coordinate system:

Code 
rmse_mv <- function(approx) {
  resid <- g_aligned - vctrs::vec_cast(approx, g_aligned)
  sqrt(mean(unlist(lapply(tf_evaluations(resid), function(m) {
    as.matrix(m[, -1L, drop = FALSE])^2
  }))))
}
n_indep <- sum(vapply(tf_components(g_b),
                      function(co) length(attr(co, "score_variance")), integer(1)))
data.frame(
  representation = c("independent FPC", "joint MFPCA"),
  stored_scores  = c(n_indep, attr(g_m, "mfpc")$npc),
  rmse           = round(c(rmse_mv(g_b), rmse_mv(g_m)), 3)
)
#>    representation stored_scores  rmse
#> 1 independent FPC            23 0.359
#> 2     joint MFPCA             7 1.228

The first 2 shared modes already capture 69% of the joint variance. New subjects can be projected onto this fitted basis with tf_rebase(), which re-scores them jointly rather than component by component.

3. Atlantic storms as 4-dimensional curves

dplyr::storms records four time-varying quantities for every Atlantic tropical storm or hurricane from 1975 to 2024, sampled every 6 hours along the storm’s life: position (longitude, latitude) and intensity (sustained wind speed, central pressure). The natural object is a single four-component vector-valued curve per storm, f_i: [0, T_i] -> R^4 – spatial trajectory and intensity life-cycle bundled together.

Two practical points:

Code 
KM_PER_DEG <- 111.32

storms_clean <- storms |>
  mutate(
    storm_id = paste(name, year),
    ts       = as.POSIXct(ISOdate(year, month, day, hour), tz = "UTC")
  ) |>
  group_by(storm_id) |>
  distinct(ts, .keep_all = TRUE) |>      # dedupe duplicate-timestamp rows
  mutate(
    t_hours = as.numeric(ts - min(ts), units = "hours"),
    phase   = if (max(t_hours) > 0) t_hours / max(t_hours) else 0,
    ref_lat = mean(lat),
    x_km    = (long - mean(long)) * KM_PER_DEG * cos(ref_lat * pi / 180),
    y_km    = (lat  - ref_lat)    * KM_PER_DEG
  ) |>
  filter(n() >= 16, !is.na(wind), !is.na(pressure)) |>  # >= 4 days
  ungroup()

# spatial (km, real time): the compact data-frame constructor.
tracks_km <- storms_clean |>
  transmute(storm_id, t_hours, x = x_km, y = y_km) |>
  tfd_mv(id = "storm_id", arg = "t_hours", value = c("x", "y"))

# Equivalent spelling if the component tfd vectors are already available.
tracks_km_from_components <- tfd_mv(list(
  x = tfd(storms_clean, id = "storm_id", arg = "t_hours", value = "x_km"),
  y = tfd(storms_clean, id = "storm_id", arg = "t_hours", value = "y_km")
))

# full 4-d, on normalised life-cycle phase
tracks4 <- tfd_mv(list(
  long = tfd(storms_clean, id = "storm_id", arg = "phase", value = "long"),
  lat  = tfd(storms_clean, id = "storm_id", arg = "phase", value = "lat"),
  wind = tfd(storms_clean, id = "storm_id", arg = "phase", value = "wind"),
  pres = tfd(storms_clean, id = "storm_id", arg = "phase", value = "pressure")
))
tracks4
#> tfd_mv<d=4>[518] (long, lat, wind, pres): [0, 1] -> [-136.9, 13.5] x [7, 70.7] x [10, 165] x [882, 1024]
#> components based on 16 to 96 evaluations each, interpolation by tf_approx_linear
#> [1]: (0.000,-79);(0.033,-79);(0.067,-79); ... | (0.000,28);(0.033,28);(0.067,30); ... | (0.000,25);(0.033,25);(0.067,25); ... | (0.000,1013);(0.033,1013);(0.067,1013); ...
#> [2]: (0.000,-68);(0.053,-70);(0.105,-70); ... | (0.000,26);(0.053,26);(0.105,26); ... | (0.000,20);(0.053,20);(0.105,25); ... | (0.000,1014);(0.053,1014);(0.105,1014); ...
#> [3]: (0.000,-70);(0.031,-71);(0.062,-72); ... | (0.000,22);(0.031,22);(0.062,22); ... | (0.000,25);(0.031,25);(0.062,25); ... | (0.000,1011);(0.031,1011);(0.062,1010); ...
#> [4]: (0.000,-46);(0.036,-48);(0.071,-48); ... | (0.000,33);(0.036,34);(0.071,34); ... | (0.000,35);(0.036,40);(0.071,40); ... | (0.000,1005);(0.036,1005);(0.071,1005); ...
#> [5]: (0.000,-55);(0.022,-56);(0.044,-57); ... | (0.000,18);(0.022,18);(0.044,18); ... | (0.000,25);(0.022,25);(0.044,25); ... | (0.000,1009);(0.022,1009);(0.044,1009); ...
#> [6]: (0.000,-57);(0.056,-57);(0.111,-58); ... | (0.000,23);(0.056,24);(0.111,24); ... | (0.000,25);(0.056,30);(0.111,35); ... | (0.000,1005);(0.056,1005);(0.111,1005); ...
#> 
#>     [....]   (512 not shown)

# storm-level metadata
peak <- storms |>
  group_by(name, year) |>
  summarise(peak_cat = suppressWarnings(max(category, na.rm = TRUE)),
            .groups = "drop") |>
  mutate(peak_cat = if_else(is.finite(peak_cat), peak_cat, 0),
         peak_cat = as.integer(peak_cat),
         storm_id = paste(name, year))

df <- tibble::tibble(
  storm_id  = names(tracks4),
  track     = tracks4,
  track_km  = tracks_km
) |>
  left_join(peak, by = "storm_id") |>
  mutate(strength = factor(
    pmin(peak_cat, 4),
    levels = 0:4,
    labels = c("TS/TD", "Cat 1", "Cat 2", "Cat 3", "Cat 4+")
  ))

Before analysing the storms, inspect the irregular object directly. Counts vary by storm because each lifetime has a different number of 6-hourly observations, while the long data-frame representation keeps the components aligned by (storm, phase).

Code 
head(tf_count(tracks4))
#>               long lat wind pres
#> Amy 1975        31  31   31   31
#> Blanche 1975    20  20   20   20
#> Caroline 1975   33  33   33   33
#> Doris 1975      29  29   29   29
#> Eloise 1975     46  46   46   46
#> Faye 1975       19  19   19   19
as.data.frame(tracks4[1:2], unnest = TRUE) |> head()
#>         id        arg component  value
#> 1 Amy 1975 0.00000000      long  -79.0
#> 2 Amy 1975 0.00000000       lat   27.5
#> 3 Amy 1975 0.00000000      wind   25.0
#> 4 Amy 1975 0.00000000      pres 1013.0
#> 5 Amy 1975 0.03333333      long  -79.0
#> 6 Amy 1975 0.03333333       lat   28.5

Movement workflow: regularize, then describe

Joo et al. (2019) describe movement analysis as a workflow built around tracking records (x, y, t): clean the fixes, regularize or reconstruct the path when needed, visualize the track, then extract descriptors such as speed, heading and turning angles. The storm data follow the same pattern. The build step above already performed basic preprocessing: duplicate timestamps were removed, very short tracks were dropped, and longitude/latitude were projected into local kilometre coordinates.

For diagnostic plots it is useful to regularize a few storms onto a common lifecycle grid. Differentiating the regularized (x, y) curves gives velocity components; their norm is forward speed, and atan2(v_y, v_x) gives heading. Turning angle is the wrapped difference between successive headings.

Code 
movement_ids <- c("Andrew 1992", "Katrina 2005", "Sandy 2012")
movement_ids <- movement_ids[movement_ids %in% df$storm_id]
movement_rows <- match(movement_ids, df$storm_id)
movement_grid <- seq(0, 1, length.out = 81)
movement_duration <- vapply(
  tf_arg(df$track_km[movement_rows]),
  \(t) max(t) - min(t),
  numeric(1)
)

movement_eval <- lapply(seq_along(movement_rows), function(k) {
  tf_evaluate(
    df$track_km[movement_rows[k]],
    arg = movement_grid * movement_duration[k]
  )[[1]]
})
movement_x <- do.call(rbind, lapply(movement_eval, `[[`, "x"))
movement_y <- do.call(rbind, lapply(movement_eval, `[[`, "y"))

movement_track <- tfd_mv(list(
  x = tfd(movement_x, arg = movement_grid, domain = c(0, 1)),
  y = tfd(movement_y, arg = movement_grid, domain = c(0, 1))
), domain = c(0, 1))
names(movement_track) <- movement_ids

movement_velocity <- tf_derive(movement_track)
movement_speed <- tf_norm(movement_velocity) / movement_duration

vx <- as.matrix(tf_component(movement_velocity, "x"),
                arg = movement_grid, interpolate = TRUE)
vy <- as.matrix(tf_component(movement_velocity, "y"),
                arg = movement_grid, interpolate = TRUE)
heading_mat <- atan2(vy, vx) * 180 / pi
heading <- tfd(heading_mat, arg = movement_grid, domain = c(0, 1))
names(heading) <- movement_ids

wrap_degrees <- function(x) ((x + 180) %% 360) - 180
turning_mat <- t(apply(heading_mat, 1, function(x) {
  c(NA_real_, wrap_degrees(diff(x)))
}))
turning <- tfd(turning_mat, arg = movement_grid, domain = c(0, 1))
names(turning) <- movement_ids

movement_cols <- c("#1b9e77", "#d95f02", "#7570b3")[seq_along(movement_ids)]
xy <- as.matrix(movement_track, interpolate = TRUE)

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(range(xy[, , "x"], na.rm = TRUE), range(xy[, , "y"], na.rm = TRUE),
     type = "n", asp = 1, xlab = "x (km)", ylab = "y (km)",
     main = "regularized tracks")
lines(movement_track, col = movement_cols, lwd = 2)
legend("topleft", bty = "n", lwd = 2, col = movement_cols,
       legend = movement_ids, cex = 0.85)
plot(movement_speed, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "km / h", main = "forward speed")
plot(heading, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "degrees", main = "heading")
plot(turning, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "degrees", main = "turning angle")

Code 
par(op)

The result is a movement-data diagnostic rather than just a map: storms with visually similar tracks can differ sharply in when they accelerate, when they turn, and how abruptly their heading changes during recurvature.

A single 4-d curve: Hurricane Katrina (2005)

For a single storm the 4-component tf_mv is naturally displayed in facet mode (type = "facet") – one panel per component, all on the same normalised cycle phase axis:

Code 
katrina <- df |> filter(storm_id == "Katrina 2005") |> pull(track)
plot(katrina, type = "facet", lwd = 2)

Wind and pressure are anti-correlated – pressure dips below 920 mbar as wind peaks above 150 knots near phase ~0.55 – and the (long, lat) panels show Katrina sliding northwest across the Gulf and turning north into Louisiana. Bundling intensity and trajectory as a single tfd_mv lets every operation downstream (subsetting, summaries, arithmetic, plotting) treat them as one object.

Track map, faceted by peak intensity

The geographic view uses the (long, lat) components only. Faceting by peak Saffir-Simpson category separates the populations cleanly, and a coastline backdrop (from maps::map, if installed) anchors the geography:

Code 
have_maps <- requireNamespace("maps", quietly = TRUE)
pal <- c("grey50", "#fed976", "#feb24c", "#fd8d3c", "#e31a1c")

draw_coast <- function(xlim, ylim) {
  if (!have_maps) return(invisible())
  m <- maps::map("world", plot = FALSE,
                 xlim = xlim, ylim = ylim, fill = FALSE)
  lines(m$x, m$y, col = "grey55", lwd = 0.6)
}
Code 
xlim <- range(unlist(tf_evaluations(tf_component(df$track, "long"))), na.rm = TRUE)
ylim <- range(unlist(tf_evaluations(tf_component(df$track, "lat"))),  na.rm = TRUE)

op <- par(mfrow = c(2, 3), mar = c(3.4, 3.4, 2, 1), mgp = c(2, 0.7, 0))
for (k in seq_along(levels(df$strength))) {
  lev   <- levels(df$strength)[k]
  trks  <- df$track[df$strength == lev]
  spatial <- tfd_mv(list(
    long = tf_component(trks, "long"),
    lat  = tf_component(trks, "lat")
  ))
  plot(range(xlim), range(ylim), type = "n",
       xlab = "long", ylab = "lat",
       main = sprintf("%s (n = %d)", lev, length(trks)))
  draw_coast(xlim, ylim)
  lines(spatial, col = pal[k], alpha = 0.4, lwd = 1)
}
plot.new()
legend("center", bty = "n", lwd = 2, col = pal, legend = levels(df$strength),
       title = "peak intensity", cex = 1.05)

Code 
par(op)

Tropical storms / depressions stay in the southern, western basin; cat-3 and cat-4+ tracks reach further north and east – they live longer, get caught by the westerlies, and recurve out to sea.

Scalar features extracted from the 4-d object

Every component of a tf_mv is a real tf vector and supports the full univariate API, so per-storm scalar features fall out as one-liners:

Code 
df <- df |> mutate(
  path_km    = tf_arclength(track_km),
  duration   = vapply(tf_arg(track_km), \(t) max(t) - min(t), numeric(1)),
  mean_speed = path_km / duration,            # km/h, lifetime average
  peak_wind  = vapply(tf_evaluations(tf_component(track, "wind")), max, numeric(1)),
  min_pres   = vapply(tf_evaluations(tf_component(track, "pres")), min, numeric(1))
)

df |>
  group_by(strength) |>
  summarise(
    n = dplyr::n(),
    median_path_km     = round(median(path_km)),
    median_speed_kmh   = round(median(mean_speed), 1),
    median_peak_wind   = median(peak_wind),
    median_min_pres    = median(min_pres),
    .groups = "drop"
  )
#> # A tibble: 5 × 6
#>   strength     n median_path_km median_speed_kmh median_peak_wind
#>   <fct>    <int>          <dbl>            <dbl>            <dbl>
#> 1 TS/TD      201           2924             20.8               50
#> 2 Cat 1      117           4266             21.4               75
#> 3 Cat 2       62           4860             22.2               90
#> 4 Cat 3       52           6035             22.9              105
#> 5 Cat 4+      86           7399             24.6              130
#> # ℹ 1 more variable: median_min_pres <dbl>

Median path length grows strongly from tropical-storm to cat-4+ storms, whereas lifetime-average forward speed increases more modestly. The canonical wind / pressure gap between categories is visible.

Intensity and forward-speed life-cycles per category

Plotting all forward-speed and intensity time courses at once is hopeless (hundreds of irregular noodles). What we actually want is the mean curve per intensity stratum on the normalised lifecycle phase axis. Evaluate each component on a common phase grid, average within each stratum, and re-wrap the result as a length-G univariate tfd so we can use the standard plot.tf / lines.tf machinery:

Code 
phase_grid <- seq(0, 1, length.out = 41)

# forward speed on normalised time: re-arg each storm's km-speed by
# t / T_i, so all storms share phase domain [0, 1].
speed_km  <- tf_speed(df$track_km)
#> ✖ Differentiating over irregular grids can be unstable.
#> ✖ Differentiating over irregular grids can be unstable.
durations <- df$duration
speed_long <- do.call(rbind, lapply(seq_along(speed_km), function(i) {
  data.frame(
    id    = names(speed_km)[i],
    phase = tf_arg(speed_km[i])[[1]] / durations[i],
    value = tf_evaluations(speed_km[i])[[1]]
  )
}))
speed_phase <- tfd(speed_long, id = "id", arg = "phase", value = "value",
                   domain = c(0, 1))

# per-stratum mean curve on the regular phase grid, packaged as length-G tfd
stratum_mean_tfd <- function(comp, grp, grid = phase_grid) {
  mat <- as.matrix(comp, arg = grid, interpolate = TRUE)
  means <- vapply(levels(grp), function(g) {
    rows <- which(grp == g)
    if (length(rows) < 2L) return(rep(NA_real_, length(grid)))
    colMeans(mat[rows, , drop = FALSE], na.rm = TRUE)
  }, numeric(length(grid)))
  out <- tfd(t(means), arg = grid, domain = c(0, 1))
  names(out) <- levels(grp)
  out
}

wind_avg  <- stratum_mean_tfd(tf_component(df$track, "wind"), df$strength)
pres_avg  <- stratum_mean_tfd(tf_component(df$track, "pres"), df$strength)
speed_avg <- stratum_mean_tfd(speed_phase,                     df$strength)
Code 
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(wind_avg, col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "wind (knots)",
     main = "mean sustained wind")
legend("topright", bty = "n", lwd = 2, col = pal, legend = levels(df$strength),
       cex = 0.9, title = "peak intensity")
plot(pres_avg,  col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "pressure (mbar)",
     main = "mean central pressure")
plot(speed_avg, col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "forward speed (km/h)",
     main = "mean forward speed")
plot.new()
text(0.5, 0.7, "wind and pressure peak\nnear lifecycle phase 0.5;",
     adj = 0.5, cex = 1.05)
text(0.5, 0.35, "forward speed peaks\nlate, during recurvature",
     adj = 0.5, cex = 1.05)

Code 
par(op)

Three clean stories on a single tf_mv object: intensity peaks mid-life-cycle (around phase ~ 0.5) and the gap between categories is roughly uniform across the cycle; central pressure mirrors wind; forward speed peaks much later (around phase ~ 0.75) and the strongest storms accelerate the most – the canonical recurvature signature.

Smoothing on normalised time

Fitting a basis representation on raw real-time data has a tractable domain-extrapolation problem – the longest storm runs 600+ hours, but most tracks end after 100-200 hours, and a shared basis over [0, 600] extrapolates wildly. On the normalised phase version this just disappears, because every storm spans [0, 1]:

Code 
top6 <- df |> arrange(desc(path_km)) |> slice(1:6) |> pull(storm_id)

# fit a per-component spline basis on (long, lat, wind, pres) over [0, 1]
tb <- tfb_mv(df$track[df$storm_id %in% top6], k = 12, verbose = FALSE)

# pull just the (long, lat) components out of the 4-d objects so plot.tf_mv
# defaults to the trajectory (long, lat) view
raw_xy <- df$track[df$storm_id %in% top6, , c("long", "lat")]
sm_xy  <- tb[, , c("long", "lat")]

op <- par(mfrow = c(1, 2), mar = c(3.6, 3.6, 2, 1), mgp = c(2, 0.7, 0))
plot(range(xlim), range(ylim), type = "n",
     xlab = "long", ylab = "lat", main = "raw observations")
draw_coast(xlim, ylim)
lines(raw_xy, col = 1:6, lwd = 1.6)
plot(range(xlim), range(ylim), type = "n",
     xlab = "long", ylab = "lat", main = "spline-smoothed")
draw_coast(xlim, ylim)
lines(sm_xy, col = 1:6, lwd = 1.6)

Code 
par(op)

The smoothed trajectories pick up the gross recurvature shape without the 6-hourly sampling jitter and – thanks to the normalised time axis – without any extrapolation pathology. The same tb object also carries clean smoothed wind / pressure curves for each storm; for example:

Code 
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(tf_component(tb, "wind"), col = 1:6, lwd = 2,
     xlab = "lifecycle phase", ylab = "wind (knots)",
     main = "smoothed wind")
plot(tf_component(tb, "pres"), col = 1:6, lwd = 2,
     xlab = "lifecycle phase", ylab = "pressure (mbar)",
     main = "smoothed pressure")

Code 
par(op)

4. Recap

The two case studies above exercised exactly the same surface:

The common pattern is simple: keep coupled signals bundled as one object, use component accessors when a scalar or univariate summary is needed, and let the multivariate methods reuse the existing univariate kernels component-wise.

References

Happ, Clara, and Sonja Greven. 2018. “Multivariate Functional Principal Component Analysis for Data Observed on Different (Dimensional) Domains.” Journal of the American Statistical Association 113 (522): 649–59. https://doi.org/10.1080/01621459.2016.1273115.
Joo, Rocío, Matthew E. Boone, Thomas A. Clay, Samantha C. Patrick, Susana Clusella-Trullas, and Mathieu Basille. 2019. “Navigating Through the r Packages for Movement.” Journal of Animal Ecology 89 (1): 248–67. https://doi.org/10.1111/1365-2656.13116.
Marron, J. S., James O. Ramsay, Laura M. Sangalli, and Anuj Srivastava. 2015. “Functional Data Analysis of Amplitude and Phase Variation.” Statistical Science 30 (4): 468–84. https://doi.org/10.1214/15-STS524.
Olshen, Richard A., Edmund N. Biden, Marilynn P. Wyatt, and David H. Sutherland. 1989. “Gait Analysis and the Bootstrap.” The Annals of Statistics 17 (4): 1419–40. https://doi.org/10.1214/aos/1176347373.
Ramsay, James O., and Bernard W. Silverman. 2005. Functional Data Analysis. 2nd ed. Springer Series in Statistics. New York: Springer.
Srivastava, Anuj, and Eric P. Klassen. 2016. Functional and Shape Data Analysis. Springer Series in Statistics. New York: Springer. https://doi.org/10.1007/978-1-4939-4020-2.
Srivastava, Anuj, Eric Klassen, Shantanu H. Joshi, and Ian H. Jermyn. 2011. “Shape Analysis of Elastic Curves in Euclidean Spaces.” IEEE Transactions on Pattern Analysis and Machine Intelligence 33 (7): 1415–28. https://doi.org/10.1109/TPAMI.2010.184.
Srivastava, Anuj, Wei Wu, Sebastian Kurtek, Eric Klassen, and J. S. Marron. 2011. “Registration of Functional Data Using Fisher-Rao Metric.” arXiv Preprint arXiv:1103.3817. https://arxiv.org/abs/1103.3817.
Tucker, J. Derek, Wei Wu, and Anuj Srivastava. 2013. “Generative Models for Functional Data Using Phase and Amplitude Separation.” Computational Statistics & Data Analysis 61: 50–66. https://doi.org/10.1016/j.csda.2012.12.001.