variation: Analyse Coefficient Variation Systématique

#' Calcul des CVS et des gains
#'
#' Tout ce qui est relatif aux CVS ainsi que les gains selon différentes hypothèses.
#'
#' \strong{pop.region} : Tableau ayant les colonnes suivantes : an, geocode, pop, pop.totale.
#'
#' @param x BD.
#' @param col.an CHR. Nom de la colonne indiquant l'année ou la période.
#' @param col.region CHR. Nom de la colonne indiquant la région.
#' @param col.code CHR. Nom de la colonne indiquant les codes de regroupement.
#' @param col.pop CHR. Nom de la colonne indiquant la population totale.
#' @param col.obs CHR. Nom de la colonne indiquant le nombre d'observations.
#' @param col.sexe CHR. Facultatif. Nom de la colonne indiquant le sexe.
#' @param col.catage CHR. Facultatif. Nom de la colonne indiquant les catégories d'âge.
#' @param col.indic CHR. Facultatif. Nom de la colonne indiquant l'indicateur social.
#' @param cvs INT. Valeur(s) du(des) cvs.
#' @param nsim INT. Nombre de simulations.
#' @param p100 NUM. Pourcentages cibles. Ex : 0.40.
#' @param pctl NUM. Percentiles cibles. Ex : 0.40.
#' @param alpha0 INT. Itérateur initial (par défaut 1).
#' @param min.zone,min.obs INT. Nombre minimum de zones géographiques (\code{min.zone}) ayant au moins \code{min.obs} observations.
#' @param pop.region DATA. Tableau des populations. Voir détails.
#'
#' @import plyr
#' @import data.table
#'
#' @export
cvs_val <- function (x,
                     col.an = "an",
                     col.region = "region",
                     col.code = "code",
                     col.pop = "pop",
                     col.obs = "obs",
                     col.sexe = NULL,
                     col.catage = NULL,
                     col.indic = NULL,
                     cvs = c(3,2,1),
                     nsim = 20000,
                     p100 = c(0.4, 0.3, 0.2, 0.1),
                     pctl = c(0.4, 0.3, 0.2),
                     alpha0 = 1,
                     min.zone = 10,
                     min.obs = 100,
                     nyear = 1) {

  dt <- as.data.table(x$BD.analyse)
  # Liste des colonnes nécessaires
  dtnames <- c(col.an, col.code, col.region, col.obs, col.pop)
  # Ajouter colonnes facultatives si non NULL
  if (!is.null(col.sexe)) dtnames <- c(dtnames, col.sexe)
  if (!is.null(col.catage)) dtnames <- c(dtnames, col.catage)
  if (!is.null(col.indic)) dtnames <- c(dtnames, col.indic)
  # Sélection des colonnes nécessaires
  dt <- dt[, ..dtnames]
  # Renommer les colonnes
    # Nécessaires
  setnames(dt,
           c(col.an, col.region, col.code, col.pop, col.obs),
           c("an", "region", "code", "pop", "obs"))
    # Facultatives
  if (!is.null(col.sexe)) setnames(dt, col.sexe, "sexe")
  if (!is.null(col.catage)) setnames(dt, col.catage, "catage")
  if (!is.null(col.indic)) setnames(dt, col.indic, "indic")
  # Enregistrer les nouveaux noms
  dtnewnames <- copy(names(dt))
  dt[, pop := NULL]

  #### Calcul des expected ####
  # Créer toutes les combinaisons possibles
  possibilities <- list(
    an = unique(x$BD.analyse$an),
    code = unique(x$BD.analyse$code),
    region = unique(x$pop.region$geocode)
  )
  if(!is.null(col.sexe)) possibilities[["sexe"]] <- unique(dt$sexe)
  if(!is.null(col.catage)) possibilities[["catage"]] <- unique(dt$catage)
  if(!is.null(col.indic)) possibilities[["indic"]] <- unique(dt$indic)
  exp <- as.data.table(expand.grid(possibilities))
  dt <- merge(dt, exp, all = T)
  dt <- dt[is.na(obs), obs := 0]

  # Calcul des ratios par strate
  by.x <- c("an", "code", "obs")
  if(!is.null(col.sexe)) by.x <- c(by.x, "sexe")
  if(!is.null(col.catage)) by.x <- c(by.x, "catage")
  if(!is.null(col.indic)) by.x <- c(by.x, "indic")
  ratio <- dt[, ..by.x]
  by.x <- by.x[-3]
  ratio <- ratio[, .(obs = sum(obs)), keyby = by.x]
  by.merge <- by.x[-2]
  ratio <- merge(ratio, x$pop.strates, by = by.merge, all.x = T)
  ratio[, ratio := obs / pop]
  ratio_col <- c(by.x, "ratio")
  dt <- merge(dt, ratio[, ..ratio_col], by = by.x)  # ajout des ratios à dt
  pop.region.strates <- copy(x$pop.region.strates)
  setnames(pop.region.strates, "geocode", "region")
  dt <- merge(dt, pop.region.strates, all.x = T, by = c(by.merge, "region"))  # ajout de la population par strate
  dt <- dt[is.na(pop), pop := 0]
  dt[, exp := ratio * pop]
  dt <- dt[
    , .(obsj = sum(obs),
        expj = sum(exp),
        pop = sum(pop)),
    .(an, code, region)
  ]
  dt[, pop.totale := sum(pop), .(an, code)]
  dt[, RTj := obsj / expj]
  dt[, Tbrute := sum(obsj) / sum(pop) / nyear, .(an, code)]


  # conserver les codes qui respectent les conditions : min.zone + min.obs
  keep <- dt[obsj >= min.obs, .N, .(an, code)][N > min.zone, .(an, code)]
  dt <- merge(dt, keep,
              by = c("an", "code"))
  # Créer tableau des codes non conservées
  exclu <- as.data.table(x$BD.analyse)
  exclunames <- c(col.an, col.code)
  exclu <- unique(exclu[, ..exclunames])
  setnames(exclu, exclunames, c("an", "code"))
  exclu <- merge(dt, exclu, by = c("an","code"), all = T)
  exclu <- exclu[is.na(region), .(an, code)]


  #### Calculs ####
  dt <- as.data.table(rbind.fill(lapply(split(dt, by = c("an", "code")), function(z) {

    if(nrow(z)){
      #### CVS & cv ####
      # CVS
      nregion <- as.numeric(z[, .N])  # nombre de régions = nbr de lignes
      # Processus itératif
      alpha1 <- mean(z[["RTj"]]) / (sum((1 + alpha0 / z[["expj"]]) * (z[["RTj"]] - mean(z[["RTj"]]))^2) / (nregion - 1))
      nu <- alpha1 * mean(z[["RTj"]])
      for (i in 1:30) {
        set(z, NULL, "thetaj", (nu + z[["obsj"]]) / (alpha1 + z[["expj"]]))
        var <- (1 + (alpha1 / z[["expj"]])) * (z[["thetaj"]] - mean(z[["thetaj"]]))^2
        alpha1 <- mean(z[["thetaj"]]) / (sum(var) / (nregion - 1))
        nu <- mean(z[["thetaj"]]) * alpha1
      }



      # cv du CVS
      dt_sim <- copy(z)
      for (sim in seq(nsim)) {  # nombre de simulations
        obs_sim <- rnbinom(nregion, nu, mu = (dt_sim[["expj"]] * nu) / alpha1)  # valeurs observées simulées
        RTj_sim <- obs_sim / dt_sim[["expj"]]  # RTj simulés
        var <- (1 + (alpha0 / dt_sim[["expj"]])) * (RTj_sim - mean(RTj_sim))^2
        alpha_sim <- mean(RTj_sim) / (sum(var)/(nregion-1))
        nu_sim <- mean(RTj_sim) * alpha_sim
        dt_sim[["obsj"]] <- obs_sim  # insérer les obs simulés pour analyse
        for (i in 1:30) {
          set(dt_sim, NULL, "theta", (nu_sim + dt_sim[["obsj"]]) / (alpha_sim + dt_sim[["expj"]]))
          var <- (1 + (alpha_sim / dt_sim[["expj"]])) * (dt_sim[["theta"]] - mean(dt_sim[["theta"]]))^2
          alpha_sim <- mean(dt_sim[["theta"]]) / (sum(var)/(nregion-1))
          nu_sim <- mean(dt_sim[["theta"]]) * alpha_sim
        }
        # Conserver les CVS simulés
        if (sim == 1) {
          CVS_sim <- nu_sim / alpha_sim^2 * 100
        } else {
          CVS_sim <- c(CVS_sim, (nu_sim / alpha_sim^2 * 100))
        }
      }
      # Résultats
      z[, `:=` (THETA = expj * thetaj,  # nécessaire au gains % et perc
                CVS = nu / alpha1^2 * 100,  # CVS
                cv = sqrt(var(CVS_sim)) / mean(CVS_sim) * 100,  # cv du CVS
                QP9010 = as.numeric(quantile(z[["thetaj"]], 0.90) / quantile(z[["thetaj"]], 0.10)))]  # UR / QP9010


      #### Gains CVS hors normes ####
      z <- cbind(z, Reduce(function(a,b) cbind(a,b), lapply(cvs, function(cvsx) {  # Ajouter les colonnes créées
        gaincvs <- z[, .(obsj, expj, RTj, THETA, thetaj)]  # colonnes nécessaires
        cvs_ref <- cvsx / 100  # CVS de référence
        alpha_ref <- mean(gaincvs[["RTj"]]) / cvs_ref
        nu_ref <- mean(gaincvs[["RTj"]]) * alpha_ref
        for (i in 1:30) {
          set(gaincvs, NULL, "theta", (nu_ref + gaincvs[["obsj"]]) / (alpha_ref + gaincvs[["expj"]]) )
          var <- (1 + (alpha_ref / gaincvs[["expj"]])) * (gaincvs[["theta"]] - mean(gaincvs[["theta"]]))^2
          alpha_ref <- mean(gaincvs[["theta"]]) / cvs_ref
          nu_ref <- mean(gaincvs[["theta"]]) * alpha_ref
        }
        # Limites
        liminf <- qgamma(0.025, nu_ref, alpha_ref)  # limite inférieures
        limsup <- qgamma(0.975, nu_ref, alpha_ref)  # limite supérieures
        gaincvs[
          # Gain par région (par défaut 0)
          , `:=` (ginf_cvs = 0,
                  gsup_cvs = 0)
          ][
            # Gain si en-dessous limite inférieure
            thetaj < liminf,
            ginf_cvs := expj * liminf - THETA
            ][
              # Gain si au-dessus limite supérieure
              thetaj > limsup,
              gsup_cvs := expj * limsup - THETA
              ]
        # Tableau des résultats
        gaincvs <- gaincvs[
          , .(liminf_cvs = liminf,  # valeur de la limite inférieure selon cvs
              ginf_cvs = ginf_cvs / length(an.analyse),  # gain inférieure par région
              gaininf_cvs = sum(ginf_cvs),  # gain inférieure total
              limsup_cvs = limsup,  # valeur de la limite supérieure selon cvs
              gsup_cvs = gsup_cvs / length(an.analyse),  # gain supérieure par région
              gainsup_cvs = sum(gsup_cvs))  # gain supérieure total
          ]
        # Renommer les colonnes selon paramètres d'analyse
        setnames(gaincvs,
                 c("liminf_cvs", "limsup_cvs", "ginf_cvs", "gaininf_cvs", "gsup_cvs", "gainsup_cvs"),
                 c(paste0("liminf_cvs",cvsx),
                   paste0("limsup_cvs",cvsx),
                   paste0("ginf_cvs",cvsx),
                   paste0("gaininf_cvs",cvsx),
                   paste0("gsup_cvs",cvsx),
                   paste0("gainsup_cvs",cvsx)),
                 skip_absent = TRUE)
        gaincvs
      })))


      #### Pourcentage ####
      z <- cbind(z, Reduce(function(a,b) cbind(a,b),lapply(p100, function(p100x) {
        # Tableau avec les résultats
        gain100 <- z[, .(obsj)]  # colonnes nécessaires
        gain100[, obs_F := obsj * (1 + p100x)]  # obs avec modification selon pourcentage
        gain100[, g100 := (obs_F - obsj) / length(an.analyse)]  # gains par région
        gain100[, gain100 := sum(g100)]  # gains totaux
        # Renommer les colonnes selon paramètres d'analyse
        gain100 <- gain100[, .(g100, gain100)]
        setnames(gain100, c("g100", "gain100"),
                 c(paste0("gain_p",p100x*100), paste0("gaintot_p",p100x*100)))
        gain100
      })))


      #### Percentiles ####
      pctl <- sort(c(pctl, sapply(pctl, function(x) 1 - x)))  # complement du percentile : 1 - pctl
      z <- cbind(z, Reduce(function(a,b) cbind(a,b),lapply(pctl, function(pctlx) {
        gainpctl <- z[, .(obsj)]
        pctlx_F <- (quantile(z$obsj, pctlx) / mean(z$obsj)) - 1  # transformer le pctl en pourcentage
        gainpctl[, obs_F := obsj * (1 + pctlx_F)]  # obs avec modification selon pourcentage
        gainpctl[, gpctl := (obs_F - obsj) / length(an.analyse)]  # gains par région
        gainpctl[, gainpctl := sum(gpctl)]  # gains totaux
        gainpctl <- gainpctl[, .(gpctl, gainpctl)]
        # Renommer les colonnes selon paramètres d'analyse
        setnames(gainpctl, c("gpctl", "gainpctl"),
                 c(paste0("gpctl",pctlx*100), paste0("gainpctl",pctlx*100)))
        gainpctl
      })))


      #### Pourcentage + CVS ####
      z <- cbind(z, Reduce(function(a,b) cbind(a,b), mapply(function(p100x, cvsx) {
        # Réduction de % selon normalisation
        gainsup <- gain_pourcent_norm(z, cvsx, -p100x, method = 1)
        setnames(gainsup, c("g100", "gain100", "tj.region"),
                 c(
                   paste0("gsup_p",p100x*100,"_cvs",cvsx),
                   paste0("gainsup_p",p100x*100,"_cvs",cvsx),
                   paste0("Tj_sup_region_p",p100x*100,"_cvs",cvsx)
                 ))
        # Accroissement de % selon normalisation
        gaininf <- gain_pourcent_norm(z, cvsx, p100x, method = 1)
        setnames(gaininf, c("g100", "gain100", "tj.region"),
                 c(
                   paste0("ginf_p",p100x*100,"_cvs",cvsx),
                   paste0("gaininf_p",p100x*100,"_cvs",cvsx),
                   paste0("Tj_inf_region_p",p100x*100,"_cvs",cvsx)
                 ))
        # Concaténer résultats
        gain <- cbind(gainsup, gaininf)

        return(gain)
      },
      # Arguments
      p100,  # pourcentages
      rep(cvs, each = length(p100)),  # CVS
      SIMPLIFY = FALSE # Résultat est une liste
      )))

    } else {
      NULL
    }

  })))


  list(
    result = dt,
    exclus = exclu
  )

}