Where we are

So far, every model in this course has been linear in its predictors. Many linguistic data are curves: formant trajectories, pitch contours, tongue contours, pupil size over a trial, looks to a target over time.

Today

1. Curves
2. Smooth terms
3. Fitting GAMs with {mgcv}
4. Comparing curves

The data

price

# A tibble: 2,200 × 6
  speaker l1      traj  duration  time    f2
  <fct>   <fct>   <fct>    <dbl> <dbl> <dbl>
1 S01     English S01.1      252     0  1199
2 S01     English S01.1      252    10  1232
3 S01     English S01.1      252    20  1174
4 S01     English S01.1      252    30  1194
# ℹ 2,196 more rows

The trajectories

ggplot(price, aes(x = time, y = f2, group = traj)) +
  geom_line(colour = kingfisher, alpha = 0.25) +
  facet_wrap(~ l1) +
  labs(x = "Time (% of vowel duration)", y = "F2 (Hz)")

Each line is one token. The research question is whether the two groups produce different F2 trajectories.

A straight line

The line misses the shape, and the residuals follow a pattern: neighbouring residuals are similar. This is autocorrelation, and it makes standard errors unreliable.

Polynomials

• A quadratic or cubic term underfits: it cannot follow the flat start and the late rise
• A high-degree polynomial overfits: it bends to follow the noise

Nothing in the data says which degree to choose.

Smooth terms

A generalised additive model can contain smooth terms next to the usual parametric terms:

f(x₂) is a smooth function of x₂ whose shape is estimated from the data. It is built as a sum of simpler functions, the basis functions, much as a polynomial is a sum of powers of x. A smoothing penalty keeps it from following the noise.

Basis functions

Each basis function is multiplied by a coefficient; the sum (salmon), plus the intercept, is the fitted curve. bs = "cr" asks for cubic regression splines, and k = 10 gives 10 basis functions, one of which is absorbed by the intercept.

How wiggly?

• k sets the maximum wiggliness: how many basis functions the smooth may use
• The smoothing parameter sets the actual wiggliness: the larger it is, the more wiggliness is penalised
• {mgcv} estimates the smoothing parameter from the data, so the strategy is to set k high enough and let the penalty do the rest
• The effective degrees of freedom (edf) describe the result: an edf of 1 is a straight line, and higher values mean more wiggly curves

The smoothing parameter

The same 10 basis functions with four smoothing parameters. A large penalty pushes the smooth toward a straight line.

gam() and bam()

• {mgcv} (Wood, 2017) fits GAMs with gam(), and with bam(), which is much faster on big data sets such as curves
• s(time) is a smooth over time; s(time, by = l1) fits a separate smooth for each group
• Smooths in a by model are centred on zero, so the group difference in height needs its own parametric term, l1
• method = "ML" is needed to compare models that differ in their fixed effects; the default of bam() is "fREML"

A first GAM

m1 <- bam(f2 ~ l1 + s(time, by = l1, k = 10), data = price, method = "ML")
round(summary(m1)$p.table, 2)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  1555.40       3.38  460.62     0.00
l1Italian      -4.78       4.78   -1.00     0.32

round(summary(m1)$s.table, 2)

                   edf Ref.df      F p-value
s(time):l1English 7.54   8.48 888.25       0
s(time):l1Italian 4.92   6.01 424.13       0

Reading the summary

• The parametric table works as in lm(): the intercept is the mean F2 of the English trajectories, and l1Italian is how much higher or lower the Italian trajectories are on average
• The smooth table has one row per smooth: edf close to 9 means a very wiggly curve, close to 1 a straight line

The F test asks whether a smooth is flat, so two identical curves would both come out significant. Testing whether two smooths differ needs a different model, later in this lecture.

Plotting with {tidygam}

m1_pred <- predict_gam(m1, series = "time", length_out = 50)
plot(m1_pred, series = "time", comparison = "l1")

predict_gam() returns a tibble of predictions with standard errors and 95% confidence intervals; its plot() method draws them with ggplot2.

Is k high enough?

k.check(m1)

                  k'      edf   k-index p-value
s(time):l1English  9 7.540370 0.5732537       0
s(time):l1Italian  9 4.924248 0.5732537       0

• k' is the number of basis functions available; edf is how many the smooth used
• A k-index below 1 with a small p-value suggests there is pattern left in the residuals that more basis functions could capture

To be sure, raise k (for example to 15) and refit. If the curves and the edf barely change, k was high enough.

Choosing k

• Start from what you expect: a formant trajectory over 11 points rarely needs more than 10 basis functions; a pitch contour over a whole sentence may need 20 or more
• Look at the fitted curves against the data, compare edf with k', and use k.check()

Reference and difference smooths

To test whether two curves differ, the model needs a term for the difference. With an ordered factor in by, s(time) becomes the reference smooth for the first level, and s(time, by = l1_o) the difference between the other level and the reference.

price <- mutate(price, l1_o = as.ordered(l1))
contrasts(price$l1_o) <- "contr.treatment"
m2 <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10),
          data = price, method = "ML")

The difference smooth

round(summary(m2)$p.table, 2)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  1555.40       3.38  460.36     0.00
l1_oItalian    -4.78       4.78   -1.00     0.32

round(summary(m2)$s.table, 2)

                     edf Ref.df      F p-value
s(time)             7.75   8.51 883.90       0
s(time):l1_oItalian 6.59   7.65 108.44       0

• s(time) is the English trajectory, and s(time):l1_oItalian is the difference smooth: the Italian trajectory minus the English one, apart from the overall difference in height captured by l1_oItalian
• m1 and m2 fit the same curves in two parametrisations

The binary smooth

price <- mutate(price, is_italian = as.numeric(l1 == "Italian"))
m2_bin <- bam(f2 ~ s(time, k = 10) + s(time, by = is_italian, k = 10),
              data = price, method = "ML")

• In a binary smooth, the by variable is numeric (0 or 1), so the smooth is not centred, so a single term covers both the difference in height and the difference in shape
• Its p-value tests both differences at once, which is often the question of interest

Which test to report?

Sóskuthy’s (2018) recommendations:

• Use the model summary only as a rough guide
• Report the significance from a model comparison
• Use prediction plots and the difference smooth for the qualitative discussion: where and how the curves differ

Summary

Slides and code

These slides

mcanzi.de/ss27_advanced_statistics/03_gams.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics