models-diagnostics.Rmd

---
title: "Models Diagnostics"
author: "Mickaël Canouil, *Ph.D.* (<mickael.canouil@cnrs.fr>)"
params:
  working_directory: !r tempdir()
  pheno_file: "internal"
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# R setup

```{r setup}
output_directory <- file.path(params[["working_directory"]], "models-diagnostics")
dir.create(output_directory, recursive = TRUE, showWarnings = FALSE, mode = "0775")

library(eggla)
library(growthcleanr)

library(data.table)
library(nlme)
library(stats)

library(grDevices)
library(ggplot2)
library(ggtext)
library(patchwork)
library(scales)

library(gt)
library(performance)
library(see)
library(qqplotr)

library(future)
library(future.apply)
library(future.callr)
```

```{r future-setup}
plan("callr", workers = availableCores())
message(sprintf("Number of workers: %d", nbrOfWorkers()))
```

```{r ggplot-setup}
options(
  ggplot2.discrete.colour = function(...) scale_colour_viridis_d(..., begin = 0.15, end = 0.85),
  ggplot2.discrete.fill = function(...) scale_fill_viridis_d(..., begin = 0.15, end = 0.85),
  ggplot2.continuous.colour = function(...) scale_colour_viridis_c(..., begin = 0.15, end = 0.85),
  ggplot2.continuous.fill = function(...) scale_fill_viridis_c(..., begin = 0.15, end = 0.85)
)
theme_set(theme_minimal())
theme_update(
  plot.title.position = "plot",
  plot.caption.position = "plot",
  plot.title = element_markdown(),
  plot.subtitle = element_markdown(face = "italic"),
  plot.caption = element_markdown(face = "italic"),
  axis.title.x = element_markdown(),
  axis.text.x = element_markdown(),
  axis.title.y = element_markdown(),
  axis.text.y = element_markdown()
)
```

# Data

```{r input-data}
if (params[["pheno_file"]] == "internal") {
  pheno_dt <- setDT(get(data("bmigrowth", package = "eggla")))
} else {
  pheno_dt <- fread(params[["pheno_file"]])
}
```

## Tidy data for Daymont's quality control 

- `age` - age in years (double).
- `agedays` - age in days (integer).
- `WEIGHTKG` - weight in kilograms (double).
- `HEIGHTCM` - height in centimetres (double).

```{r format-data}
is_male_zero <- FALSE
pheno_dt[
  j = `:=`(
    "agedays" = floor(age * 365.25), # convert to age in days and as integers ...
    "WEIGHTKG" = as.numeric(weight),
    "HEIGHTCM" = as.numeric(height),
    "sex_daymont" = {
      if (is_male_zero) {
        as.character(sex)
      } else {
        c("0" = "1", "1" = "0")[as.character(sex)]
      }
    }
  )
]
```

## Add Daymont's quality control tags

```{r add-daymont-data}
visits_long <- melt(
  data = pheno_dt,
  id.vars = c("ID", "age", "sex", "agedays", "sex_daymont"),
  measure.vars = c("WEIGHTKG", "HEIGHTCM"),
  variable.name = "param",
  value.name = "measurement"
)[
  j = clean := cleangrowth(
    subjid = ID,
    param = param,
    agedays = agedays,
    sex = sex_daymont,
    measurement = measurement,
    quietly = FALSE
  )
]
```

## Compute/Update BMI

- BMI based on raw measures

    ```{r visits-raw}
    visits_raw <- dcast(
      data = visits_long,
      formula = ... ~ param,
      value.var = "measurement"
    )[
      j = "bmi" := WEIGHTKG / (HEIGHTCM / 100)^2 # recompute bmi based on QC variables
    ][
      !is.na(bmi) # exclude missing BMI related to measurements exclusion
    ]
    ```

- BMI based on Clean measures

    ```{r visits-qc}
    visits_clean <- dcast(
      data = visits_long[clean %in% "Include"], # Exclude all flags
      formula = ... ~ param,
      value.var = "measurement"
    )[
      j = "bmi" := WEIGHTKG / (HEIGHTCM / 100)^2 # recompute bmi based on QC variables
    ][
      !is.na(bmi) # exclude missing BMI related to measurements exclusion
    ]
    ```

# Model diagnostics 

*With Daymont's code for sex, Female: 1 and Male: 0.*

## Setup all models to be tested

```{r models-setup}
random_effect <- c(
  "~ gsp(age, knots = c(5.5, 11), degree = rep(1, 3), smooth = rep(0, 2)) | ID",
  "~ 1 | ID",
  "~ gsp(age, knots = c(2, 8, 12), degree = rep(3, 4), smooth = rep(2, 3)) | ID",
  "~ gsp(age, knots = c(2, 8, 12), degree = rep(3, 4), smooth = rep(2, 3))[, 1:3] | ID",
  "~ gsp(age, knots = c(2, 8, 12), degree = rep(2, 4), smooth = rep(2, 3)) | ID",
  "~ gsp(age, knots = c(2, 8, 12), degree = rep(1, 4), smooth = rep(2, 3)) | ID",
  "~ 1 | ID",
  "~ poly(age, degree = 3) | ID",
  "~ poly(age, degree = 2) | ID",
  "~ age | ID",
  "~ 1 | ID"
)

fixed_effect <- c(
  rep("log(bmi) ~ gsp(age, knots = c(5.5, 11), degree = rep(1, 3), smooth = rep(0, 2))", 2),
  rep("log(bmi) ~ gsp(age, knots = c(2, 8, 12), degree = rep(3, 4), smooth = rep(2, 3))", 5),
  rep("log(bmi) ~ poly(age, degree = 3)", 4)
)

all_models <- setDT(
  expand.grid(
    sex = c(0, 1),
    dataset = c(
      "visits_raw",
      "visits_clean",
      "visits_clean[between(age, 1, 18)]"
    ),
    correlation_structure = c("NULL", "corCAR1(form = ~ 1 | ID)")[2], # Use AR(1) correlation structure (default is without)
    models = sprintf("fixed = %s, random = %s", fixed_effect, random_effect),
    default = paste0(
      c(
        "na.action = na.omit",
        "method = \"ML\"",
        "control = lmeControl(opt = \"optim\", maxIter = 500, msMaxIter = 500)"
      ),
      collapse = ", "
    ),
    lme = list(NULL),
    results = list("Good"),
    perf_check_model = list(NULL),
    perf_model_performance = list(NULL),
    stringsAsFactors = FALSE
  )
)
```

## Run all models

*In parallel using `future`.*

```{r models-try}
all_models[["lme"]] <- future_mapply(
  FUN = function(models, dataset, sex, correlation_structure, default) {
    m <- sprintf(
      "list(lme(%s, data = %s[sex_daymont == %s], correlation = %s, %s))",
      models,
      dataset,
      sex,
      correlation_structure,
      default
    )
    tryCatch(
      expr = eval(parse(text = m)),
      error = function(e) sprintf("Error: %s", e$message),
      warning = function(w) sprintf("Warning: %s", w$message)
    )
  },
  models = all_models[["models"]],
  dataset = all_models[["dataset"]],
  sex = all_models[["sex"]],
  correlation_structure = all_models[["correlation_structure"]],
  default = all_models[["default"]],
  future.globals = c("visits_raw", "visits_clean"),
  future.packages = c("eggla", "nlme", "data.table"),
  USE.NAMES = FALSE
)

all_models[["is_good"]] <- future_mapply(
  FUN = function(lme) inherits(summary(lme), "lme"),
  lme = all_models[["lme"]],
  future.globals = FALSE,
  future.packages = c("nlme"),
  USE.NAMES = FALSE
)

all_models[
  i = !(is_good),
  j = results := list(lme)
]

all_models[
  i = (is_good),
  j = `:=`(
    perf_check_model = lapply(
      X = lme,
      FUN = function(x) {
        devnull <- capture.output(
          out <- check_model(
            x = x,
            panel = FALSE,
            check = c("normality", "linearity", "qq"),
            verbose = FALSE
          )
        )
        out
      }
    ),
    perf_model_performance = lapply(lme, model_performance)
  )
]

saveRDS(
  object = all_models,
  file = file.path(output_directory, "summary_models.rds"),
  compress = FALSE
)
```

## Models results

```{r best-model}
best_models <- paste(
  "fixed = log(bmi) ~ gsp(age, knots = c(2, 8, 12), degree = rep(3, 4), smooth = rep(2, 3)),",
  "random = ~ gsp(age, knots = c(2, 8, 12), degree = rep(3, 4), smooth = rep(2, 3)) | ID"
)
best_dataset <- "visits_clean"
```

```{r get-best-model-results}
good_res_models <- readRDS(file.path(output_directory, "summary_models.rds"))[
  correlation_structure != "NULL" &
    models %in% best_models &
    dataset %in% best_dataset
]
```

```{r models-diagnostics-plot, dev.args = list(scaling = 0.50)}
wrap_plots(
  lapply(
    X = good_res_models[j = perf_check_model],
    FUN = function(.x) {
      if (is.null(.x)) {
        wrap_plots(rep(list(ggplot() + geom_blank()), 3), nrow = 1, tag_level = "new")
      } else {
        wrap_plots(plot(.x), nrow = 1, tag_level = "new")
      }
    }
  ),
  nrow = 2
) +
  plot_annotation(
    title = "Cubic Splines (Fixed Effects = Random Effects)",
    subtitle = "With A) sex = 0, and with B) sex = 1",
    tag_levels = c("A", "1")
  )
```

```{r model-predict-plot}
ggplot(
  data = good_res_models[
    j = rbindlist(list(
      (function(x) {
        data.table(age = seq(0, 16, 0.1))[
          j = bmi := if (inherits(x[[1]], "lme")) {
            exp(predict(x[[1]], .SD, level = 0))
          } else NA_real_,
          .SDcols = "age"
        ]
      })(lme)
    )),
    by = list(sex, dataset, models)
  ]
) +
  aes(
    x = age,
    y = bmi,
    colour = c("0" = "Male", "1" = "Female")[as.character(sex)]
  ) +
  geom_path() +
  labs(
    title = "Cubic Splines (Fixed Effects = Random Effects)",
    x = "AGE (years)",
    y = "BMI (kg/m\u00B2)",
    colour = NULL
  )
```

```{r models-diagnostics}
fmt_models <- all_models[
  j = list(
    sex = c("0" = "Male", "1" = "Female")[as.character(sex)],
    dataset = c(
      "visits_raw" = "Raw",
      "visits_clean" = "Clean",
      "visits_clean[between(age, 1, 18)]" = "Clean<br><i style='font-size: 6pt'>(1 ≤ age ≤ 18)</i>"
    )[dataset],
    correlation_structure,
    model = setNames(
      object = c(
        "<b style='color: #463480FF'>Linear Splines</b> - Random Linear Splines",
        "<b style='color: #463480FF'>Linear Splines</b> - Random Intercepts",
        "<b style='color: #21908CFF'>Cubic Splines</b> - Random Cubic Splines",
        "<b style='color: #21908CFF'>Cubic Splines</b> - Random Cubic Splines (1:3)",
        "<b style='color: #21908CFF'>Cubic Splines</b> - Random Quadratic Splines",
        "<b style='color: #21908CFF'>Cubic Splines</b> - Random Linear Splines",
        "<b style='color: #21908CFF'>Cubic Splines</b> - Random Intercepts",
        "<b style='color: #9AD93CFF'>Cubic Slope</b> - Random Cubic Slopes",
        "<b style='color: #9AD93CFF'>Cubic Slope</b> - Random Quadratic Slopes",
        "<b style='color: #9AD93CFF'>Cubic Slope</b> - Random Linear Slopes",
        "<b style='color: #9AD93CFF'>Cubic Slope</b> - Random Intercepts"
      ),
      nm = sprintf("fixed = %s, random = %s", fixed_effect, random_effect)
    )[models],
    results = sub(", block [0-9]+$| = [0-9.e-]+$", "", results),
    perf = perf_model_performance
  )
][
  j = c("model", "random") := tstrsplit(model, " - ")
][
  j = `:=`(
    y = paste0(
      model,
      fifelse(correlation_structure == "NULL", "<br>", " -- C-AR1<br>"),
      "<i style='font-size: 8pt'>", random, "</i>"
    )
  )
][
  j = `:=`(
    x = factor(x = dataset, levels = unique(dataset)),
    y = factor(x = y, levels = rev(unique(y))),
    fill = factor(
      x = sub(":.*", "", results),
      levels = c("Good", "Warning", "Error")
    )
  )
]

names_before_unlist <- names(fmt_models)

fwrite(
  x = fmt_models[
    j = unlist(perf, recursive = FALSE),
    by = list(
      sex,
      dataset = gsub("<b [^>]*>|</b>|<i [^>]*>|</i>", "", x),
      model = gsub("<br>", " -- ", gsub("<b [^>]*>|</b>|<i [^>]*>|</i>", "", y)),
      status = sub("TRUE", "Good", results)
    )
  ],
  file = file.path(output_directory, "models-performance.csv")
)
```

## Models status

```{r models-success-plot, fig.asp = 1}
ggplot(data = fmt_models) +
  aes(x = x, y = y, fill = fill) +
  facet_grid(cols = vars(sex)) +
  geom_tile(colour = "white") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  scale_fill_manual(
    values = c(
      "Good" = "#00b200",
      "Warning" = "#0000b2",
      "Error" = "#b20000"
    )
  ) +
  labs(x = NULL, y = NULL, fill = NULL) +
  theme(legend.position = "top")
```

```{r models-success-table}
if (fmt_models[fill != "Good", .N > 0]) {
  diag_tab <- gt(
    data = fmt_models[
      i = fill != "Good",
      j = list(
        Model = gsub("<i [^>]*>|</i>", "", gsub("<br>", " -- ", gsub("8pt", "0.75em", as.character(y)))),
        Dataset = gsub("<br>", " ", gsub("|<i [^>]*>|</i>", "", gsub("6pt", "0.75em", as.character(x)))),
        Result = results
      )
    ]
  ) %>%
    fmt(columns = everything(), fns = identity) %>%
    tab_header(title = "Models Diagnostics") %>%
    data_color(
      columns = "Result",
      colors = function(x) {
        c(
          "Good" = "#00b200",
          "Warning" = "#0000b2",
          "Error" = "#b20000"
        )[sub("TRUE", "Good", sub(":.*", "", x))]
      },
      apply_to = "text"
    ) %>%
    tab_options(
      table.font.size = "10pt",
      data_row.padding = "2pt"
    ) %>%
    opt_all_caps() %>%
    opt_row_striping()

  gtsave(diag_tab, file.path(output_directory, "models-diagnostics.html"))

  diag_tab
}
```

## Models performance

```{r models-performances-table}
perf_tab <- gt(
  data = fmt_models[
    i = fill %in% "Good",
    j = list(
      Sex = sex,
      Model = gsub("<i [^>]*>|</i>", "", gsub("<br>", " -- ", gsub("8pt", "0.75em", as.character(y)))),
      Dataset = gsub("<br>", " ", gsub("|<i [^>]*>|</i>", "", gsub("6pt", "0.75em", as.character(x)))),
      rbindlist(perf, fill = TRUE)
    )
  ],
  groupname_col = c("Sex", "Dataset")
) %>%
  fmt(columns = everything(), fns = identity) %>%
  fmt_number(
    columns = c("AIC", "BIC"),
    decimals = 1,
    drop_trailing_zeros = FALSE
  ) %>%
  fmt_number(
    columns = c("R2_conditional", "R2_marginal"),
    n_sigfig = 3,
    drop_trailing_zeros = FALSE
  ) %>%
  fmt_number(
    columns = c("RMSE", "Sigma"),
    n_sigfig = 3,
    drop_trailing_zeros = FALSE
  ) %>%
  fmt_number(
    columns = "ICC",
    rows = ICC >= 0.10,
    decimals = 3,
    n_sigfig = 3,
    drop_trailing_zeros = FALSE
  ) %>%
  fmt_scientific(
    columns = "ICC",
    rows = ICC < 0.10,
    decimals = 2,
    drop_trailing_zeros = FALSE
  ) %>%
  tab_header(title = "Models Performance") %>%
  data_color(
    columns = c("AIC", "BIC", "R2_conditional", "R2_marginal", "ICC", "RMSE", "Sigma"),
    colors = function(x) {
      sx <- rescale(x, to = c(0, 1))
      col_numeric(palette = "plasma", domain = c(0, 1))(sx)
    }
  ) %>%
  tab_options(
    table.font.size = "10pt",
    data_row.padding = "2pt"
  ) %>%
  opt_all_caps() %>%
  opt_row_striping()

gtsave(perf_tab, file.path(output_directory, "models-performance.html"))

perf_tab
```