Une fonction qu'un autre chargé peut exécuter — R et le tidyverse pour les données de programme — Analyse de données

Le test de la transmission

Une analyse est terminée quand quelqu’un d’autre peut l’exécuter sur l’export du trimestre prochain sans modifier le code. Non pas « sans trop le modifier » — sans le modifier.

Ce critère élimine les trois habitudes que tout script R accumule : un chemin écrit pour une machine, un seuil changé à la main entre deux exécutions, et une étape qui ne fonctionne que si l’on a déjà exécuté les lignes précédentes dans le bon ordre.

Du script à la fonction

Le nettoyage des leçons précédentes, en une fonction :

# R/clean_register.R

clean_register <- function(muac, gam_mm = 125, sam_mm = 115) {
  stopifnot(is.data.frame(muac))
  required <- c("child_id", "commune", "muac_mm", "oedema")
  missing <- setdiff(required, names(muac))
  if (length(missing) > 0) {
    stop("Colonnes absentes de l'entree : ", paste(missing, collapse = ", "))
  }

  muac |>
    dplyr::mutate(
      # Un PB sous 40 est en centimetres dans un champ en millimetres. Applique
      # seulement sous un seuil qu'aucune mesure plausible n'atteint.
      muac_mm = dplyr::if_else(muac_mm < 40, muac_mm * 10L, muac_mm),
      muac_mm = dplyr::if_else(dplyr::between(muac_mm, 80L, 220L), muac_mm, NA_integer_),
      oedema  = dplyr::recode(
        tolower(trimws(oedema)),
        "true" = TRUE, "y" = TRUE, "yes" = TRUE,
        "false" = FALSE, "n" = FALSE, "no" = FALSE,
        .default = NA
      ),
      assessed = !is.na(muac_mm) | !is.na(oedema),
      gam = assessed & (muac_mm < gam_mm | oedema),
      sam = assessed & (muac_mm < sam_mm | oedema)
    )
}

Quatre aspects de cette forme sont délibérés.

Elle prend un tableau de données et en renvoie un. Aucune lecture de fichier, aucune écriture, aucun <<-. C’est ce qui la rend testable et ce qui permet à un notebook et à un script de la partager.

Les seuils sont des arguments avec valeurs par défaut. Un seuil qu’on ne change qu’en éditant une ligne finira par être changé et non remis. Un seuil qui est un argument est visible dans l’appel qui a produit la sortie.

Elle valide son entrée. Un stop() nommant les colonnes absentes vaut mieux qu’un NULL se propageant dans un mutate et refaisant surface en erreur sans rapport deux fonctions plus loin.

Rien n’est groupé en sortie. Un tableau groupé qui s’échappe d’une fonction est l’erreur de la leçon sur le regroupement, arrivant par une porte nouvelle.

`stopifnot()` pour ce qui ne peut pas arriver

stopifnot(
  "chaque ligne exige un identifiant" = !any(is.na(muac$child_id)),
  "le PB doit etre plausible apres nettoyage" =
    all(is.na(muac$muac_mm) | dplyr::between(muac$muac_mm, 80, 220))
)

Les expressions nommées deviennent le message d’erreur, ce qui fait la différence entre un échec utile et Error: ... is not TRUE.

Employez stop() pour « cette entrée n’est pas celle qu’on m’avait promise » et stopifnot() pour « ceci ne peut pas arriver à moins que le code soit faux ». La distinction compte pour qui lira la défaillance à sept heures du matin.

Passer un nom de colonne : l’évaluation tidy

C’est le seul point de R qui surprenne les personnes venant de Python, et il n’y a qu’une chose à apprendre.

# Ne fonctionne pas.
rate_by <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = group_col)
}

rate_by(muac, commune)
#> Error: object 'commune' not found

Les verbes dplyr évaluent leurs arguments à l’intérieur du tableau de données. Le group_col de votre fonction est une variable de la fonction, non une colonne : dplyr cherche donc une colonne nommée group_col et échoue.

Le correctif tient en un opérateur :

rate_by <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = {{ group_col }})
}

rate_by(muac, commune)

{{ }} — l’« embrassement » — signifie prends ce que l’appelant a écrit et évalue-le dans les données. C’est la réponse pour presque toute fonction que vous écrirez autour de dplyr.

Pour un nom de colonne arrivant sous forme de chaîne, ce que donne un argument de ligne de commande :

rate_by_name <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = dplyr::all_of(group_col))
}

rate_by_name(muac, "commune")

all_of() lève une erreur si la colonne n’existe pas ; any_of() la saute en silence. Pour un script dont la sortie est rapportée, all_of() est celui qu’il vous faut.

Pour nommer une colonne de sortie à partir d’un argument :

count_as <- function(df, group_col, out_name) {
  df |> dplyr::summarise("{out_name}" := dplyr::n(), .by = {{ group_col }})
}

La forme "{...}" := est la syntaxe glue à l’intérieur d’un verbe du tidyverse. Elle mérite d’être reconnue ; vous en aurez rarement besoin.

La fonction d’indicateur

indicator_table <- function(muac, gam_mm = 125, sam_mm = 115) {
  clean_register(muac, gam_mm, sam_mm) |>
    dplyr::summarise(
      screened  = dplyr::n(),
      assessed  = sum(assessed, na.rm = TRUE),
      gam_cases = sum(gam, na.rm = TRUE),
      sam_cases = sum(sam, na.rm = TRUE),
      .by = commune
    ) |>
    dplyr::mutate(
      assessment_rate = round(assessed / screened, 3),
      gam_rate = round(gam_cases / assessed, 3),
      sam_rate = round(sam_cases / assessed, 3)
    ) |>
    dplyr::arrange(dplyr::desc(gam_rate))
}

Notez ce qu’elle ne fait pas : elle ne supprime pas les communes à faible taux d’évaluation. En supprimer une changerait le dénominateur du district et rien dans la sortie ne le dirait. Le signalement revient à l’appelant, qui peut voir le taux.

Tester la partie qui calcule

Une fois le calcul devenu une fonction prenant un tableau, un test tient en trois lignes :

# tests/testthat/test-indicator-table.R
test_that("le denominateur exclut les enfants non mesures", {
  muac <- tibble::tibble(
    child_id = c("a", "b", "c"),
    commune  = "X",
    muac_mm  = c(110L, 130L, NA_integer_),
    oedema   = c("false", "false", NA_character_)
  )

  out <- indicator_table(muac)

  expect_equal(out$screened, 3)
  expect_equal(out$assessed, 2)     # l'enfant non mesure n'est pas au denominateur
  expect_equal(out$gam_rate, 0.5)
})

Ce test encode la décision de la leçon sur le regroupement — quels enfants entrent au dénominateur — sous une forme qui échoue si quelqu’un la change. À écrire pour tout indicateur dont la définition a fait débat.

usethis::use_testthat() met en place le répertoire ; devtools::test() l’exécute.

Le script en ligne de commande

#!/usr/bin/env Rscript
# scripts/indicator_table.R

suppressPackageStartupMessages({
  library(optparse)
  library(readr)
})

source(here::here("R", "clean_register.R"))
source(here::here("R", "indicator_table.R"))

option_list <- list(
  make_option("--input", type = "character", help = "CSV du registre de depistage."),
  make_option("--output", type = "character", help = "Repertoire des tableaux."),
  make_option("--gam-threshold", type = "integer", default = 125L,
              help = "Seuil de PB en mm pour la MAG [defaut %default].")
)

opt <- parse_args(OptionParser(option_list = option_list))
if (is.null(opt$input) || is.null(opt$output)) {
  stop("--input et --output sont requis.")
}

muac <- read_register(opt$input)
message("Lu ", nrow(muac), " lignes depuis ", opt$input)

table <- indicator_table(muac, gam_mm = opt$`gam-threshold`)

thin <- table$commune[table$assessment_rate < 0.8]
if (length(thin) > 0) {
  warning("Taux d'evaluation sous le plancher : ", paste(thin, collapse = ", "))
}

dir.create(opt$output, recursive = TRUE, showWarnings = FALSE)
destination <- file.path(opt$output, "gam_by_commune.csv")
write_csv(table, destination)
message("Ecrit ", destination, " (", nrow(table), " communes)")

Rscript scripts/indicator_table.R \
  --input data/raw/muac-screening-artibonite-2024.v1.csv \
  --output outputs/tables

message() plutôt que print() : les messages partent sur la sortie d’erreur, si bien que la narration du script ne se mêle pas à une sortie redirigée. warning() plutôt que message() pour les communes minces, car c’est une condition qu’un appelant peut vouloir durcir avec options(warn = 2).

Consigner ce qui a produit la sortie

run_metadata <- function(opt) {
  list(
    generated_at  = format(Sys.time(), "%Y-%m-%dT%H:%M:%S%z"),
    input         = opt$input,
    gam_threshold = opt$`gam-threshold`,
    r_version     = as.character(getRversion()),
    dplyr         = as.character(packageVersion("dplyr"))
  )
}

jsonlite::write_json(run_metadata(opt), file.path(opt$output, "run.json"),
                     auto_unbox = TRUE, pretty = TRUE)

Six mois plus tard, « quel seuil a produit ce tableau » a une réponse qui ne dépend de la mémoire de personne.

Puis le véritable test, issu de la première leçon :

rm -rf outputs/
Rscript scripts/indicator_table.R --input data/raw/muac.csv --output outputs/tables
git status

Si cela ne restaure pas les sorties à l’identique, quelque chose n’est pas reproductible.

Où ce cours s’arrête

Vous savez construire un projet qui se rouvre un an plus tard avec ses versions de paquets consignées, lire n’importe quel export sans le corrompre, conserver les étiquettes qu’un fichier d’enquête transporte, mettre les catégories dans l’ordre qu’exige un rapport, agréger vers un numérateur et un dénominateur issus d’un seul appel, restructurer un groupe répété aplati, et remettre le résultat à quelqu’un d’autre sous forme de script plutôt que de service rendu.

Ce que ce cours n’a délibérément pas enseigné, c’est quoi calculer — quel indicateur, quel dénominateur, quel seuil, et comment le défendre. C’est l’objet de Fondamentaux de l’analyse de données, et les cours sectoriels du module Analyses sectorielles de la feuille de route le poussent plus loin, jusqu’aux classifications que votre cluster vous impose.

Si votre équipe travaille en Python plutôt qu’en R, Python pour les données de programme couvre le même terrain dans ce langage — et le tableau des renversements de la leçon 6 est le moyen le plus rapide de passer de l’un à l’autre.