polymedic: Calculs de la polypharmacie

Documented in polysimul

#' Polypharmacie Simultanée
#'
#' \code{polysimul} analyse quel est le nombre de médicaments différents consommés par un usager au jour \emph{T}. Analyse pour tous les jours compris dans l'intervalle [\code{t1}, \code{t2}].
#'
#' \strong{\code{t2}} est facultatif. S'il n'est pas défini, \code{t2} prend la valeur de \code{t1} pour faire l'analyse d'une journée, sinon l'analyse se fera pour tous les jours allant de \code{t1} à \code{t2}.
#'
#' @param x Base de données créée par \code{polyBD()}.
#' @param t1 "AAAA-MM-JJ". Jour à analyser.
#' @param t2 Facultatif. "AAAA-MM-JJ". Voir détails.
#' @param catage Borne(s) inférieure(s) des groupes d'âge.
#'
#' @import data.table
#' @import plyr
#' @export
polysimul <- function(x, t1, t2, catage = c(60,70,80,90)){

  ## Noms des colonnes
  # Conserver
  xcols_ <- copy(names(x))
  # Modifier
  names(x)[[1]] <- "BenBanls"
  names(x)[[2]] <- "DIN"
  # Liste des ID
  xID <- unique(x$BenBanls)
  # BD Age et Sexe
  BDage <- attributes(x)$Naissance
  BDsex <- attributes(x)$Sexe

  # T2 est manquant
  if(missing(t2)){t2 <- t1} #prend la valeur de t1 : analyse de 1 journée

  ## Afficher le nombre de médicaments pour chaque personne pour chaque jour entre t1 et t2
  DT <- lapply(
    as.Date(t1):as.Date(t2), #les jours à analyser
    function(j){
      # Sélection du jour x où la personne consomme le médicament
      simul <- x[PerDebut <= j & PerFinAjust >= j]
      # Indique le nombre de médicament consommé ce jour-là
      simul <- simul[, .(n = .N), .(BenBanls)]
      # Changer le nom de la colonne pour le jour analysé
      setnames(simul, "n", as.character(as.Date(j, "1970-01-01")))

      # Résultat fonction
      simul
    }
  )
  # Joindre tous les jours analysés à la suite les uns les autres pour une même personne
  DT <- Reduce(function(x,y){merge(x, y, all = TRUE)}, DT)

  # Changer NA par 0
  for(j in names(DT)[names(DT) != "BenBanls"]){ #pour toutes les colonnes sauf BenBanls
    set(DT, which(is.na(DT[[j]])), j, 0)
  }
  # Modifier les colonnes en lignes
  DT <- melt(DT, id.vars = 1, variable.name = "Date", value.name = "n")
  setkey(DT, BenBanls, Date)


  ## Ajout des ID manquants - sans consommation
  # Liste des ID manquants
  xID <- xID[!xID %in% DT$BenBanls]
  # Tableau des ID manquants avec leur consommation du jour = 0
  if(length(xID) != 0){
    xID <- data.table(BenBanls = rep(xID, each = length(as.Date(t1):as.Date(t2))),
                      Date = as.Date(t1):as.Date(t2),
                      n = 0L)
    # Convertir les variables pour fusion des tableaux
    xID[, Date := as.character(as.Date(Date, origin = "1970-01-01"))]
    DT[, Date := as.character(Date)]
    # Fusion des tableaux
    DT <- rbind(DT, xID)
  }

  ## Nombre d'usagers
  nID <- length(unique(DT$BenBanls))


  ## Résultats
  # simul1 summary - statistiques
  simul1_summary <- DT[, .(Moyenne = mean(n),
                           `Écart-type` = sd(n),
                           Min = min(n),
                           P5 = quantile(n, 0.05),
                           P10 = quantile(n, 0.10),
                           Q1 = quantile(n, 0.25),
                           `Médiane` = median(n),
                           Q3 = quantile(n, 0.75),
                           P90 = quantile(n, 0.90),
                           P95 = quantile(n, 0.95),
                           Max = max(n)),
                       .(BenBanls)]
  simul1_summary <- simul1_summary[, c(.SD[, 1], lapply(.SD[, 2:12], as.double))]
  # BD d'analyse
  simul1_BD <- copy(DT)
  
  # Nombre de jours avec au moins x médicaments consommés
  simul_nmed_BD <- as.data.table(rbind.fill(lapply(1:max(simul1_BD$n), function(x) {
    dt <- copy(simul1_BD)
    dt[, Medicaments := paste(">=", x)]
    dt[, medic := 0][n >= x, medic := 1][, .(n = sum(medic)), keyby = .(Medicaments, BenBanls)]
  })))
  # Statistiques pour tous les individus
  simul_nmed <- 
    simul_nmed_BD[
    , .(Moyenne = mean(n),
        `Écart-type` = sd(n),
        Min = min(n),
        P5 = quantile(n, 0.05),
        P10 = quantile(n, 0.10),
        Q1 = quantile(n, 0.25),
        `Médiane` = median(n),
        Q3 = quantile(n, 0.75),
        P90 = quantile(n, 0.90),
        P95 = quantile(n, 0.95),
        Max = max(n),
        n = nID),
    ,.(Medicaments)
  ]
  simul_nmed[["Medicaments"]] <- factor(simul_nmed[["Medicaments"]],
                                        levels = paste(">=", 1:200))
  setorder(simul_nmed, Medicaments)
  # Distribution des valeurs de simul_nmed
  dist_nmed <- as.data.table(rbind.fill(lapply(split(simul_nmed_BD, by = "Medicaments"), function(x) {
    seq_nmed <- 0:max(x$n)
    dt <- x[, .(Freq = .N), keyby = .(Medicaments, n)]
    if (length(seq_nmed[!seq_nmed %in% unique(x[, n])]) != 0) {
      dt <- rbind(dt, data.table(Medicaments = unique(x$Medicaments),
                                 n = seq_nmed[!seq_nmed %in% unique(x[, n])],
                                 Freq = 0))
    }
    setnames(dt, "n", "nbJours")
    setorder(dt, nbJours)
    dt[, Pourcentage := Freq / sum(Freq) * 100]
    dt[, Cumul := cumsum(Pourcentage)]
  })))
  
  # Statistiques par Strate
  # Ajout de la date de naissance
  stats_group <- merge(simul_nmed_BD, unique(BDage),
                       by.x = "BenBanls", by.y = xcols_[[1]])
  # Convertir date de naissance en Age (selon la date d'analyse de départ)
  stats_group[
    , Age := as.integer(floor(difftime(t1, DatNais, units = "days") / 365.25))
  ]
  # Inclure tous les âges dans 'catage'
  catage <- c(catage, max(stats_group[, Age]) + 1)
  # Création de la catégorie d'age
  stats_group[
    , `:=` (CatAge = cut(Age, catage, right = FALSE),
            # Supprimer les colonnes non nécessaires
            Age = NULL,
            DatNais = NULL)
  ]
  # Ajout du Sexe
  stats_group <- merge(stats_group, unique(BDsex),
             by.x = "BenBanls", by.y = xcols_[[1]])
  # Stats par groupe d'âge et Sexe
  stats_group_AS <- stats_group[
    , .(Moyenne = mean(n),
        `Écart-type` = sd(n),
        Min = min(n),
        P5 = quantile(n, 0.05),
        P10 = quantile(n, 0.10),
        Q1 = quantile(n, 0.25),
        `Médiane` = median(n),
        Q3 = quantile(n, 0.75),
        P90 = quantile(n, 0.90),
        P95 = quantile(n, 0.95),
        Max = max(n),
        n = stats_group[, .N, .(CatAge, Sexe)]$N),
    .(CatAge, Sexe)
  ]
  # Stats par groupe d'âge
  stats_group_A <- stats_group[
    , .(Moyenne = mean(n),
        `Écart-type` = sd(n),
        Min = min(n),
        P5 = quantile(n, 0.05),
        P10 = quantile(n, 0.10),
        Q1 = quantile(n, 0.25),
        `Médiane` = median(n),
        Q3 = quantile(n, 0.75),
        P90 = quantile(n, 0.90),
        P95 = quantile(n, 0.95),
        Max = max(n),
        n = stats_group[, .N, .(CatAge)]$N),
    .(CatAge)
  ]
  # Stats par Sexe
  stats_group_S <- stats_group[
    , .(Moyenne = mean(n),
        `Écart-type` = sd(n),
        Min = min(n),
        P5 = quantile(n, 0.05),
        P10 = quantile(n, 0.10),
        Q1 = quantile(n, 0.25),
        `Médiane` = median(n),
        Q3 = quantile(n, 0.75),
        P90 = quantile(n, 0.90),
        P95 = quantile(n, 0.95),
        Max = max(n),
        n = stats_group[, .N, .(Sexe)]$N),
    .(Sexe)
  ]
  
  names(DT)[[1]] <- xcols_[[1]]  # colname initial
  names(simul_nmed_BD)[[2]] <- xcols_[[1]]  # colname initial

  
  
  # Distribution des valeurs
  dist_values <- function(x,  # BD
                          col,  # nom de la colonne
                          method = "int"  # 'int' si les valeurs sont des entiers, sinon 'double' -> arrondissement en entier
                          ){

    setnames(x, col, "xcolx")  # Les valeurs à afficher -> 0 à valeur maximale
    if(method == "double"){
      seq_values <- 0:round(max(x$xcolx))
      dist <- copy(x[, xcolx := round(xcolx)])
      dist <- dist[, .(Freq = .N), xcolx]
    } else {
      seq_values <- 0:max(x$xcolx)  # toutes les valeurs à afficher -> 0 à valeur max
      dist <- x[, .(Freq = .N), xcolx]  # fréquence pour chaque valeur
    }
    if(length(seq_values[!seq_values %in% dist$xcolx]) != 0){
      dist <- rbind(dist, data.table(xcolx = seq_values[!seq_values %in% dist$xcolx],
                                     Freq = 0))  # ajouter les valeurs manquante -> Freq = 0
    }
    setorder(dist, xcolx)  # tri croissant par valeur
    dist[, Pourcentage := Freq / sum(Freq) * 100]  # pourcentage
    dist[, Cumul := cumsum(Pourcentage)] #pourcentage cumulé
    # Arranger noms des colonnes
    setnames(dist, "xcolx", "Valeur") # dist
    setnames(x, "xcolx", col) # x

    # Résultat
    dist

  }  # fonction
  # Min
  dist_min <- dist_values(simul1_summary, "Min")
  # Moyenne
  dist_moy <- dist_values(simul1_summary, "Moyenne", "double")
  # Médiane
  dist_med <- dist_values(simul1_summary, "Médiane")
  # Max
  dist_max <- dist_values(simul1_summary, "Max")


  # Statistiques sur chaque colonne
  simul2_summary <- list(
    Moyenne = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], mean)],
    `Écart-type` = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], sd)],
    Min = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], min)],
    P5 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.05)],
    P10 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.10)],
    Q1 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.25)],
    `Médiane` = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], median)],
    Q3 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.75)],
    P90 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.90)],
    P95 = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], quantile, probs = 0.95)],
    Max = simul1_summary[, lapply(.SD[, 2:ncol(simul1_summary)], max)],
    n = nID
  )
  # Statistiques groupé par Sexe ou CatAge
  # Par catage et Sexe
  # Ajout de la date de naissance
  simul1_summary <- merge(simul1_summary, unique(attributes(x)$Naissance),
                          by.x = "BenBanls", by.y = xcols_[[1]])
  # Convertir date de naissance en Age (selon la date d'analyse de départ)
  simul1_summary[
    , Age := as.integer(floor(difftime(t1, DatNais, units = "days") / 365.25))
  ]
  # Création de la catégorie d'age
  simul1_summary[
    , `:=` (CatAge = cut(Age, catage, right = FALSE),
            # Supprimer les colonnes non nécessaires
            Age = NULL,
            DatNais = NULL)
  ]
  # Ajout du Sexe
  simul1_summary <- merge(simul1_summary, unique(attributes(x)$Sexe),
                          by.x = "BenBanls", by.y = xcols_[[1]])
  setnames(simul1_summary, "BenBanls", xcols_[[1]])

  # Statistiques descriptive
  # Sur CatAge et Sexe
  simul2_summary_AS <- list(
    Moyenne = simul1_summary[, lapply(.SD[, 2:12], mean), keyby = .(CatAge, Sexe)],
    `Écart-type` = simul1_summary[, lapply(.SD[, 2:12], sd), keyby = .(CatAge, Sexe)],
    Min = simul1_summary[, lapply(.SD[, 2:12], min), keyby = .(CatAge, Sexe)],
    P5 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.05), keyby = .(CatAge, Sexe)],
    P10 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.10), keyby = .(CatAge, Sexe)],
    Q1 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.25), keyby = .(CatAge, Sexe)],
    `Médiane` = simul1_summary[, lapply(.SD[, 2:12], median), keyby = .(CatAge, Sexe)],
    Q3 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.75), keyby = .(CatAge, Sexe)],
    P90 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.90), keyby = .(CatAge, Sexe)],
    P95 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.95), keyby = .(CatAge, Sexe)],
    Max = simul1_summary[, lapply(.SD[, 2:12], max), keyby = .(CatAge, Sexe)],
    n = simul1_summary[, .N, .(CatAge, Sexe)]$N
  )
  # Sur Age
  simul2_summary_A <- list(
    Moyenne = simul1_summary[, lapply(.SD[, 2:12], mean), keyby = .(CatAge)],
    `Écart-type` = simul1_summary[, lapply(.SD[, 2:12], sd), keyby = .(CatAge)],
    Min = simul1_summary[, lapply(.SD[, 2:12], min), keyby = .(CatAge)],
    P5 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.05), keyby = .(CatAge)],
    P10 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.10), keyby = .(CatAge)],
    Q1 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.25), keyby = .(CatAge)],
    `Médiane` = simul1_summary[, lapply(.SD[, 2:12], median), keyby = .(CatAge)],
    Q3 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.75), keyby = .(CatAge)],
    P90 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.90), keyby = .(CatAge)],
    P95 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.95), keyby = .(CatAge)],
    Max = simul1_summary[, lapply(.SD[, 2:12], max), keyby = .(CatAge)],
    n = simul1_summary[, .N, .(CatAge)]$N
  )
  # Sur Sexe
  simul2_summary_S <- list(
    Moyenne = simul1_summary[, lapply(.SD[, 2:12], mean), keyby = .(Sexe)],
    `Écart-type` = simul1_summary[, lapply(.SD[, 2:12], sd), keyby = .(Sexe)],
    Min = simul1_summary[, lapply(.SD[, 2:12], min), keyby = .(Sexe)],
    P5 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.05), keyby = .(Sexe)],
    P10 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.10), keyby = .(Sexe)],
    Q1 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.25), keyby = .(Sexe)],
    `Médiane` = simul1_summary[, lapply(.SD[, 2:12], median), keyby = .(Sexe)],
    Q3 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.75), keyby = .(Sexe)],
    P90 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.90), keyby = .(Sexe)],
    P95 = simul1_summary[, lapply(.SD[, 2:12], quantile, probs = 0.95), keyby = .(Sexe)],
    Max = simul1_summary[, lapply(.SD[, 2:12], max), keyby = .(Sexe)],
    n = simul1_summary[, .N, .(Sexe)]$N
  )


  ## Résultat final
  # Liste des éléments calculés
  simul <- list(id_BD = simul1_BD,
                id_stats = simul1_summary,
                pop_stats = simul2_summary,
                pop_stats_group = list(`CatAge & Sexe` = simul2_summary_AS,
                                       CatAge = simul2_summary_A,
                                       Sexe = simul2_summary_S),
                dist_values = list(Min = dist_min,
                                   Moy = dist_moy,
                                   Med = dist_med,
                                   Max = dist_max),
                stats_nmed = list(BD = simul_nmed_BD,
                                  stats = simul_nmed),
                stats_nmed_group = list(`CatAge & Sexe` = stats_group_AS,
                                        CatAge = stats_group_A,
                                        Sexe = stats_group_S),
                dist_nmed = dist_nmed)
  # Trier les catégories d'âge
  attributes(simul)$catage <- simul$pop_stats_group$CatAge$Moyenne$CatAge
  # Afficher
  simul
  
}