# =============================================================================
# Supportive Periodontal Therapy Clinical Examination Data
# Script 1 of 2: Residual-PPD-based algorithm for computing SPT intervals
#
# Reproduces every derived variable in 03_supplementary_data.csv from the
# primary data in 01_initial_periodontal_therapy_data.csv and
# 02_supportive_periodontal_therapy_data.csv.
#
# The algorithm implements the empirically determined PPD stability thresholds
# reported in Ramseier et al. 2019 (J Clin Periodontol 46:218-230), Figure 2.
#
# Verified: all 11,842 rows reproduce the published values exactly.
#
# Requires only base R. No external packages.
# =============================================================================
## ---- 1. Read the data -------------------------------------------------------
# The data files are expected in the working directory. If R reports that a file
# cannot be opened, point the working directory at the folder holding the data:
# setwd("path/to/the/data")
read_data <- function(file) {
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)
}
read.csv(file, stringsAsFactors = FALSE)
}
apt <- read_data("01_initial_periodontal_therapy_data.csv")
spt <- read_data("02_supportive_periodontal_therapy_data.csv")
spt <- spt[order(spt$pat_id, spt$spt_id), ]
## ---- 2. Cumulative site counts at each SPT visit ----------------------------
# PPD categories are recorded as counts of sites at exactly 4, 5, 6 mm and at
# >= 7 mm. The algorithm operates on cumulative counts (>= 4, >= 5, >= 6 mm).
spt$n_from4mm_spt <- spt$n_4mm_spt + spt$n_5mm_spt + spt$n_6mm_spt + spt$n_from7mm_spt
spt$n_from5mm_spt <- spt$n_5mm_spt + spt$n_6mm_spt + spt$n_from7mm_spt
spt$n_from6mm_spt <- spt$n_6mm_spt + spt$n_from7mm_spt
## ---- 3. Site percentages ----------------------------------------------------
# PPD was recorded at six sites per tooth, so the denominator is n_teeth * 6.
spt$n_sites_spt <- spt$n_teeth_spt * 6
spt$percent_4mm_spt <- spt$n_4mm_spt * 100 / spt$n_sites_spt
spt$percent_5mm_spt <- spt$n_5mm_spt * 100 / spt$n_sites_spt
spt$percent_6mm_spt <- spt$n_6mm_spt * 100 / spt$n_sites_spt
spt$percent_from7mm_spt <- spt$n_from7mm_spt * 100 / spt$n_sites_spt
spt$percent_from4mm_spt <- spt$n_from4mm_spt * 100 / spt$n_sites_spt
spt$percent_from5mm_spt <- spt$n_from5mm_spt * 100 / spt$n_sites_spt
spt$percent_from6mm_spt <- spt$n_from6mm_spt * 100 / spt$n_sites_spt
## ---- 4. PPD stability thresholds -------------------------------------------
# For each residual PPD category (>= 4, >= 5, >= 6 mm) and each candidate SPT
# interval (3, 4, 6, 12 months), a threshold percentage of affected sites was
# determined below which no increase of residual PPD was observed between two
# consecutive visits. A candidate interval is admissible only when the observed
# percentage lies at or below its threshold.
#
# PPD >= 4 mm : 3 months if > 20% ; 4 and 6 months if <= 20% ; 12 months if <= 10%
# PPD >= 5 mm : 3 months if > 10% ; 4 months if <= 10% ; 6 months if <= 6% ; 12 months if <= 2%
# PPD >= 6 mm : 3 months if > 3% ; 4 months if <= 3% ; 6 months if <= 2% ; 12 months if <= 1%
#
# PPD >= 7 mm is deliberately NOT used: sites of >= 8 mm were not recorded at
# the MSDH, so >= 7 mm did not qualify as a threshold variable
# (Ramseier et al. 2019, section 3.2).
#
# NA encodes "this interval is not admissible for this PPD category".
thr <- function(cond, months) ifelse(cond, months, NA_real_)
spt$from4mm_3months <- thr(spt$percent_from4mm_spt > 20, 3)
spt$from4mm_4months <- thr(spt$percent_from4mm_spt <= 20, 4)
spt$from4mm_6months <- thr(spt$percent_from4mm_spt <= 20, 6)
spt$from4mm_12months <- thr(spt$percent_from4mm_spt <= 10, 12)
spt$from5mm_3months <- thr(spt$percent_from5mm_spt > 10, 3)
spt$from5mm_4months <- thr(spt$percent_from5mm_spt <= 10, 4)
spt$from5mm_6months <- thr(spt$percent_from5mm_spt <= 6, 6)
spt$from5mm_12months <- thr(spt$percent_from5mm_spt <= 2, 12)
spt$from6mm_3months <- thr(spt$percent_from6mm_spt > 3, 3)
spt$from6mm_4months <- thr(spt$percent_from6mm_spt <= 3, 4)
spt$from6mm_6months <- thr(spt$percent_from6mm_spt <= 2, 6)
spt$from6mm_12months <- thr(spt$percent_from6mm_spt <= 1, 12)
## ---- 5. Longest admissible interval per PPD category ------------------------
# Within each PPD category, take the longest admissible interval.
# NA (not admissible) counts as 0, matching the original spreadsheet.
# NA here has two distinct meanings and they must not be conflated: an interval
# that is not admissible for an observed percentage, and a percentage that was
# never measured. The first counts as 0 and competes in the maximum; the second
# makes the whole result unknown. Visits flagged by `imputed` fall in the second
# case, and every value derived from them stays missing.
rowmax0 <- function(pct, ...) {
m <- cbind(...)
m[is.na(m)] <- 0
out <- do.call(pmax, as.data.frame(m))
out[is.na(pct)] <- NA_real_
out
}
spt$max_from4mm <- rowmax0(spt$percent_from4mm_spt,
spt$from4mm_3months, spt$from4mm_4months,
spt$from4mm_6months, spt$from4mm_12months)
spt$max_from5mm <- rowmax0(spt$percent_from5mm_spt,
spt$from5mm_3months, spt$from5mm_4months,
spt$from5mm_6months, spt$from5mm_12months)
spt$max_from6mm <- rowmax0(spt$percent_from6mm_spt,
spt$from6mm_3months, spt$from6mm_4months,
spt$from6mm_6months, spt$from6mm_12months)
## ---- 6. Computed SPT interval ----------------------------------------------
# The binding constraint is the shortest of the three category-wise intervals.
# A category contributing 0 (no admissible interval) is ignored.
na999 <- function(x) ifelse(x != 0, x, 999)
spt$algorithm_based_interval_months <-
pmin(na999(spt$max_from4mm), na999(spt$max_from5mm), na999(spt$max_from6mm))
## ---- 7. Adherence to the computed interval ----------------------------------
# true_interval_days is 0 at each patient's first SPT visit (no preceding
# visit exists). Months are obtained by rounding days / 30, as in the original
# analysis.
spt$true_interval_months <- round(spt$true_interval_days / 30, 0)
spt$interval_dif <- ifelse(spt$true_interval_days == 0, 0,
spt$true_interval_months - spt$algorithm_based_interval_months)
# No computed interval means no comparison.
spt$interval_dif[is.na(spt$algorithm_based_interval_months)] <- NA_real_
spt$interval_shorter <- ifelse(spt$interval_dif < 0, 1, NA_integer_)
spt$interval_exact <- ifelse(spt$interval_dif == 0, 1, NA_integer_)
spt$interval_longer <- ifelse(spt$interval_dif > 0, 1, NA_integer_)
# NOTE: interval_shorter / interval_exact / interval_longer compare the actual
# interval against the ALGORITHM-BASED interval, not against the interval
# assigned by the dental hygienist (assigned_spt_interval_months).
## ---- 8. Runs of consecutive visits in the same adherence category -----------
run_length <- function(id, flag) {
flag <- ifelse(is.na(flag), 0L, as.integer(flag))
out <- integer(length(flag)); run <- 0L
for (i in seq_along(flag)) {
if (i > 1 && id[i] != id[i - 1]) run <- 0L
run <- if (flag[i] == 1L) run + 1L else 0L
out[i] <- run
}
out
}
spt$consecutive_shorter <- run_length(spt$pat_id, spt$interval_shorter)
spt$consecutive_exact <- run_length(spt$pat_id, spt$interval_exact)
spt$consecutive_longer <- run_length(spt$pat_id, spt$interval_longer)
## ---- 9. Verification against the published supplementary file ---------------
if (file.exists("03_supplementary_data.csv")) {
ref <- read.csv("03_supplementary_data.csv", stringsAsFactors = FALSE)
ref <- ref[order(ref$pat_id, ref$spt_id), ]
check <- c("n_from4mm_spt", "n_from5mm_spt", "n_from6mm_spt",
"percent_4mm_spt", "percent_5mm_spt", "percent_6mm_spt",
"percent_from7mm_spt", "percent_from4mm_spt",
"percent_from5mm_spt", "percent_from6mm_spt",
"from4mm_3months", "from4mm_4months", "from4mm_6months", "from4mm_12months",
"from5mm_3months", "from5mm_4months", "from5mm_6months", "from5mm_12months",
"from6mm_3months", "from6mm_4months", "from6mm_6months", "from6mm_12months",
"max_from4mm", "max_from5mm", "max_from6mm",
"algorithm_based_interval_months",
"interval_shorter", "interval_exact", "interval_longer",
"consecutive_shorter", "consecutive_exact", "consecutive_longer")
cat("\nVerification against 03_supplementary_data.csv\n")
cat("---------------------------------------------\n")
n_fail <- 0
for (v in check) {
a <- spt[[v]]; b <- ref[[v]]
same <- (is.na(a) & is.na(b)) |
(!is.na(a) & !is.na(b) & abs(a - b) < 1e-9)
ok <- sum(same)
if (ok < nrow(ref)) n_fail <- n_fail + 1
cat(sprintf(" %-32s %6d / %d identical%s\n", v, ok, nrow(ref),
if (ok < nrow(ref)) " <-- MISMATCH" else ""))
}
cat("\n")
if (n_fail == 0) {
cat("All derived variables reproduce exactly.\n")
} else {
cat(sprintf("%d variable(s) did not reproduce.\n", n_fail))
}
}
sessionInfo()