Analysis of Uni- and Multi-Dimensional Contours with GAMs and Functional PCA: an Application to Ultrasound Tongue Imaging
Author
Affiliation
Michele Gubian
Institute for Phonetics and Speech Processing, Ludwig Maximilian University, Munich
Published
April 24, 2025
Loading libraries and a few helper scripts.
Code
# data processing and plottinglibrary(tidyverse)# GAMslibrary(mgcv)library(itsadug)# FPCAlibrary(funData)library(MFPCA)# install.packages("devtools")# devtools::install_github("uasolo/landmarkregUtils")library(landmarkregUtils)library(abind)# LMERlibrary(lme4)library(emmeans)# plottinglibrary(RColorBrewer)library(gridExtra)# the following scripts are available at: https://github.com/uasolo/UTI-GAM-FPCAsource("polar_utils.R")# artificial UTI data generatorsource("parab_tongue.R")# cosmeticsCategory.colors <-c(A ="black", B ="limegreen")theme_set(theme_light())
1 Static analysis
Static refers to absence of time. The typical scenario is the analysis of UTI where only one representative contour is used for each speech sound, e.g. the onset or the point in time where the articulatory target is reached.
1.1 Data generation
This tutorial is based on artificial data. The helper function parab_tongue() generates curves that resemble tongue contours. The shape consists of two branches of parabolas joined at the vertex. The shape is parametrised so that it can vary in various ways. The output is a tibble with three columns: knot, x and y, where knot simulates the (11) knots obtained with DeepLabCut (Wrench & Balch-Tomes 2022). See the source code parab_tongue.R for more detail.
The \((x, y)\) coordinates are usually not suitable for a classic \(y = f(x)\) model, as the tongue contour in some regions has more than one \(y\) value corresponding to the same \(x\), which is incompatible with a standard function, e.g. with a GAM y ~ s(x). The polar coordinate representation mitigates this problem, as most of the time \(radius = f(angle)\) is well behaved.
In order to transform \((x, y)\) coordinates into polar ones we need to specify an origin point. Here the location of the origin will be chosen by eye inspection. The helper function cart2polar() performs the transformation, while radial_plot() is a plotting helper.
The dataset for the static analysis consists of 60 artificial tongue contours, 30 for each Category, A and B. This is the simplest possible scenario, with only one categorical predictor with two levels and no random factors. Hints on adaptations for more complex and realistic scenarios will be provided throughout the exposition.
Here we fit a GAM to the dataset, using the polar representation. This means that we predict radius given angle and Category. As angle is a continuous variable, this will be modelled with smooth terms s().
GAM fitting is carried out with a simple, standard procedure. Essentially, we first fit a minimal model without AR1 (autoregressive) term, from which we estimate an AR1 term that is used for a second fit. The rest of the model parameters are set using default values. A scaled-t residual family is used whenever possible. Model fitting, diagnostics and inference are generally based on Wieling (2018), to which the reader is referred for a general introduction on GAMs. For more complex predictor structure, including random terms, see Wieling (2018) and the tutorial materials on this GitHub.
Method: fREML Optimizer: perf chol
$grad
[1] -1.339989e-07 -1.995888e-06
$hess
[,1] [,2]
[1,] 2.826731e+00 -2.113013e-18
[2,] 4.882741e-31 2.598469e+00
Model rank = 20 / 20
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(angle):CategoryA 9.00 6.78 0.87 0.005 **
s(angle):CategoryB 9.00 7.09 0.87 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf_resid(GAM)
angle_grid_plot <-seq(min(curves$angle), max(curves$angle), length.out =20)pred <-get_predictions(GAM,cond =list(angle = angle_grid_plot, Category =c("A", "B")),print.summary =FALSE)gam_polar_plot <-ggplot(pred) +aes(x = angle, y = fit, group = Category ) +geom_line(mapping =aes(color = Category)) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +geom_ribbon(mapping =aes(ymin = fit - CI, ymax = fit + CI, fill = Category), alpha =0.5) +radial_plot(60) +theme(legend.position ="bottom")
Code
gam_polar_plot
Figure 3: GAM predicted contours
The contour shapes capture the characteristics of the original dataset fairly well. A few observations:
The original data presents high variation at the left (tongue root), which is not visible in the predicted contours. This is because the confidence bands reflect the uncertainty in the estimation of the mean contours, not the extra unexplained noise (the same principle applies to ordinary linear regression estimates, e.g. lm())
The right extreme of the B contour has a wide confidence band. This is because the B contours end early at the right (tongue tip) in the original dataset, hence there are not enough data points to perform accurate estimations.
To overcome the second issue, two simple approaches can be applied:
linear angle normalization
Procrustean fit
where Procrustean means that we establish a fixed angle range for all contours, cut off the points that are beyond that range and extrapolate points for contours that do not have values up to the range limits.
Both approaches are quite simple, yet can introduce problems. Cutting extremes eliminates information, e.g. part of the tongue tip is removed from some contours. Linear angle normalization distorts the shape of the contour.
1.4.1 Linear angle normalization
A simple way to establish the reference angle range is to take the median angle of the first/last point from all contours. The contours are then reinterpolated on a regular grid within the common angle range. This of course introduces approximations, including the fact that now knots lose their original interpretation. However, a regular grid offers computational advantages, especially for FPCA (see below).
Below we apply linear angle normalization, then re-fit the GAM in the same way as above.
Linear angle normalization
angle_range <- raw_curves %>%filter(knot %in%range(knot)) %>%# first and last knotgroup_by(knot) %>%summarise(median_angle =median(angle)) %>%pull(median_angle) %>%sort()angle_grid <-seq(angle_range[1], angle_range[2], length.out=11)curves <- raw_curves %>%group_by(Category, frame_id) %>%# linear angle normalization on the angle_range intervalmutate(angle = (angle -min(angle)) *diff(angle_range) /diff(range(angle)) + angle_range[1]) %>%# linear interpolation on angle_grid (optional for GAM) reframe(bind_cols(knot =seq_along(angle_grid), # fake knots, needed for AR1 term in GAMapprox(angle, radius, angle_grid) %>%as_tibble() )) %>%rename(angle = x, radius = y)
radius ~ angle: GAM fit linear angle normalization
Method: fREML Optimizer: perf chol
$grad
[1] 9.204291e-10 2.899667e-09
$hess
[,1] [,2]
[1,] 3.088541e+00 -1.197073e-18
[2,] -1.503855e-19 3.610893e+00
Model rank = 20 / 20
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(angle):CategoryA 9.00 6.68 0.95 0.12
s(angle):CategoryB 9.00 7.82 0.95 0.15
acf_resid(GAM)
angle_grid_plot <- angle_gridpred <-get_predictions(GAM,cond =list(angle = angle_grid_plot, Category =c("A", "B")),print.summary =FALSE)gam_polar_lin_plot <-ggplot(pred) +aes(x = angle, y = fit, group = Category ) +geom_line(mapping =aes(color = Category)) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +geom_ribbon(mapping =aes(ymin = fit - CI, ymax = fit + CI, fill = Category), alpha =0.5) +radial_plot(60) +theme(legend.position ="bottom")
Code
gam_polar_lin_plot
Figure 4: GAM predicted contours, linear angle normalization
In this case, Category B’s tongue tip reaches the same angle as A, which is not the case in the real data set. This is an artefact introduced by the linear stretching along the angle axis.
1.4.2 Procrustean angle range
As an alternative to linear angle normalization, we apply a Procrustean angle interpolation using the same angle range estimated above.
Procrustean angle range
curves <- raw_curves %>%group_by(Category, frame_id) %>%# interpolationreframe(bind_cols(knot =seq_along(angle_grid), # fake knots, needed for AR1 term in GAMapprox(angle, radius, angle_grid, rule =2) %>%as_tibble() )) %>%rename(angle = x,radius = y)
Method: fREML Optimizer: perf chol
$grad
[1] -1.423750e-12 7.857182e-11
$hess
[,1] [,2]
[1,] 2.823996e+00 -3.510719e-19
[2,] 9.333367e-20 3.211993e+00
Model rank = 20 / 20
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(angle):CategoryA 9.00 6.29 0.94 0.020 *
s(angle):CategoryB 9.00 7.76 0.94 0.035 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf_resid(GAM)
angle_grid_plot <- angle_gridpred <-get_predictions(GAM,cond =list(angle = angle_grid_plot, Category =c("A", "B")),print.summary =FALSE)gam_polar_proc_plot <-ggplot(pred) +aes(x = angle, y = fit, group = Category ) +geom_line(mapping =aes(color = Category)) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +geom_ribbon(mapping =aes(ymin = fit - CI, ymax = fit + CI, fill = Category), alpha =0.5) +radial_plot(60) +theme(legend.position ="bottom")
Code
gam_polar_proc_plot
Figure 5: GAM predicted contours, Procrustean fit
Some of the information related to the tongue tip has been lost. Perhaps this is the reason why the GAM’s explained deviance is reduced. However, the predicted contours have a closer resemblance to the original dataset than in the linear angle normalization version, which introduced distortions in the contour shape.
In the rest of the tutorial, we will always use the Procrustean method. It is up to the user to choose the method most appropriate for the data at hand.
1.4.3 Testing for significance: binary smooths
Here we aim at showing if and where on the angle axis the difference between CategoryA and B is significant, and in which direction it goes. A graphical command is available for this (itsadug::plot_diff). Here instead we use a more general technique based on binary smooths (Wieling 2018), which consists in explicitly defining a curve as the difference between two given factor levels (in our case B - A), and then use its confidence interval to determine whether the difference is significantly different from zero. This approach can be applied also in more complex inference situations, e.g. difference of differences, see Wieling (2018) and GAMInteractions for more detail.
Binary smooths
# define a binary variable to indicate class B (as in Wieling 2018)curves <- curves %>%mutate(IsB = (Category =='B') *1)# the GAM equation contains a difference smooth representing B - Agam_formula <-"radius ~ s(angle) + s(angle, by = IsB)"GAM <-bam(gam_formula %>%as.formula(),discrete =TRUE,family ='scat',rho = AR1.rho,AR.start = curves %>%pull(knot) ==1,data = curves)summary(GAM)
diff_plot_GAM_angle <-get_modelterm(GAM, select =2, print.summary =FALSE) %>%# the second smooth term in the summaryggplot() +aes(x = angle, y = fit) +geom_line() +geom_ribbon(mapping =aes(ymin = fit - se.fit, ymax = fit + se.fit), fill ="grey", alpha =0.5) +scale_x_reverse() +geom_hline(yintercept=0, linetype="dashed", color ="red", linewidth=1.5) +ylab("Radius difference")
Code
diff_plot_GAM_angle
Figure 6: Difference curve B - A
To plot the difference contour B - A we have abandoned the radial plot, as differences in radius can be negative, but a negative radius cannot be represented in polar coordinates.
1.5 radius ~ angle: FPCA
Here we apply FPCA followed by linear regression to the same dataset we use to estimate GAMs. This procedure is an alternative to GAMs. The procedure is in two steps: First a functional PCA model is estimated on the basis of the dataset, then the model parameters, called scores, are used as response variable in a linear regression model, i.e. lm(score ~ Category). The FPCA model describes the independent variations in the entire curve dataset without taking into account the categorical predictors. This provides an overview of all the sources of variation in the data. Then it is our task to identify those sources that are correlated to the relevant predictors (here Category). For more information, see Gubian, Torreira & Boves (2015) and the tutorial materials on this GitHub.
In the code that follows, we use the the same curves as the last radius ~ angle GAM example, obtained by Procrustean fit. The code is based on the ex1D.R script.
radius ~ angle: FPCA fit
# build a funData objectcurvesFun <-long2irregFunData(curves, id ="frame_id", time ="angle", value ="radius") %>%as.funData()curvesFun
Functional data with 60 observations of 1-dimensional support
argvals:
1.41687 1.516304 ... 2.411207 (11 sampling points)
X:
array of size 60 x 11
FPCA is computed by the function MFPCA::PACE(). As default, it computes a number of PCs covering 99% of the variance. In this case it computed 4 PCs.
We will look at the first two PCs, which explain 79.8% and 12.3% of the variance.
radius ~ angle: FPCA curves
nPC <-2# scores st. dev.sdFun <- fpca$values %>%sqrt()# PC curves to be plottedPCcurves <-expand_grid(PC =1:nPC,fractionOfStDev =seq(-1, 1, by=.25)) %>%group_by(PC, fractionOfStDev) %>%reframe(funData2long1(fpca$mu + fractionOfStDev * sdFun[PC] * fpca$functions[PC],time ="angle", value ="radius") )# PlotFPCA_angle_curves_plot <-ggplot(PCcurves) +aes(x = angle, y = radius, group = fractionOfStDev, color = fractionOfStDev) +geom_line() +scale_color_gradient2(low ="blue", mid ="grey", high ="orangered",breaks =c(-1, 0 , 1)) +facet_wrap(~ PC, nrow =1,labeller =labeller(PC =~str_glue("PC{.x}"))) +labs(color =expression(frac(s[k], sigma[k]))) +geom_line(data = PCcurves %>%filter(fractionOfStDev ==0), color ='black', linewidth =1.2) +radial_plot(60) +theme(legend.position ="bottom")
Code
FPCA_angle_curves_plot
Figure 7: First 2 PCs: effect of scores on contour shape
Each PC captures a different shape variation in the data set, parametrized by a score, (s1 on the left, s2 on the right). The black contour is the mean shape, color-coded contours show the relation between scores (normalized by their respective standard deviation) and the shape modelled by the corresponding PC.
In this case we see that PC1 models the large variation towards the left (tongue root), while PC2 models the smaller variation in the middle (tongue dorsum). The two variations are independent, i.e. knowing the value of s1 does not help guessing the value of s2.
This model is estimated without making use of the Category values. When we associate back the values of Category to scores, we can determine which score correlates with Category.
radius ~ angle: FPCA scores
# collect PC scoresPCscores <- fpca$scores %>%`colnames<-`( paste0("s", 1:fpca$npc)) %>%as_tibble() %>%bind_cols(curves %>%distinct(frame_id, Category), .)# boxplots PC scores by CategoryFPCA_angle_scores_plot <- PCscores %>%pivot_longer(cols = s1:all_of(str_glue("s{fpca$npc}")), names_to ="PC", values_to ="score") %>%filter(PC %in%str_c("s", 1:nPC)) %>%ggplot() +aes(x = Category, y = score, color = Category) +geom_boxplot() +facet_wrap(~ PC) +scale_color_manual(values=Category.colors) +theme(legend.position ="bottom")
Code
FPCA_angle_scores_plot
Figure 8: Distribution of the first 2 PC scores by Category
Clearly s2 is the score that correlates with Category. Note that it is s2 and not s1 the score that is more interesting to our analysis, although s1 is associated to PC1, which by definition explains most of the variance. But this variance is largely independent of the variation correlated to Category. In general, it is not guaranteed that there will be exactly one score correlated with a variable of interest: there could be none, one or more than one.
Let us estimate a simple linear model predicting s2 from Category, i.e. lm(s2 ~ Category). Note that this is the simplest possible model, as we do not have multiple predictor or random terms. In realistic cases where these are present, we have to apply a correspondingly more complex model, e.g. using lme4::lmer(), but the idea is the same: predict a score based on predictors not used in the estimation of the FPCA model.
radius ~ angle: FPCA lm
s <-2# score indexmodel_eq <-str_glue("s{s} ~ Category") %>%as.formula()mod <-lm(model_eq, data = PCscores)mod %>%summary()
Call:
lm(formula = model_eq, data = PCscores)
Residuals:
Min 1Q Median 3Q Max
-1.47535 -0.42984 -0.02173 0.40532 2.58482
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1219 0.1386 8.097 4.22e-11 ***
CategoryB -2.4740 0.1960 -12.625 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7589 on 58 degrees of freedom
Multiple R-squared: 0.7332, Adjusted R-squared: 0.7286
F-statistic: 159.4 on 1 and 58 DF, p-value: < 2.2e-16
emmeans(mod, pairwise ~ Category)
$emmeans
Category emmean SE df lower.CL upper.CL
A 1.12 0.139 58 0.845 1.40
B -1.35 0.139 58 -1.629 -1.07
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B 2.47 0.196 58 12.625 <.0001
# reconstruct predicted curvespredCurves <-emmeans(mod, pairwise ~ Category)$emmeans %>%as_tibble() %>%group_by(Category) %>%reframe(bind_cols(funData2long1(fpca$mu + emmean * fpca$functions[s], time ='angle', value ="radius"),funData2long1(fpca$mu + lower.CL * fpca$functions[s], time ='angle', value ="yl") %>%select(yl),funData2long1(fpca$mu + upper.CL * fpca$functions[s], time ='angle', value ="yu") %>%select(yu) ))FPCA_angle_pred_plot <-ggplot(predCurves) +aes(angle, radius, color = Category) +geom_line() +geom_ribbon(aes(x = angle, ymin = yl, ymax = yu, fill = Category),alpha =0.3, inherit.aes =FALSE) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +radial_plot(60) +theme(legend.position ="bottom")# diff curve B - A diffCurve <-emmeans(mod, pairwise ~ Category)$contrasts %>%as_tibble() %>%reframe(bind_cols(# emmeans computes A - B, so we need to invert the sign of all curvesfunData2long1(-estimate * fpca$functions[s], time ="angle", value ="radius"),funData2long1(-(estimate -1.96* SE) * fpca$functions[s], time ="angle", value ="radius_l") %>%select(radius_l),funData2long1(-(estimate +1.96* SE) * fpca$functions[s], time ="angle", value ="radius_u") %>%select(radius_u) ))FPCA_angle_diff_plot <-ggplot(diffCurve) +aes(angle, radius) +geom_line() +geom_ribbon(aes(x = angle, ymin = radius_l, ymax = radius_u),alpha =0.3, inherit.aes =FALSE) +geom_hline(yintercept =0, linetype ="dashed", color ='red') +scale_x_reverse() +ylab("Radius difference")
Code
FPCA_angle_pred_plot
Figure 9: Reconstructed contours by Category, based on regression model: s2 ~ Category
Code
FPCA_angle_diff_plot
Figure 10: Difference curve B - A, based on regression model: s2 ~ Category
The shape of the predicted curves is similar to those obtained with the GAM. Comparing the difference curve obtained with FPCA with the one in Figure 6 obtained with the GAM, we see that the shape is similar, while the confidence intervals are very different. The confidence interval of the GAM-based difference curve reflects local similarities, i.e. measured at the same point on the angle axis, while the confidence interval of the FPCA-based difference curve reflects the uncertainty in the estimation of the difference between levels A and B in the lm(s2 ~ Category) model modulated by the shape of the PC2 curve. It is probably not advisable to use the FPCA-based confidence intervals to perform localized inference, i.e. determine whether contours A and B differ significantly at a specific angle location, as these do not reflect local properties directly. Difference curves obtained from FPCA will not be shown further on.
1.6 (x, y) ~ knot: GAM
Here we consider the alternative geometric parametrization consisting in using the knots as independent axis and setting up a multivariate (or multiresponse) model that predicts two values, x and y, given knot. This approach is also showcased in Coretta & Sakr (2025) on real UTI data, applying both GAMs and FPCA.
1.6.1 Multivariate GAM
The canonical way to construct a multivariate GAM is to set up a system of two parallel models, one predicting x and one predicting y. This is possible using mgcv::gam().
Multivariate GAM fit
GAM <-gam(list(x ~ Category +s(knot, by = Category, k =5), y ~ Category +s(knot, by = Category, k =5) ),family=mvn(d =2),data = raw_curves)AR1.rho <-acf_resid(GAM)[2]
GAM <-gam(list(x ~ Category +s(knot, by = Category, k =5), y ~ Category +s(knot, by = Category, k =5) ),family=mvn(d =2),rho = AR1.rho,AR.start = raw_curves %>%pull(knot) ==1,data = raw_curves)summary(GAM)
Method: REML Optimizer: outer newton
full convergence after 17 iterations.
Gradient range [-0.0007176559,0.0009222604]
(score 1855.407 & scale 1).
Hessian positive definite, eigenvalue range [1.328535,1.92555].
Model rank = 23 / 23
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(knot):CategoryA 4.00 3.22 0.99 0.38
s(knot):CategoryB 4.00 3.56 0.99 0.30
s.1(knot):CategoryA 4.00 3.98 0.99 0.34
s.1(knot):CategoryB 4.00 3.97 0.99 0.29
acf_resid(GAM)
The model diagnostic shows that the residuals are clearly non-Gaussian. However, there is no remedy to this within the mgcv library, as multivariate GAMs are implemented only for the Gaussian case (family=mvn(d = 2)). Besides that, there are no clear advantages of setting up this type of model with mgcv since there are no functions that allow to estimate any multivariate inference, i.e. taking x and y into consideration jointly and not separately.
There are other R libraries that implement multivariate GAMs: mvgam and gmfamm, both Bayesian. These are not explored here.
1.6.2 Univariate trick
An alternative approach is to apply the trick used in Wieling et al. (2016). The trick consists in creating a factor that encodes the coordinate axis, i.e. factor Dimension with two levels: x and y. Then use a regular, univariate GAM where the response variable is Position, which means x or y depending on the value of the ‘predictor’ Dimension. The advantage is that the model architecture is simpler, as it avoids the complexity of the double response variable, and that non-Gaussian models are available (here we will use scat as in the previous examples). The downside is that this is a sub-optimal model in that correlations across x and y dimension are not taken into account, and that adding a factor complicates the model equation (more on how to handle multiple factors in GAMs can be found in Wieling (2018) and in GAMInteractions).
Univariate trick fit
# Create a compound factor Cat_Dim for all combos of Category and Dimensionraw_curves_long <- raw_curves %>%pivot_longer(c(x, y), names_to ='Dimension', values_to ='Position') %>%unite("Cat_Dim", c(Category, Dimension)) %>%mutate(Cat_Dim =factor(Cat_Dim)) %>%# order to capture AR along knotsarrange(frame_id, Cat_Dim, knot)GAM <-bam(Position ~ Cat_Dim +s(knot, by = Cat_Dim),data = raw_curves_long,family ='scat',discrete =TRUE)AR1.rho <-acf_resid(GAM)[2]
The predictions are faithful to the original dataset. The advantage of the multivariate parametrization, even when implementing the workaround to avoid fitting an actual multivariate GAM, is that knots (indicated in red) preserve their original position when used to fit the model, while in the polar coordinate approach the knots were re-interpolated and lost their interpretable value.
The difference curves for CategoryB - A are computed separately by dimension x and y:
Code
plot_diff(GAM, view ="knot", comp =list(Cat_Dim =c("B_x", "A_x")),print.summary =FALSE, main ="B - A, Dimension x", ylab ="x: B - A")plot_diff(GAM, view ="knot", comp =list(Cat_Dim =c("B_y", "A_y")),print.summary =FALSE, main ="B - A, Dimension y", ylab ="y: B - A")
Figure 12: Difference curves B - A, separated by dimension x, y
Figure 13: Difference curves B - A, separated by dimension x, y
Two considerations for the interpretation of results in this setting.
For a given knot, if at least one difference curve shows a significant difference, then we can say the contours are different at that knot. For example, the difference between A and B on the y axis is not significant for knots 4 to 7, while it is on the x axis. In fact, those corresponding knots are at about the same height (y), while A knots are clearly more fronted (x) than B knots, which means they are in a different location.
If two knots with a different index are close to each other, this is not revealed by the difference curves. For example, knot 9 A and knot 11 B are very close to each other, but since they have different indices this does not show in the difference curves. Compare this to the difference curve Figure 6 obtained from the angle-based approach: In the knot-based case A and B are indicated as separated towards high values of knot (tongue tip), while in the angle-based case A and B are indicated overlapping towards low values of angle (tongue tip).
1.7 (x, y) ~ knot: FPCA
Here we implement the multivariate approach with FPCA. The advantage here is that while the FPCA model is multivariate, PC scores control multiple dimensions jointly, here x and y. This means that the following linear model lm(score ~ Category) is not any more complicated than in the univariate case. The code that follows is based on the ex2D.R script.
An object of class "multiFunData"
[[1]]
Functional data with 60 observations of 1-dimensional support
argvals:
1 2 ... 11 (11 sampling points)
X:
array of size 60 x 11
[[2]]
Functional data with 60 observations of 1-dimensional support
argvals:
1 2 ... 11 (11 sampling points)
X:
array of size 60 x 11
# Compute FPCAmfpca <-MFPCA(curvesFun2D,# computing more PCs than we will use:# a quick way to get a good approx for the % explained var per PCM =5, uniExpansions =list(list(type ="uFPCA"),list(type ="uFPCA")))# Prop of explained varPropExpVar <-round(mfpca$values /sum( mfpca$values) , digits =3)
Multivariate FPCA is computed by the function MFPCA::MFPCA(). The first two PCs explain 69.1% and 23.2% of the variance.
(x, y) ~ knot: FPCA curves
nPC <-2# scores st. dev.sdFun <- mfpca$values %>%sqrt()# PC curves to be plottedPCcurves <-expand_grid(PC =1:nPC,Dim =1:2, fractionOfStDev =seq(-1, 1, by=.25)) %>%group_by(PC, Dim, fractionOfStDev) %>%reframe(funData2long1( mfpca$meanFunction[[Dim]] + fractionOfStDev * sdFun[PC] * mfpca$functions[[Dim]][PC],time="knot", value ="Pos") ) %>%mutate(Dim =factor(Dim, levels =c(1,2), labels =c('x', 'y')))# PlotFPCA_knot_curves_plot <-ggplot(PCcurves) +aes(x = knot, y = Pos, group = fractionOfStDev, color = fractionOfStDev) +geom_line() +scale_color_gradient2(low ="blue", mid ="grey", high ="orangered",breaks =c(-1, 0 , 1)) +facet_grid(Dim ~ PC,scales ="free_y",labeller =labeller(PC =~str_glue("PC{.x}"))) +labs(color =expression(frac(s[k], sigma[k]))) +theme(legend.position ="right",axis.title.y =element_blank())FPCA_knot_traj_plot <-ggplot(PCcurves %>%pivot_wider(names_from = Dim, values_from = Pos)) +aes(x = x, y = y, group = fractionOfStDev, color = fractionOfStDev) +geom_path() +geom_point(color ='black') +scale_color_gradient2(low ="blue", mid ="grey", high ="orangered",breaks =c(-1, 0 , 1)) +facet_grid(. ~ PC,# scales = "free_y",labeller =labeller(PC =~str_glue("PC{.x}"))) +labs(color =expression(frac(s[k], sigma[k]))) +theme(legend.position ="right")
Code
FPCA_knot_curves_plot
Figure 14: First 2 PCs: effect of scores on contour shape, separated by dimension x, y
Code
FPCA_knot_traj_plot
Figure 15: First 2 PCs: effect of scores on contour shape, (x,y) trajectory version
The first plot is an explicit representation of the (x, y) ~ knot model, while the second one eliminates the knot axis, with knots marked as black points.
Comparing this result with the one obtained with the angle-based approach (Figure 7), we see that the first two PCs are similar but swapped, this time the variation in the middle part of the contour explains more variance than the variation on the left. It is reassuring that both geometrical parametrizations lead to the same qualitative conclusion that these two types of variation are present in the data and that they are independent. The swapping is probably a consequence of the different parametrization, which results in a different computation of variance.
Given what we know from the angle-based analysis, we expect now that PC1 is the component that describes the variation connected to Category, which means that s1 is the score that will be used in the linear regression.
(x, y) ~ knot: FPCA lm and predictions
# collect PC scoresPCscores <- mfpca$scores[, 1:nPC] %>%`colnames<-`( paste0("s", 1:nPC)) %>%as_tibble() %>%bind_cols(raw_curves %>%distinct(frame_id, Category), .)# boxplots PC scores by CategoryFPCA_knot_scores_plot <- PCscores %>%pivot_longer(cols =paste0("s", 1:nPC), names_to ="PC", values_to ="score") %>%ggplot() +aes(x = Category, y = score, color = Category) +geom_boxplot() +facet_wrap(~ PC) +scale_color_manual(values=Category.colors) +theme(legend.position ="bottom")s <-1# score indexmodel_eq <-str_glue("s{s} ~ Category") %>%as.formula()mod <-lm(model_eq, data = PCscores)mod %>%summary()
Call:
lm(formula = model_eq, data = PCscores)
Residuals:
Min 1Q Median 3Q Max
-9.8873 -3.3157 -0.5061 3.4827 8.8519
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.4100 0.8715 -15.39 <2e-16 ***
CategoryB 27.2800 1.2325 22.13 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.773 on 58 degrees of freedom
Multiple R-squared: 0.8941, Adjusted R-squared: 0.8923
F-statistic: 489.9 on 1 and 58 DF, p-value: < 2.2e-16
emmeans(mod, pairwise ~ Category)
$emmeans
Category emmean SE df lower.CL upper.CL
A -13.4 0.872 58 -15.2 -11.7
B 13.9 0.872 58 12.1 15.6
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B -27.3 1.23 58 -22.134 <.0001
# reconstruct predicted curvespredCurves <-emmeans(mod, pairwise ~ Category)$emmeans %>%as_tibble() %>%expand_grid(Dim =1:2) %>%group_by(Category, Dim) %>%reframe(bind_cols(funData2long1(mfpca$meanFunction[[Dim]] + emmean * mfpca$functions[[Dim]][s], time ="knot", value ="Pos"),funData2long1(mfpca$meanFunction[[Dim]] + lower.CL * mfpca$functions[[Dim]][s], time ="knot", value ="Pos_l") %>%select(Pos_l),funData2long1(mfpca$meanFunction[[Dim]] + upper.CL * mfpca$functions[[Dim]][s], time ="knot", value ="Pos_u") %>%select(Pos_u) )) %>%mutate(Dim =factor(Dim, levels =c(1,2), labels =c('x', 'y'))) %>%pivot_wider(names_from = Dim, values_from =starts_with("Pos"))FPCA_knot_pred_plot <-ggplot(predCurves) +aes(group = Category, color = Category) +geom_path(aes(x = Pos_x, y = Pos_y)) +geom_path(aes(x = Pos_l_x, y = Pos_l_y), linetype ="dotted") +geom_path(aes(x = Pos_u_x, y = Pos_u_y), linetype ="dotted") +scale_color_manual(values=Category.colors) +xlab("x") +ylab("y") +theme(legend.position ="bottom")
Code
FPCA_knot_scores_plot
Figure 16: Distribution of the first 2 PC scores by Category
Code
FPCA_knot_pred_plot
Figure 17: Reconstructed contours by Category, based on regression model: s1 ~ Category
2 Dynamic analysis
With a dynamic analysis, we are adding one extra independent axis to the picture, the time axis. All considerations regarding the difference between the angle-based and the knot-based geometrical parametrization remain the same. What changes is that the model support becomes two-dimensional: angle * time or knot * time. This means that data are organized on a two-dimensional grid, while in the static analysis data were organized on a uni-dimensional grid along angle or knot. The response variables are the same as in the static case. In the angle-based case, the response variable is radius, while in the knot-based case, the response variables are x and y, i.e. a multivariate model.
In the following, we will focus mainly on the practical aspects concerning the change in the support dimensionality.
2.1 Dataset
The dataset for the dynamic analysis consists of 60 artificial tongue contour sequences, or clips, 30 for each Category, A and B. Each clip is sampled along the normalized time axis normTime on a 6-points regular grid : 0, 0.2, 0.4, 0.6, 0.8, 1.0. Below two example clips, one for each Category.
Figure 18: Dynamic dataset, one example clip per Category
The trends are clear: the tongue tip moves downward in Category A and upward in B, while the tongue root is tilted differently in the two categories but it does not change in time. These easy patterns are useful to get familiar with the complexity of the models that will be built in this section.
2.2 radius ~ angle * time: GAM
Like in the static case, we compute polar coordinates and we apply Procrustean interpolation to the angle axis. The time axis (normTime) is sampled at regular intervals and it does not need any treatment. In a realistic scenario, one may need to carry out some form of normalization or even time warping to the time axis as well, prior starting the analysis.
Below two example clips, one for each Category, after computation of polar coordinates. The first plot is a representation faithful to the data form, i.e. at each location on the two-dimensional grid support (angle * normTime) corresponds one response value (radius). The second plot (polar coordinates) reflects the geometry of the problem.
Figure 19: Dynamic dataset, one example clip per Category, flat version of polar coordinates
Code
dataset_clip_ex_polar
Figure 20: Dynamic dataset, one example clip per Category, polar coordinates
The GAM we build is similar to the one in Section 1.4, the only difference being that the uni-dimensional support smooth term s(angle, by = Category) is substituted by a two-dimensional support smooth term te(angle, normTime, by = Category).
Let us take two sections, one at a fixed angle (\(\pi/2\), around tongue tip) along time, another at a fixed time (at 0.5, half way through the gesture) along the angle axis.
Cross-sections
sec_angle <- pi/2pred <-get_predictions(GAM,cond =list(angle = sec_angle,normTime = normTime_grid,Category =c("A", "B")),print.summary =FALSE)sec_angle_plot <-ggplot(pred) +aes(x = normTime, y = fit, group = Category ) +geom_line(mapping =aes(color = Category)) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +geom_ribbon(mapping =aes(ymin = fit - CI, ymax = fit + CI, fill = Category), alpha =0.5) +theme_light() +ylab("Radius") +theme(legend.position ="bottom")sec_time <-0.5pred <-get_predictions(GAM,cond =list(angle = angle_grid_plot,normTime = sec_time,Category =c("A", "B")),print.summary =FALSE)sec_time_plot <-ggplot(pred) +aes(x = angle, y = fit, group = Category ) +geom_line(mapping =aes(color = Category)) +scale_color_manual(values=Category.colors) +scale_fill_manual(values=Category.colors) +geom_ribbon(mapping =aes(ymin = fit - CI, ymax = fit + CI, fill = Category), alpha =0.5) +theme_light() +radial_plot(60) +theme(legend.position ="bottom")
Code
sec_angle_plot
Figure 23: Cross-section of GAM predictions at angle \(\pi/2\)
Code
sec_time_plot
Figure 24: Cross-section of GAM predictions at normalized time 0.5
2.2.2 Testing for significance: binary surfaces
Similarly to the static case (Section 1.4.3), we use binary smooths, in this case a binary surface, to quantify statistical significance between CategoryA and B. Below we plot the difference surface, where the dashed contours indicate where the difference is significant.
The plot shows that B is always lower than A on the left, i.e. tongue root has a smaller radius for B than for A, while on the right (tongue tip) the radius difference varies along time, according to the patterns in Figure 21 and Figure 22.
2.3 radius ~ angle * time: FPCA
The same geometrical setup as in Section 2.2 is applied here to PCA followed by linear regression. As it is always the case, the FPCA model changes with the geometrical setup, while linear regression remains in the geometry-agnostic form lm(score ~ Category).
The FPCA code consists in first organising the dataset into a two-dimensional support funData object, then computing FPCA with MFPCA::MFPCA(). The twist here is that while this is a univariate model, we have to use MFPCA() and not PACE(). The reason is that PACE() does not accept multidimensional supports. The model is then a ‘multi’-variate one with only one response.
Functional data with 60 observations of 2-dimensional support
argvals:
0 0.2 ... 1 (6 sampling points)
1.545106 1.64118 ... 2.505843 (11 sampling points)
X:
array of size 60 x 6 x 11
# FCP_TPAfpca <-MFPCA(curvesFun %>%multiFunData(),M =5,# FCP_TPA is one of the possible multidim. support parametrizationsuniExpansions =list(list(type ="FCP_TPA",npc =5,alphaRange =list(v =c(10^-2, 10^2), w =c(10^-3, 10^3)))))# Prop of explained varPropExpVar <-round(fpca$values /sum( fpca$values) , digits =3)nPC <-2# scores st. dev.sdFun <- fpca$values %>%sqrt()# PC curves to be plottedPCcurves <-expand_grid(PC =1:nPC,fractionOfStDev =c(-1, 0, 1)) %>%group_by(PC, fractionOfStDev) %>%reframe( (fpca$meanFunction[[1]]@X[1,,] + fractionOfStDev * sdFun[PC] * fpca$functions[[1]]@X[PC,,]) %>%`colnames<-`(angle_grid) %>%as_tibble() %>%bind_cols(normTime = normTime_grid) ) %>%pivot_longer(cols =-c(PC, fractionOfStDev, normTime), names_to ="angle", values_to ="radius") %>%mutate(angle =as.numeric(angle))# PlotFPCA_angle_time_curves_plot <-ggplot(PCcurves %>%rename(`score/std`= fractionOfStDev)) +aes(x = angle, y = radius, group = normTime, color = normTime) +geom_path() +facet_grid(PC ~`score/std`, labeller =labeller(PC =~str_glue("PC{.x}"), `score/std`= label_both), ) +radial_plot(60) +theme(legend.position ="bottom")
Code
FPCA_angle_time_curves_plot
Figure 26: First 2 PCs: effect of scores on contour evolution in time
PC1 explains 98.4% of the variance, PC2 explains 0.7% of the variance. To avoid cluttering, the variability modulated by the scores is represented only by three values of standard deviation-normalized scores: -1, 0 and 1. For example, the top left quadrant shows the clip reconstructed by taking the mean clip and subtracting the first PC multiplied by one standard deviation of s1.
Clearly, PC1 is capturing the variability in the right part of the contour (tongue tip), which is the one correlated with Category. This is confirmed by the PC scores boxplots below.
radius ~ angle * time: FPCA predictions
# collect PC scoresPCscores <- fpca$scores[, 1:nPC] %>%`colnames<-`(paste0("s", 1:nPC)) %>%as_tibble() %>%bind_cols(curves %>%distinct(clip_id, Category), .)# boxplots PC scores by CategoryFPCA_angle_time_scores_plot <- PCscores %>%pivot_longer(cols =paste0("s", 1:nPC), names_to ="PC", values_to ="score") %>%filter(PC %in%str_c("s", 1:nPC)) %>%ggplot() +aes(x = Category, y = score, color = Category) +geom_boxplot() +facet_wrap(~ PC) +scale_color_manual(values=Category.colors) +theme(legend.position ="bottom")s <-1# score indexmodel_eq <-str_glue("s{s} ~ Category") %>%as.formula()mod <-lm(model_eq, data = PCscores)mod %>%summary()
Call:
lm(formula = model_eq, data = PCscores)
Residuals:
Min 1Q Median 3Q Max
-0.38541 -0.18281 -0.00389 0.16257 0.55832
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.13439 0.03898 80.41 <2e-16 ***
CategoryB -6.26877 0.05513 -113.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2135 on 58 degrees of freedom
Multiple R-squared: 0.9955, Adjusted R-squared: 0.9955
F-statistic: 1.293e+04 on 1 and 58 DF, p-value: < 2.2e-16
emmeans(mod, pairwise ~ Category)
$emmeans
Category emmean SE df lower.CL upper.CL
A 3.13 0.039 58 3.06 3.21
B -3.13 0.039 58 -3.21 -3.06
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B 6.27 0.0551 58 113.710 <.0001
Figure 27: Distribution of the first 2 PC scores by Category
A linear regression model predicting s1 from Category is then fitted, which provides the predicted clips below.
Code
FPCA_angle_time_pred_plot
Figure 28: Reconstructed clips by Category, based on regression model: s1 ~ Category
2.4 (x, y) ~ knot * time: GAM
Left as exercise. In short, any term that in Section 1.6 looks like s(knot, by = Category) should be turned into te(knot, normTime, by = Category).
2.5 (x, y) ~ knot * time: FPCA
This is the most complex model in this tutorial. Data objects have a two-dimensional support, (knot by normTime), the FPCA model is multivariate, as it takes values in x and y. The linear model following FPCA remains simple and geometry-agnostic.
(x, y) ~ knot * time: FPCA fit and predict
# build a funData objectnClips <- raw_curves$clip_id %>%n_distinct()argvals <-list(normTime=normTime_grid, knot=1:11)# curvesFun is a multiFunData with two response outcomes (x, y), # each of which is an array of dim c(nClips, length(normTime_grid), number of knots)curvesFun <-multiFunData( lapply(c('x', 'y'), function(Dim) {funData(argvals, raw_curves %>%pivot_wider(id_cols =c(clip_id, normTime), names_from = knot, values_from = {{Dim}}) %>%select(-normTime) %>%nest(.by = clip_id) %>%pull(data) %>%lapply(as.matrix) %>%abind(along =0) )}))curvesFun
An object of class "multiFunData"
[[1]]
Functional data with 60 observations of 2-dimensional support
argvals:
0 0.2 ... 1 (6 sampling points)
1 2 ... 11 (11 sampling points)
X:
array of size 60 x 6 x 11
[[2]]
Functional data with 60 observations of 2-dimensional support
argvals:
0 0.2 ... 1 (6 sampling points)
1 2 ... 11 (11 sampling points)
X:
array of size 60 x 6 x 11
# FCP_TPA# the same base expansion applied to x and y uniExpansion <-list(type ="FCP_TPA",npc =5,alphaRange =list(v =c(10^-2, 10^2), w =c(10^-3, 10^3)))fpca <-MFPCA(curvesFun,M =5,uniExpansions =list(uniExpansion, uniExpansion))# Prop of explained varPropExpVar <-round(fpca$values /sum( fpca$values) , digits =3)nPC <-2# scores st. dev.sdFun <- fpca$values %>%sqrt()# PC curves to be plottedPCcurves <-expand_grid(PC =1:nPC,fractionOfStDev =c(-1, 0, 1),Dim =1:2) %>%group_by(PC, fractionOfStDev, Dim) %>%reframe( (fpca$meanFunction[[Dim]]@X[1,,] + fractionOfStDev * sdFun[PC] * fpca$functions[[Dim]]@X[PC,,]) %>%`colnames<-`(1:11) %>%as_tibble() %>%bind_cols(normTime = normTime_grid) ) %>%pivot_longer(cols =-c(PC, fractionOfStDev, normTime, Dim), names_to ="knot", values_to ="Pos") %>%mutate(knot =as.integer(knot),Dim =factor(Dim, levels =1:2, labels =c('x','y'))) %>%pivot_wider(names_from = Dim, values_from = Pos)# PlotFPCA_knot_time_curves_plot <-ggplot(PCcurves %>%rename(`score/std`= fractionOfStDev)) +aes(x = x, y = y, group = normTime, color = normTime) +geom_path() +facet_grid(PC ~`score/std`, labeller =labeller(PC =~str_glue("PC{.x}"), `score/std`= label_both), ) +theme(legend.position ="bottom")# collect PC scoresPCscores <- fpca$scores[, 1:nPC] %>%`colnames<-`(paste0("s", 1:nPC)) %>%as_tibble() %>%bind_cols(raw_curves %>%distinct(clip_id, Category), .)# boxplots PC scores by CategoryFPCA_knot_time_scores_plot <- PCscores %>%pivot_longer(cols =paste0("s", 1:nPC), names_to ="PC", values_to ="score") %>%filter(PC %in%str_c("s", 1:nPC)) %>%ggplot() +aes(x = Category, y = score, color = Category) +geom_boxplot() +facet_wrap(~ PC) +scale_color_manual(values=Category.colors) +theme(legend.position ="bottom")s <-1# score indexmodel_eq <-str_glue("s{s} ~ Category") %>%as.formula()mod <-lm(model_eq, data = PCscores)mod %>%summary()
Call:
lm(formula = model_eq, data = PCscores)
Residuals:
Min 1Q Median 3Q Max
-1.14144 -0.35568 0.05615 0.34425 0.81853
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.42065 0.09094 -92.6 <2e-16 ***
CategoryB 16.84130 0.12860 131.0 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4981 on 58 degrees of freedom
Multiple R-squared: 0.9966, Adjusted R-squared: 0.9966
F-statistic: 1.715e+04 on 1 and 58 DF, p-value: < 2.2e-16
emmeans(mod, pairwise ~ Category)
$emmeans
Category emmean SE df lower.CL upper.CL
A -8.42 0.0909 58 -8.60 -8.24
B 8.42 0.0909 58 8.24 8.60
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B -16.8 0.129 58 -130.954 <.0001
Figure 29: First 2 PCs: effect of scores on contour evolution in time
PC1 explains 97% of the variance, PC2 explains 1.1% of the variance.
Code
FPCA_knot_time_scores_plot
Figure 30: Distribution of the first 2 PC scores by Category
A linear regression model predicting s1 from Category is then fitted, which provides the predicted clips below.
Code
FPCA_knot_time_pred_plot
Figure 31: Reconstructed clips by Category, based on regression model: s1 ~ Category
Citation information
For attribution, please cite this work as:
Michele Gubian (2025) Analysis of Uni- and Multi-Dimensional Contours with GAMs and Functional PCA: an Application to Ultrasound Tongue Imaging [Version 1.0.0]. https://uasolo.github.io/UTI-GAM-FPCA
References
Coretta, Stefano & Georges Sakr. 2025. Multivariate Analyses of Tongue Contours from Ultrasound Tongue Imaging. Draft V0.2. https://stefanocoretta.github.io/mv_uti/.
Gubian, Michele, Francisco Torreira & Lou Boves. 2015. Using functional data analysis for investigating multidimensional dynamic phonetic contrasts. Journal of Phonetics 49. Elsevier, 16–40.
Wieling, Martijn. 2018. Analyzing dynamic phonetic data using generalized additive mixed modeling: A tutorial focusing on articulatory differences between L1 and L2 speakers of english. Journal of Phonetics 70. Elsevier, 86–116.
Wieling, Martijn, Fabian Tomaschek, Denis Arnold, Mark Tiede, Franziska Bröker, Samuel Thiele, Simon N Wood & R Harald Baayen. 2016. Investigating dialectal differences using articulography. Journal of Phonetics 59. Elsevier, 122–143.
Wrench, Alan & Jonathan Balch-Tomes. 2022. Beyond the edge: Markerless pose estimation of speech articulators from ultrasound and camera images using DeepLabCut. Sensors 22(3). MDPI, 1133.
