2

I want to model the relationship between individuals' weekly mean speed and the timing in the year (represented by the week number), to then extract the dates where we observe a change in the animals' behaviour. I have several ID, accross multiple years (but IDs are not necessarily repeated every year).

The relationship looks like that (it is just an example and my dataset is much larger): enter image description here

I want to extract the breakpoints, i.e. the weeks (on the x-axis) where we can observe a change, here for example around weeks 18, 37 and 50.

I tried with the package segmented to extract the breakpoints, the results and residuals are pretty good but I can't account for any random effects or information about cyclic cubic spline or temporal autocorrelation, as segmented only deals with lm or glm object. However I need to model the relationship with those parameters.

I was wondering if it is possible to extract those breakpoints from a mgcv:gamm object? Like with the following model:

model <- gamm(Speed ~ s(Week, bs = "cc"), correlation = corAR1(form = ~ Week | ID/Year), random = list(ID = ~ 1, Year = ~ 1), data = dat)

I saw some threads talking about the model$gam$smooth$xp arguments in the model output, but those values do not seem to correspond to the breakpoints I see in the graphic output.

I also tried to apply a segmented model on the predicted values from the gamm, but the residuals do not really look good and I think it multiplies the methods which is maybe not a good idea.

Here is a reproducible example:

dat <- structure(list( 
  ID = c(4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,6,6,6,
         6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7),
  Year = c(2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2013,2013,2013,2013,2013,2013,2013,2013,2013,
           2013,2013,2013,2013,2013,2013,2013,2013,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,
           2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2012,2013,2013,2013,2013,2013,2013,2013,2013,
           2013,2013,2013,2013,2013,2013,2013,2013,2013,2012,2012),
  Speed = c(7.577964,9.670290,11.219523,14.657209,13.350411,14.340927,14.203197,14.334321,14.033588,14.761327,14.328223,14.031016,14.347175,14.515277,
               13.863412,14.026990,26.776336,32.723147,39.243111,39.934073,40.435426,39.799271,39.519325,38.298101,38.247134,37.781086,37.431256,38.415292,
               30.570140,14.686493,14.679004,24.997272,25.025281,11.598623,18.378834,15.933284,14.078759,14.197417,35.862854,32.265989,30.287343,18.903351,
               24.378491,22.726181,22.749912,20.870110,14.752451,13.344393,13.461780,30.726149,26.708165,29.985766,26.163577,22.178019,20.518806,21.453945,
               19.909988,19.130840,21.175718,22.096469,22.131889,21.793827,15.320547,21.300952,21.358867,22.028204,21.447345,21.283998,21.704764,21.742116,
               22.099805,21.872716,22.239127,21.422524,22.006441,22.425995,31.680002,35.859048,38.537170,36.753729,39.346547,39.337961,41.574394,41.378615,
               41.417221,41.175806,41.540436,41.515605,21.708493,30.846135,34.337545,30.533737,28.536433,22.454174,24.001463,30.824520,29.026101,22.555618,
               21.763984,28.320257,30.267227),
  Week = c(24,26,27,28,29,30,31,32,33,36,37,38,39,40,41,42,43,45,46,47,48,49,51,52,53,1,2,3,4,5,6,8,9,11,13,14,16,17,7,8,9,10,11,12,13,14,15,16,17,8,
            9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,25,26,27,28,30,31,33,35,36,40,41,42,43,44,45,46,47,49,50,51,52,53,1,5,6,7,8,9,10,11,12,13,14,17,8,9)))

Thanks a lot for your help!

airC
  • 21
  • 2

1 Answers1

1

One way to do this is via the derivatives of the estimated smooth, and look where the uncertainty band around the derivative doesn't include 0 to identify periods of "change" and then breakpoints are where the derivative shifts from indistinguishable from 0 to distinguishable from 0 or vice versa.

Here's an example using derivatives() from my {gratia}

Fit a model so we have something to work with

library('gratia')
library('dplyr')
library('ggplot2')
load_mgcv()
set.seed(42)
op <- options(cli.unicode = FALSE)
dat <- gamSim(1, n = 400, dist = "normal", scale = 2, verbose = FALSE)
mod <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, method = "REML")

Compute derivatives for a selected term

fx2_d <- derivatives(mod, type = "central")

Use the credible interval on the derivative to decide where the function is changing:

fx2_d <- fx2_d %>%
    mutate(change = if_else(lower < 0 & upper > 0, NA_real_, derivative))

Now change contains the value of the derivative where there is change and NA when 0 is within the credible interval. Plotting all this gives:

draw(fx2_d) +
  geom_line(data = fx2_d, aes(x = data, y = change), lwd = 2)

enter image description here

where the thick lines indicate the periods of change.

The finer the grid you evaluate the derivative on the more precise you can make the estimate of the change point, up to a point.

I'll leave extracting the break points to you as best suits your application. One way to do it would be to as if change is NA, convert to that logical to a numeric and then diff() that numeric vector

with(fx2_d, diff(is.na(change) * 1))

which gives

> with(fx2_d, diff(is.na(change) * 1))
  [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
 [51]  0 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 [76]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0
[101]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
[126]  0  0  0  0  0  0 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
[151]  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0
[176]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

where we can see for example that at value 46 of x2 the derivative flipped from increasing to including 0, and then a few values later (x2[52]) the derivative flips negative, indicating that this changepoing was at a peak of function. The function looks like this, by the way:

enter image description here

And you can use these indices to find the value of x2 at the change points:

> fx2_d[c(46:47, 52:53), ]
# A tibble: 4 x 9
  smooth var    data derivative    se  crit  lower upper change
  <chr>  <chr> <dbl>      <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>
1 s(x2)  x2    0.229      10.1   4.09  1.96   2.07 18.1    10.1
2 s(x2)  x2    0.234       7.01  4.16  1.96  -1.15 15.2    NA  
3 s(x2)  x2    0.259      -7.91  4.67  1.96 -17.1   1.24   NA  
4 s(x2)  x2    0.264     -10.6   4.74  1.96 -19.9  -1.36  -10.6
Gavin Simpson
  • 37,567
  • 5
  • 110
  • 153