# =============================================================================
# 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()