In exercise testing (spiroergometry), participants' oxygen uptake is measured continuously while they engange in physical exercise that progressively intensifies. The measurements are taken at regular intervals, say every 10 seconds. It is useful to know whether the oxygen uptake reaches a plateau at the end of the measurements or not. We're talking about the last minute or so of the testing. Basically, it comes down to whether the oxygen uptake is still increasing at the end or if it is flattening. The following graph shows the measurement of three participants. In the left panel, no plateau is visible. The middle panel shows a clear plateau and the right panel is somewhat ambiguous.


There exist multiple heuristics of how to detect a plateau. All of them are somewhat arbitrary as there exists no ground truth or gold standard. Nevertheless, I'd like to know if there is a principled way to detect plateaus in time series.

My own idea was the following:

  1. Smooth the measurements using LOESS (black curves in the figure above). As the picture shows, the measurements are somewhat noisy. The amount of smoothing was chosen using a procedure explained here.
  2. If the standard deviation of the last $x$ (e.g. 6) measurements of the predicted smoothed data is below a certain threshold, a plateau is detected.

My thinking was that a small standard deviation signifies measurements that change little/are roughly constant. This still begs the question of how to chose the threshold for the standard deviation of the last smoothed measurements below which a plateau is declared.

My questions: Does the algorithm above make sense and how could it be improved? What other algorithms/heuristics could be used to detect plateaus at the end of the measurements?

Things I tried but didn't work well:

  • Calculate the confidence interval for the first derivative of a smooth function (e.g. calculated with packages mgcv and gratia). The confidence limits are just too large and almost always include $0$.
  • Some kind of changepoint detection. Again, the uncertainty is simply too large and the trajectory often can't be approximated well by linear functions.
  • 1
    $\begingroup$ I’m not certain but I think CUSUM can be used for this purpose. $\endgroup$
    – Sycorax
    Commented May 4 at 16:07
  • $\begingroup$ @Sycorax That's a great suggestion which I had as well. I got lost in the details of how to set the detection limits and which measurements to take etc. So I'd be grateful for any details and further thoughts. $\endgroup$ Commented May 4 at 16:13
  • 2
    $\begingroup$ I have to say that I read your examples differently. To me "clear plateau" is ambiguous and ambiguous is "no plateau". Dependence on some smoothing method may be inescapable but that's a long slippery slope in which everything hinges on which method, how much to smooth and what is to be done with the results. $\endgroup$
    – Nick Cox
    Commented May 5 at 15:20
  • 2
    $\begingroup$ If your interest is in whether a maximum oxygen consumption rate has been achieved, I'm not sure that the answer lies solely in the statistical analysis. See Sports Med. 2021; 51: 1815 for a recent review of methods used to estimate such plateaus and the physiological issues that can further complicate the interpretation of a plateau (or lack of one) in terms of an individual's maximum oxygen consumption. Also, Front Physiol. 2013; 4: 203. $\endgroup$
    – EdM
    Commented May 5 at 20:09

1 Answer 1


I would use a Generalized Additive Model (GAM) to tackle this question. The penalized splines fitted with GAMs will be less likely to overfit than your LOESS algorithm, and you can use hierarchical effects to help learn these patterns from multiple time series in a single joint model (see this nice tutorial paper by Pedersen et al for guidance). Below I use some simulated data to show how you can fit such a model using the {mgcv} package. I then show how you can use the powerful {marginaleffects} package to calculate the 1st derivative (slope) at each point along a fine grid of covariate values. This can be useful to test the hypothesis that the slope is different from zero, which can provide a quantitatively more rigorous way to determine whether a plateau has been reached.


# Simulate a nonlinear curve that plateaus
x <- sort(runif(80) * 4 - 1)
f <- exp(4 * x) / (1 + exp(4 * x))
y <- f + rnorm(80) * 0.1
dat <- data.frame(y, x)

# Plot the simulated points
ggplot(dat, aes(x, y)) +

# Model the curve with a GAM
fit <- gam(y ~ 
             # A b-spline is useful for estimating derivatives near
             # the boundaries
             s(x, bs = "bs", k = 10),
           # Extending the penalty beyond the boundaries can help 
           # estimate derivatives near the bounds
           knots = list(x = c(-2, 0, 3, 4)),
           data = dat)
#> Warning in smooth.construct.bs.smooth.spec(object, dk$data, dk$knots): there is
#> *no* information about some basis coefficients

# View the estimated function, on the response scale
plot_predictions(fit, condition = 'x',
                 points = 0.5)

# View th etimated slope (1st derivative) 
# (using a finer grid of points for a smoother image)
            newdata = datagrid(x = seq(-1, 3, by = 0.1)),
            variables = 'x', by = 'x') +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(y = '1st derivative (slope)',
       x = 'x')

# Use a hypothesis test to determine whether the slope is 
# significantly different from zero at each point on the grid
                  newdata = datagrid(x = seq(-1, 3, by = 0.1)),
                  variables = 'x', by = 'x'))
#>  Term    Contrast    x Estimate Std. Error      z Pr(>|z|)    S  2.5 %  97.5 %
#>     x mean(dY/dX) -1.0  -1.1386     0.4727 -2.409   0.0160  6.0 -2.065 -0.2122
#>     x mean(dY/dX) -0.9  -0.6795     0.3402 -1.997   0.0458  4.4 -1.346 -0.0127
#>     x mean(dY/dX) -0.8  -0.2777     0.2299 -1.208   0.2271  2.1 -0.728  0.1729
#>     x mean(dY/dX) -0.7   0.0669     0.1449  0.462   0.6440  0.6 -0.217  0.3508
#>     x mean(dY/dX) -0.6   0.3544     0.0934  3.793   <0.001 12.7  0.171  0.5375
#> --- 31 rows omitted. See ?print.marginaleffects --- 
#>     x mean(dY/dX)  2.6   0.0231     0.1578  0.147   0.8835 0.2 -0.286  0.3323
#>     x mean(dY/dX)  2.7   0.0614     0.1510  0.407   0.6843 0.5 -0.235  0.3574
#>     x mean(dY/dX)  2.8   0.1061     0.2481  0.428   0.6687 0.6 -0.380  0.5923
#>     x mean(dY/dX)  2.9   0.1573     0.4632  0.340   0.7341 0.4 -0.751  1.0653
#>     x mean(dY/dX)  3.0   0.2150     0.7660  0.281   0.7790 0.4 -1.286  1.7164
#> Columns: rowid, term, contrast, x, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

See the help pages from {marginaleffects} for more guidance on the types of hypotheses you can test. If you find that you need some more flexibility in your smooth functions, for example if you need to ensure the function is monotonically increasing (which sounds sensible for your data), you can see my package {mvgam}. This was designed for fitting Dynamic GAMs to time series and can provide you more options (see for example the help page on monotonic smooths)

Created on 2024-05-06 with reprex v2.1.0


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.