R/mod_03_trends_fct_ts_decomp.R

In fBedecarrats/bikeMonitoR: Suivi des comptages vélo à Nantes

# Cette fonction écarte du calcul de tendance les variables dès lors 
# qu'on a trois valeurs nulles (0) consécutives. Le 0 est alors remplacé
# par une valeur manquante. Ceci afin d'éviter de biaiser l'estimation de 
# tendance lors de pannes prolongées. Les valeurs écartées sont bien sûr
# conservées dans les données brutes.
anom_3nuls <- function(x) {
  x %>%
    arrange(ds) %>%
    mutate(y_orig  = y,
           nul_3_consec = ifelse((y == 0 & lead(y) == 0 & lag(y) == 0) |
                                   (y == 0 & lead(y) == 0 & lead(y, 2) == 0) |
                                   (y == 0 & lag(y) == 0 & lag(y, 2) == 0),
                                 1, 0),
           y = ifelse(nul_3_consec == 1, NA, y))
}

# Cette fonction applique l'algorithme prophet aux données de comptage
mouline_prophet <- function(x, 
                            ylim_min = TRUE, 
                            ecarte_3nuls = TRUE,
                            negatives_to_0 = TRUE,
                            jours_feries = prep_feries(),
                            titre_graph = "Variations récurrentes et détection des anomalies") {
  
  if (ecarte_3nuls == TRUE) {
    x <- anom_3nuls(x)
  }
  
  m <- x %>%
    prophet(holidays = jours_feries, yearly.seasonality = TRUE,
            interval.width = 0.90, 
            changepoints =  c('2015-01-01', '2016-01-01', '2017-01-01',
                              '2018-01-01', '2019-01-01', '2020-01-01',
                              '2020-03-17', '2020-05-11', '2020-01-01'),
            # n.changepoints = 50,
            changepoint.prior.scale = 0.5,
            changepoint.range = 1)
  
  future <- make_future_dataframe(m, periods = 365) %>%
    filter(ds >= min(x$ds[!is.na(x$y)], na.rm = TRUE), 
           ds <= max(x$ds[!is.na(x$y)], na.rm = TRUE))
  
  forecast <- predict(m, future) %>%
    filter(ds >= min(x$ds, na.rm = TRUE), ds <= max(x$ds, na.rm = TRUE)) %>%
    mutate(ds2 = as_date(ds)) %>% # prophet change le format en POSIXct
    left_join(x, by = c("ds2" = "ds")) %>%
    mutate(anom = ifelse(y > yhat_upper, 1, ifelse(y < yhat_lower, 1, NA)),
           anom_value = ifelse(is.na(anom), NA, y),
           anom_magnitude = ifelse(anom == 1, (y - yhat_upper)/y,
                                   ifelse(anom == -1, (yhat_lower-y)/y, NA)),
           anom = abs(anom),
           anom_magnitude = abs(anom_magnitude),
           `Probabilité\nd'anomalie\n(indice)` = log(anom_magnitude+1))
  
  # dans certains cas, le modèle cumulatif de forecast peut produire
  # des estimations négatives. Il est recommandé de les passer à 0 quand
  # la variable estimée ne peut pas être négative (comme c'est le cas ici)
  # cf. https://github.com/facebook/prophet, issues #1468, #1454, #1214...
  if (negatives_to_0) {
    forecast$yhat = ifelse(forecast$yhat < 0, 0, forecast$yhat)
  }
  
  
  # On enlève le y de forecast
  forecast$y <- NULL
  # On prépare le graph (un peu complexe, mais c'est plus joli comme ça)
  graph <- m$history %>%
    select(ds, y) %>%
    right_join(forecast, by = "ds") %>%
    arrange(ds) 
  # Si la valeur a été remplacée par NA car 3 valeurs nulles (O) consécutifs,
  # on réintègre la valuer initiale, on marque la valeur comme anormale mais
  # on n'inclut pas de "halo" proportionnel à la proba d'erreur.
  if (ecarte_3nuls) {
    graph <- graph %>%
      mutate(y = y_orig,
             anom = ifelse(nul_3_consec == 1, 1, anom),
             anom_value = ifelse(nul_3_consec == 1, y_orig, anom_value))
  }
  graph <- graph %>%
    ggplot(aes(x = ymd(paste("2000", month(ds), day(ds), sep = "-")), 
               y = y, shape = "Mesure")) +
    labs(y = "nombre de passages à vélo") +
    geom_point(size = 0.3, na.rm=TRUE) +
    scale_shape_manual("", values = 19) +
    geom_line(aes(y = yhat, color = "Modèle"),
              size = 0.2, na.rm = TRUE) +
    geom_ribbon(aes(ymin = yhat_lower, ymax = yhat_upper,
                    fill = "Marge\nd'erreur"),
                alpha = 0.5,
                na.rm = TRUE) +
    scale_fill_manual("", values = "#0072B2") +
    scale_colour_manual("",values="blue") +
    facet_grid(year(ds)~.) +
    geom_point(aes(y = anom_value, alpha = "Anomalie\npossible"),
               colour = "red", size = 0.3, shape = 16) +
    scale_alpha_manual("",values = 1) +
    geom_point(aes(size = `Probabilité\nd'anomalie\n(indice)`),
               colour = "red", alpha = 0.2) +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          legend.position = "bottom",
          legend.title = element_text(size = 9),
          plot.title = element_text(size = 11, face = "bold")) + #, suppr ')+
    # legend.position = "none") +
    scale_x_date(date_breaks = "1 month", date_labels = "%b") +
    guides(shape = guide_legend(order = 1),
           color = guide_legend(order = 2),
           fill = guide_legend(order = 3),
           alpha = guide_legend(order = 4,
                                override.aes=list(size = 3)),
           size = guide_legend(nrow = 2, order = 5),
           override.aes=list(size = 3))  +
    ggtitle(titre_graph)
  
  if (ylim_min) {
    graph <- graph + ylim(0, NA)
  }
  output <- list(m, forecast, graph)
  names(output) <- c("model", "forecast", "graph")
  return(output)
}

# Une fonction pour calculer la croissance annuelle
cr_an <- function(pr) {
  m_first <- pr[["forecast"]]$trend[1]
  m_last <- pr[["forecast"]]$trend[nrow(pr[["forecast"]])]
  d_first <- pr[["forecast"]]$ds[1]
  d_last <- pr[["forecast"]]$ds[nrow(pr[["forecast"]])]
  duree <- time_length(d_last - d_first, "year")
  c_tot <- round((m_last - m_first)/m_first * 100)
  c_an <- round(((m_last/m_first)^(1/duree)-1)*100,1)
  return(c_an)
}

# Une fonction pour appliquer le modèle aux diférents compteurs
process_boucle <- function(base, nom, ylim_min = FALSE) {
  # Filtre
  db <- base %>%
    filter(name == nom) %>%
    # retirer el dernier jour (décompte partiel)
    # filter(date != max(date)) %>%
    rename(ds = date, y = count)
  # Mouline en pr
  pr <- mouline_prophet(db, ylim_min, 
                        titre_graph = paste0(
                          "Variations récurrentes et détection des anomalies (",
                          nom, ")"))
  decomp <- prophet_plot_components(pr[["model"]], pr[["forecast"]], render_plot = FALSE)
  croiss_an <- cr_an(pr)
  output <- list(pr, decomp, croiss_an)
  names(output) <- c("modele", "composantes", "stat")
  return(output)
}

merge_current_historic <- function(current, historic) {
  current <- current %>%
    select(-weekday) %>%
    mutate(id = as.numeric(id),
           name = recode(name,
                         "Bonduelle vers sud" = "Pont A. Briand vers Sud",
                         "Bonduelle vers Nord" = "Pont A. Briand vers Nord",
                         "pont Anne de Bretagne vers Sud" = "Pont Anne de Bretagne vers Sud",
                         "pont Anne de Bretagne vers Nord" = "Pont Anne de Bretagne vers Nord",
                         "50 Otages Vers Sud" = "Cours des 50 Otages Vers Sud",
                         # Madeleine était cumulé en 2014-2019, on fait pareil ici
                         "Madeleine vers Sud" = "Chaussée de la Madeleine",
                         "Madeleine vers Nord" = "Chaussée de la Madeleine"),
           # idem pour les identifiant de compteur de la Madeleine
           id = ifelse(id == 880, 881, id)) %>%
    pivot_longer(-c("id", "date", "name",  "anom"),
                 names_to = "hour",
                 values_to = "count") %>%
    group_by(id, date, name, anom) %>%
    summarise(count = sum(count, na.rm = TRUE)) %>%
    mutate(anom = ifelse(anom == "Fonctionnel", 0, 1),
           id = as.numeric(id)) %>%
    filter(id %in% historic$id)
  
  
  current_names <- current %>%
    ungroup() %>%
    select(id, name) %>%
    unique()
  
  historic <- historic %>%
    select(-name) %>%
    left_join(current_names, by = "id")
  
  bind_rows(current, historic) %>%
    filter(!is.na(name))
  
}


loop_counters <- function(x) {
  # On récupère les noms de chaque compteur
  counters <- unique(x$name)
  # On enlève Tabarly car problème de date
  counters <- counters[counters != "Pont Tabarly vers Sud"]
  counters <- counters[counters != "Bd Malakoff vers Ouest"]
  counters <- counters[counters != "Bd Malakoff vers Est"]
  # On crée une liste vide
  out <- vector(mode = "list", length = length(counters))
  # On prépare les données pour chaque compteur
  for (i in 1:length(out)) {
    print(counters[i])
    out[[i]] <- process_boucle(x, counters[i])
  }
  names(out) <- counters
  return(out)
}


update_trends <- function() {
  dwn_bikecount_files(current = TRUE, historic = TRUE)
  trends <- prep_current() %>%
    merge_current_historic(prep_historic()) %>%
    loop_counters()
  save(trends, file = "trends.Rdata")
  return(trends)
}