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