Grouping/Selected Layers using a threshold

Question

I have a raster dataset with 360 layers representing monthly data from January 1981 to December 2010. I want to calculate the drought events, defined as consecutive months with values less than or equal to -1. A drought event begins when the value stays below -1 for at least 2 consecutive months and ends when the value becomes greater than -1. The duration of a drought event, from its start (Dstart) to its end (Dend), is termed as the drought duration. I have already created a function to calculate the average drought duration per year. Below is the previous code:

drought_dur <- function(x, na.rm = TRUE) {
  drought_lengths <-  with(rle(x <= -1), lengths * values)
  print(drought_lengths[drought_lengths>=2])
  dd <- drought_lengths[drought_lengths >= 2]
  return(sum(dd))
}
SPI <- c(-1, -1, -1, -0.5, -1, 2, 1, -1, -1, -0.8, 0.5, 1, -1, -1, 1, -1,-1,0.8)
> drought_dur(SPI)
 [1] 3 2 2 2 
[1] 9

# So, we can see that it returns to the sum (total duration).

I've also calculate the raster data using the same function.

> model_test
class      : RasterBrick 
dimensions : 48, 48, 2304, 360  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : 94.125, 106.125, -5.875, 6.125  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : r_tmp_2024-03-18_180102.778979_85605_65317.grd 
names      :     layer.1,     layer.2,     layer.3,     layer.4,     layer.5,     layer.6,     layer.7,     layer.8,     layer.9,    layer.10,    layer.11,    layer.12,    layer.13,    layer.14,    layer.15, ... 
min values :          NA,          NA, -1.98604915, -1.62200621, -1.01830699, -0.56920946, -1.36586824, -2.63141197, -1.24827948, -0.88913461, -0.42232346, -1.25739376, -2.48757464, -2.82657113, -1.91978022, ... 
max values :          NA,          NA,  0.79760924,  0.93382244,  1.46647213,  1.25711339,  2.00315189,  1.17795110,  1.81855801,  1.83599187,  2.27933524,  1.30173316, -0.14318272,  0.68711265,  1.34114837, ..

And, we've got 30 layers of total duration for each year because it returns to the sum.

duration_yearly <-stackApply(model_test, rep(seq(1,nlayers(model_test)/12),
                                       each = 12), drought_dur)
>duration_yearly
class      : RasterBrick 
dimensions : 48, 48, 2304, 30  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : 94.125, 106.125, -5.875, 6.125  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : index_1, index_2, index_3, index_4, index_5, index_6, index_7, index_8, index_9, index_10, index_11, index_12, index_13, index_14, index_15, ... 
min values :      NA,       0,       0,       0,       0,       0,       0,       0,       0,        0,        0,        0,        0,        0,        0, ... 
max values :      NA,       9,       7,       0,       6,       3,       7,       2,       3,        6,        5,        4,        2,        6,        4, .

But, how do we calculate not the total duration for each year, but the duration for every drought events ? Means it returns to [1] 3 2 2 2 (example per year) not sum it.