# =============================================================================
# Supportive Periodontal Therapy Clinical Examination Data
# Script 3 of 3: Reproduction of Figure 2 (PPD stability thresholds)
#
# Reproduces the empirically determined thresholds of no change of residual
# probing depth reported in Figure 2 of
#   Ramseier CA, Nydegger M, Walter C, Fischer G, Sculean A, Lang NP, Salvi GE.
#   J Clin Periodontol. 2019;46(2):218-230. doi:10.1111/jcpe.13041
#
# These twenty thresholds are the empirical basis of the SPT interval algorithm
# implemented in 01_compute_spt_algorithm.R.
#
# Verified: all 20 published thresholds are reproduced exactly from the
# published data files.
#
# Requires only base R (graphics optional).
# =============================================================================

## ---- 1. Data ----------------------------------------------------------------

# The data file is expected in the working directory. If R reports that it
# cannot be opened, point the working directory at the folder holding the data:
#   setwd("path/to/the/data")

file <- "02_supportive_periodontal_therapy_data.csv"
if (!file.exists(file)) {
  stop(sprintf(paste0("'%s' not found in the working directory (%s).\n",
                      "  Set the working directory to the folder holding the ",
                      "data files, e.g. setwd(\"path/to/the/data\")."),
               file, getwd()), call. = FALSE)
}

spt <- read.csv(file, stringsAsFactors = FALSE)
spt <- spt[order(spt$pat_id, spt$spt_id), ]

spt$n_sites_spt <- spt$n_teeth_spt * 6

# Cumulative percentages of affected sites at each visit
spt$pct4 <- (spt$n_4mm_spt + spt$n_5mm_spt + spt$n_6mm_spt + spt$n_from7mm_spt) * 100 / spt$n_sites_spt
spt$pct5 <- (                spt$n_5mm_spt + spt$n_6mm_spt + spt$n_from7mm_spt) * 100 / spt$n_sites_spt
spt$pct6 <- (                                spt$n_6mm_spt + spt$n_from7mm_spt) * 100 / spt$n_sites_spt
spt$pct7 <-                                                  spt$n_from7mm_spt  * 100 / spt$n_sites_spt


## ---- 2. Change relative to the preceding visit ------------------------------
# NOTE on the grouping variable.
# The figure caption states that visits are grouped by the residual PPD category
# "recorded at the previous SPT visit". The published values are reproduced only
# when visits are grouped by the percentage recorded at the RESPECTIVE (current)
# visit, with the change measured backwards to the preceding visit. The code
# below follows the computation that reproduces the published figure.

lagv <- function(x, id) ave(x, id, FUN = function(v) c(NA, head(v, -1)))

for (k in c("4", "5", "6", "7")) {
  cur <- spt[[paste0("pct", k)]]
  spt[[paste0("chg", k)]] <- cur - lagv(cur, spt$pat_id)
}


## ---- 3. Interval categories -------------------------------------------------
# Real time between visits, categorised into 3, 4, 6, 9 and 12+ months.
# true_interval_days is 0 at each patient's first visit (no predecessor).

spt$ivl_months <- round(spt$true_interval_days / 30, 0)

spt$ivl_cat <- cut(spt$ivl_months,
                   breaks = c(0, 3, 5, 8, 11, Inf),
                   labels = c("3", "4", "6", "9", "12+"),
                   right = TRUE)
spt$ivl_cat[spt$true_interval_days == 0] <- NA


## ---- 4. Percentage categories per panel -------------------------------------
# Upper bounds of the bins as printed on the x axis of Figure 2.

bins <- list("4" = c(10, 20, 30, 40, 100),
             "5" = c(2, 4, 6, 8, 10, 20, 30, 100),
             "6" = c(1, 2, 3, 4, 10, 100),
             "7" = c(1, 2, 3, 4, 10, 100))

ivls <- c("3", "4", "6", "9", "12+")


## ---- 5. Cell means and threshold extraction ---------------------------------
# The threshold is the largest bin whose mean change is still <= 0, i.e. the
# highest percentage of affected sites at which no increase of residual PPD is
# observed for that interval length.

panel_table <- function(k, exclude_zero = FALSE) {
  edges <- bins[[k]]
  cur <- spt[[paste0("pct", k)]]
  chg <- spt[[paste0("chg", k)]]
  keep <- !is.na(chg) & !is.na(spt$ivl_cat)
  if (exclude_zero) keep <- keep & cur > 0

  m <- matrix(NA_real_, nrow = length(edges), ncol = length(ivls),
              dimnames = list(paste0("<=", edges), ivls))
  n <- m
  for (j in seq_along(edges)) {
    lo <- if (j == 1) -1 else edges[j - 1]
    for (i in seq_along(ivls)) {
      s <- keep & cur > lo & cur <= edges[j] & spt$ivl_cat == ivls[i]
      m[j, i] <- mean(chg[s]); n[j, i] <- sum(s)
    }
  }
  list(mean = m, n = n, edges = edges)
}

threshold <- function(tab) {
  out <- setNames(rep(NA_real_, length(ivls)), ivls)
  for (i in seq_along(ivls)) {
    last_safe <- NA_real_
    for (j in seq_along(tab$edges)) {
      v <- tab$mean[j, i]
      if (is.na(v)) next
      if (v <= 0) last_safe <- tab$edges[j] else break
    }
    out[i] <- last_safe
  }
  out
}


## ---- 6. Verification against the published thresholds -----------------------

published <- rbind("4" = c(30, 20, 20, 10, 10),
                   "5" = c(20, 10,  6,  4,  2),
                   "6" = c( 4,  3,  2,  1,  1),
                   "7" = c( 2,  1,  1,  1,  1))
colnames(published) <- ivls

cat("Thresholds of no change of residual PPD (percentage of affected sites)\n")
cat("=====================================================================\n\n")
for (k in c("4", "5", "6", "7")) {
  tab <- panel_table(k)
  rec <- threshold(tab)
  cat(sprintf("PPD >= %s mm\n", k))
  cat("  published   :", sprintf("%5s", published[k, ]), "\n")
  cat("  reconstructed:", sprintf("%5s", rec), "\n")
  cat("  match       :", all(rec == published[k, ]), "\n\n")
}


## ---- 7. Sensitivity: excluding visits with zero affected sites --------------
# Visits with no residual pocket in a given category contribute a change of
# approximately zero and dominate the lowest bin, in particular for PPD >= 6 mm
# (50.6% of visits) and PPD >= 7 mm (73.4%). Re-deriving the thresholds without
# them tests whether the algorithm depends on those zeros.

cat("\nSensitivity: thresholds with and without visits at 0% affected sites\n")
cat("====================================================================\n\n")
for (k in c("4", "5", "6", "7")) {
  cat(sprintf("PPD >= %s mm\n", k))
  cat("  with zeros   :", sprintf("%5s", threshold(panel_table(k, FALSE))), "\n")
  cat("  without zeros:", sprintf("%5s", threshold(panel_table(k, TRUE))),  "\n\n")
}


## ---- 8. Cell means, for inspection or plotting ------------------------------

for (k in c("4", "5", "6", "7")) {
  tab <- panel_table(k)
  cat(sprintf("\nMean change, PPD >= %s mm (n in brackets)\n", k))
  disp <- matrix(sprintf("%+.2f (%d)", tab$mean, tab$n),
                 nrow = nrow(tab$mean), dimnames = dimnames(tab$mean))
  print(disp, quote = FALSE)
  write.csv(tab$mean, sprintf("figure2_panel_ppd%s.csv", k))
}

sessionInfo()