R/polycumulpond.R

In INESSSQC/polymedic: Calculs de la polypharmacie

#' Polypharmacie Cumulée Pondérée
#'
#' Pour chaque médicament consommé par une personne, on donne une valeur entre 0 et 1. Cette valeur représente le ratio entre le nombre de jours avec consommation et nombre de jours total. Le résultat final pour chaque personne est l'addition des ratios.
#'
#' @param x Base de données créée par \code{polyBD()}.
#' @param t1,t2 "AAAA-MM-JJ". Intervalle de temps [t1, t2] servant d'analyse.
#' @param catage Borne(s) inférieure(s) des groupes d'âge.
#' @param classes Classes pour la distribution des maximums. Par défaut 5. Avec \code{classes = 5}, les classes seront [0,5[; [5,10[; [10,15[, etc.
#' @param maxclasses Par défaut NULL, sinon INTEGER. Valeur du minimum de la dernière classe à afficher. Regroupe toutes les valeurs plus grandes ou égales à \code{maxclasses}.
#'
#' @import data.table
#'
#' @export
polycumulpond <- function(x, t1, t2,
                          catage = c(60, 70, 80, 90),
                          classes = 5, maxclasses = NULL){

  ## Noms des colonnes
  # Conserver
  xcols_ <- names(x)
  # Modifier noms
  names(x)[[1]] <- "BenBanls"
  names(x)[[2]] <- "DIN"
  # Liste des ID
  xID <- unique(x[, BenBanls])
  # Convertir t1 et t2 en Date
  t1 <- as.Date(t1)
  t2 <- as.Date(t2)
  # BD Age et Sexe
  BDage <- attributes(x)$Naissance
  BDsex <- attributes(x)$Sexe

  # Sélection des observations inclus dans l'interval [t1, t2]
  x <- x[PerDebut <= t2 & PerFinAjust >= t1]

  # Calculer le ratio de consommation sur l'interval [t1, t2]
  # Debut de l'interval :
  x[, Deb := PerDebut][Deb < t1, Deb := t1] #si PerDebut est avant l'interval d'analyse, on prend t1; sinon PerDebut
  # Fin de l'interval
  x[, Fin := PerFinAjust][Fin > t2, Fin := t2] #si PerFinAjust est après l'interval d'analyse, on prend t2; sinon PerFinAjust
  # Proportion consommation sur interval
  x[, frac := as.numeric((Fin - Deb + 1) / as.numeric(t2 - t1 + 1))]
  # Additionner les fractions pour un même médicament si plusieurs périodes différentes dans l'interval d'analyse
  x <- x[, .(nfrac = sum(frac)), .(BenBanls)]



  ## Ajouter les BenBanls manquants - sans consommation
  xID <- xID[!xID %in% x$BenBanls]
  # Tableau indiquant leur consommation = 0
  if(length(xID) != 0){
    xID <- data.table(BenBanls = xID,
                      nfrac = 0L)
    x <- rbind(x, xID)
  }
  # Fusion des tableaux
  setorder(x, BenBanls)


  ## Nombre d'usager
  nID <- nrow(x)

  ## Résultats
  # Statistiques
  cumulpond_stats <- x[, .(Moyenne = mean(nfrac),
                           `Écart-type` = sd(nfrac),
                           Min = min(nfrac),
                           P5 = quantile(nfrac, 0.05),
                           P10 = quantile(nfrac, 0.10),
                           Q1 = quantile(nfrac, 0.25),
                           `Médiane` = median(nfrac),
                           Q3 = quantile(nfrac, 0.75),
                           P90 = quantile(nfrac, 0.90),
                           P95 = quantile(nfrac, 0.95),
                           Max = max(nfrac),
                           n = nID)]
  # BD d'analyse
  cumulpond_BD <- copy(x)
  names(cumulpond_BD)[[1]] <- xcols_[[1]]

  # Ajout de la date de naissance
  x <- merge(x, unique(BDage),
             by.x = "BenBanls", by.y = xcols_[[1]])
  # Convertir date de naissance en Age (selon la date d'analyse de départ)
  x[
    , Age := as.integer(floor(difftime(t1, DatNais, units = "days") / 365.25))
  ]
  # Inclure tous les âges dans 'catage'
  catage <- c(catage, max(x[, Age]) + 1)
  # Création de la catégorie d'age
  x[
    , `:=` (CatAge = cut(Age, catage, right = FALSE),
            # Supprimer les colonnes non nécessaires
            Age = NULL,
            DatNais = NULL)
  ]
  # Ajout du Sexe
  x <- merge(x, unique(BDsex),
             by.x = "BenBanls", by.y = xcols_[[1]])

  # Statistiques descriptives
  # Par Age et Sexe
  nAS <- x[, .(n = .N), keyby = .(CatAge, Sexe)]  # Nombre d'observation par groupe
  pond_AS <- x[, .(Moyenne = mean(nfrac),
                   `Écart-type` = sd(nfrac),
                   Min = min(nfrac),
                   P5 = quantile(nfrac, 0.05),
                   P10 = quantile(nfrac, 0.10),
                   Q1 = quantile(nfrac, 0.25),
                   `Médiane` = median(nfrac),
                   Q3 = quantile(nfrac, 0.75),
                   P90 = quantile(nfrac, 0.90),
                   P95 = quantile(nfrac, 0.95),
                   Max = max(nfrac)),
               keyby = .(CatAge, Sexe)]
  pond_AS <- merge(pond_AS, nAS)
  # Par Age
  nA <- x[, .(n = .N), keyby = .(CatAge)]  # Nombre d'observation par groupe
  pond_A <- x[, .(Moyenne = mean(nfrac),
                  `Écart-type` = sd(nfrac),
                  Min = min(nfrac),
                  P5 = quantile(nfrac, 0.05),
                  P10 = quantile(nfrac, 0.10),
                  Q1 = quantile(nfrac, 0.25),
                  `Médiane` = median(nfrac),
                  Q3 = quantile(nfrac, 0.75),
                  P90 = quantile(nfrac, 0.90),
                  P95 = quantile(nfrac, 0.95),
                  Max = max(nfrac)),
              keyby = .(CatAge)]
  pond_A <- merge(pond_A, nA)
  # Par Sexe
  nS <- x[, .(n = .N), keyby = .(Sexe)]  # Nombre d'observation par groupe
  pond_S <- x[, .(Moyenne = mean(nfrac),
                  `Écart-type` = sd(nfrac),
                  Min = min(nfrac),
                  P5 = quantile(nfrac, 0.05),
                  P10 = quantile(nfrac, 0.10),
                  Q1 = quantile(nfrac, 0.25),
                  `Médiane` = median(nfrac),
                  Q3 = quantile(nfrac, 0.75),
                  P90 = quantile(nfrac, 0.90),
                  P95 = quantile(nfrac, 0.95),
                  Max = max(nfrac)),
              keyby = .(Sexe)]
  pond_S <- merge(pond_S, nS)


  # Distribution des max
  # Valeurs possibles
  seq_class <- 0:(max(cumulpond_BD$nfrac) + classes)
  # Création et affichage des classes
  if(is.null(maxclasses)){
    seq_class <- unique(cut(seq_class,
                            seq(0, max(cumulpond_BD$nfrac) + classes, classes),  # bornes des classes
                            right = FALSE))  # inclu gauche, droit exclu
  } else {
    seq_class <- unique(cut(seq_class,
                            c(seq(0, maxclasses, classes), ceiling(max(cumulpond_BD$nfrac))),  # bornes des classes
                            right = FALSE))
  }
  seq_class <- seq_class[!is.na(seq_class)]  # supprimer NA
  dist_max_arrondi <- cumulpond_BD[, nfrac := round(nfrac)]
  dist_max_arrondi <- dist_max_arrondi[, .(Freq = .N), nfrac]
  names(dist_max_arrondi) <- c("Valeur", "Freq")
  seq_arrondi <- 0:ceiling(max(cumulpond_BD$nfrac))
  if(length(seq_arrondi[!seq_arrondi %in% dist_max_arrondi$Valeur]) != 0){
    dist_max_arrondi <- rbind(dist_max_arrondi,
                              data.table(Valeur = seq_arrondi[!seq_arrondi %in% dist_max_arrondi$Valeur],
                                         Freq = 0))  # Ajout des valeurs de max manquantes
  }
  setorder(dist_max_arrondi, Valeur)  # Tri par classes croissant
  dist_max_arrondi[, Pourcentage := Freq / sum(Freq) * 100]  # pourcentage
  dist_max_arrondi[, Cumul := cumsum(Pourcentage)]  # pourcentage cumulé
  # Avec classes
  dist_max <- cumulpond_BD[, n_class := cut(nfrac,
                                            {if(is.null(maxclasses)){
                                                seq(0, max(cumulpond_BD$nfrac) + classes, classes)
                                              } else {
                                                c(seq(0, maxclasses, classes), ceiling(max(cumulpond_BD$nfrac)))
                                              }},  # bornes des classes
                                            right = FALSE)]  # gauche inclu, droite exclu
  dist_max <- dist_max[, .(Freq = .N), n_class]  # tableau de fréquence
  names(dist_max) <- c("Valeur", "Freq")  # Noms des colonnes
  if(length(seq_class[!seq_class %in% dist_max$Valeur]) != 0){
    dist_max <- rbind(dist_max, data.table(Valeur = seq_class[!seq_class %in% dist_max$Valeur],
                                           Freq = 0))  # Ajout des valeurs de max manquantes
  }
  setorder(dist_max, Valeur)  # Tri par classes croissant
  dist_max[, Pourcentage := Freq / sum(Freq) * 100]  # pourcentage
  dist_max[, Cumul := cumsum(Pourcentage)]  # pourcentage cumulé

  ## Résultat final
  # Liste des éléments
  cumulpond <- list(BD = cumulpond_BD,
                    stats = cumulpond_stats,
                    stats_group = list(`CatAge & Sexe` = pond_AS,
                                       CatAge = pond_A,
                                       Sexe = pond_S),
                    dist_max = dist_max,
                    dist_max_arrondi = dist_max_arrondi)

}