# =============================================================================
# Supportive Periodontal Therapy Clinical Examination Data
# Linear mixed-effects model
#
# Reconstruction of the model reported in Table 2 of
# Ramseier CA, Nydegger M, Walter C, Fischer G, Sculean A, Lang NP, Salvi GE.
# Time between recall visits and residual probing depths predict long-term
# stability in patients enrolled in supportive periodontal therapy.
# J Clin Periodontol. 2019;46(2):218-230. doi:10.1111/jcpe.13041
#
# A Python version of this script, 02_linear_mixed_effects_model.py, produces
# the same numbers to four decimal places.
#
#
# PROVENANCE
# ----------
# The source code written for the 2019 publication is no longer available. This
# script is a RECONSTRUCTION, derived from the Materials and Methods (section
# 2.7), from Table 2, from the statistician's working documents, and from the
# data file that was sent to the statistician in November 2017. It runs on the
# published data files alone.
#
#
# ANALYSIS POPULATION
# -------------------
# Table 2 rests on the five-year subsample of 445 patients, not on all 883.
# The paper does not say so, and this had to be established from the data.
#
# The severity of periodontal disease at the initial examination is one of the
# covariates, and it was recorded only for the 445 patients who were
# systemically healthy and attended supportive periodontal therapy for at least
# five years. A complete-case fit is therefore confined to that subsample. The
# scale of the published estimates confirms it: the outcome is rank
# transformed, so every coefficient is proportional to the number of
# observations entering the fit. Fitted on the subsample, the coefficient for
# the residual probing depth at the visit is 146.78 against the published
# 145.18; fitted on all patients it is 186.54, twenty-eight per cent too large,
# and the F values are inflated by a third.
#
# The model is therefore fitted twice below. Section 8 gives the subsample fit,
# which corresponds to Table 2. Section 9 repeats it on all 883 patients, using
# the disease severity derived for the remaining 438 from their baseline probing
# depth profile; that fit is not the published one but shows how the findings
# behave in the full cohort.
#
#
# HOW CLOSE THE RECONSTRUCTION COMES
# ----------------------------------
# Subsample fit, 8,221 visits from 445 patients, against Table 2:
#
# Covariate Published Reconstructed
# ---------------------------------------------------------------
# PPD% at respective SPT visit 145.1780 146.7790
# Time between two SPT visits 0.4302 0.3824
# BOP% at respective SPT visit 0.2982 -0.6998
# PPD% following active therapy -46.6178 -38.9592
# PPD% at initial examination -5.5292 -9.6748
# Number of teeth at baseline 24.9403 20.1734
# Gender 152.3959 171.0283
# Age at initial examination -0.5756 2.9788
#
# Standard errors agree to within eight per cent for every covariate except the
# two baseline probing depth terms. Eleven of the twelve estimates lie within
# two standard errors of the published value; the exception is severe
# periodontal disease.
#
# Sequential F values, against the column printed in Table 2:
#
# PPD% at respective SPT visit 1460.32 -> 1494.43 (1.02)
# BOP% at respective SPT visit 283.15 -> 311.52 (1.10)
# Time between two SPT visits 37.59 -> 33.79 (0.90)
# PPD% following active therapy 3.70 -> 2.61 (0.71)
# Periodontal disease 0.63 -> 0.77 (1.21)
# Gender 0.22 -> 0.25 (1.12)
# Number of teeth at baseline 0.07 -> 0.08 (1.23)
# PPD% at initial examination 0.19 -> 0.16 (0.88)
# Smoking status 0.01 -> 0.01 (1.10)
# Age at initial examination 0.10 -> 0.22 (2.28)
#
# The two categorical rows are compared against the statistician's working
# document rather than against the printed table, which carries a decimal error
# in both of them; see the note at the end of this header.
#
#
# THE F COLUMN OF TABLE 2 IS SEQUENTIAL
# -------------------------------------
# This resolves an apparent contradiction in the published table. Bleeding on
# probing is reported with a 95% confidence interval spanning zero,
# (-3.8351, 4.3178), alongside F = 283.15 and p < 0.0001.
#
# Both figures are correct and they answer different questions. Under a
# sequential (type I) decomposition, in which bleeding on probing enters before
# the residual probing depth measured at the same visit, it carries the variance
# the two share and its F is large. Adjusted for that probing depth, as the
# confidence interval is, its independent contribution is negligible. The two
# variables correlate at r = 0.39.
#
# The statistician noticed this himself. In a message of 3 December 2017 he
# wrote that the confidence interval contained zero, which would mean a
# non-significant t test and therefore contradicted the F test he was reporting,
# and advised not to raise the point unless a reader asked.
#
# The script therefore reports BOTH decompositions. The type I tests correspond
# to the published F column. The type III tests correspond to the published
# confidence intervals, and they are the ones to build on: the type I column
# depends on the order in which terms were entered and is not a test of an
# independent contribution.
#
#
# BLEEDING ON PROBING WAS RECORDED AT FOUR SITES PER TOOTH
# --------------------------------------------------------
# Probing depth was recorded at six sites per tooth, bleeding on probing at
# four. The denominator for the BOP percentage is therefore n_teeth_spt * 4.
# This reproduces the BOP percentages in the file sent to the statistician
# exactly, for all 11,842 visits, and the maximum of n_bop_pos_spt / n_teeth_spt
# in the data is exactly 4.
#
#
# TWO DECIMAL ERRORS IN THE PRINTED TABLE
# ---------------------------------------
# The statistician's working document of 3 December 2017 gives F = 0.6330 for
# periodontal disease and F = 0.0123 for smoking status. Table 2 prints 0.06 and
# 0.12. The reconstruction returns 0.77 and 0.01, which follows the working
# document in both rows; the printed values are decimal slips in opposite
# directions. Both rows are far from significance in every version, so nothing
# in the paper depends on them.
#
# The p values in those same two rows of the working document do not match their
# own F values either: 0.9387 corresponds to F = 0.0633 and 0.9978 to F = 0.0022.
# The p values of all eight single-degree-of-freedom rows are exactly consistent
# with their F values. Where the two categorical p values came from cannot be
# determined from the surviving material.
#
#
# LABELLING OF THE GENDER ROW
# ---------------------------
# Table 2 labels the gender row "Gender (male)" and reports 152.3959. In the
# data file used for the analysis the gender variable was coded 1 for female,
# and the reconstruction returns +171.03 for female, -171.03 for male. The
# published estimate is therefore the effect of being female under a label
# reading male. Both analyses agree that women showed the larger increase.
#
#
# WHAT WAS NOT DONE
# -----------------
# No search was made for a specification that reproduces the published F column
# more closely. Tuning a model until it matches a published table, without the
# original code to justify each choice, produces agreement without evidence.
#
# The algorithm derived from these data is a separate matter and is fully
# reproducible: see 01_compute_spt_algorithm.R, which reproduces all 11,842 rows
# of the derived data exactly, and 03_reproduce_figure_2.R, which reproduces all
# twenty published stability thresholds exactly.
#
# Requires: lme4, lmerTest
# =============================================================================
library(lme4)
library(lmerTest)
## ---- 1. Read the published data files ---------------------------------------
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. Patient-level covariates from the initial examination ---------------
apt$n_from4mm_baseline <- apt$n_4mm_baseline + apt$n_5mm_baseline +
apt$n_6mm_baseline + apt$n_from7mm_baseline
apt$n_from4mm_follow_up <- apt$n_4mm_follow_up + apt$n_5mm_follow_up +
apt$n_6mm_follow_up + apt$n_from7mm_follow_up
# "PPD% at initial examination" and "PPD% following APT" are the percentages of
# sites with PPD >= 4 mm. Six sites per tooth were probed.
apt$percent_from4mm_baseline <- apt$n_from4mm_baseline * 100 / (apt$n_teeth_baseline * 6)
apt$percent_from4mm_follow_up <- apt$n_from4mm_follow_up * 100 / (apt$n_teeth_follow_up * 6)
## ---- 3. Visit-level covariates ----------------------------------------------
spt$percent_from4mm_spt <- (spt$n_4mm_spt + spt$n_5mm_spt + spt$n_6mm_spt +
spt$n_from7mm_spt) * 100 / (spt$n_teeth_spt * 6)
# Bleeding on probing was recorded at FOUR sites per tooth, probing depth at six.
spt$bop_percent <- spt$n_bop_pos_spt * 100 / (spt$n_teeth_spt * 4)
## ---- 4. Outcome: change in residual PPD >= 4 mm between consecutive visits --
# Table 2 models "residual periodontal probing depth changes between two
# consecutive visits". The first visit of each patient has no predecessor and
# therefore no outcome value.
spt$ppd_change <- spt$percent_from4mm_spt -
ave(spt$percent_from4mm_spt, spt$pat_id, FUN = function(x) c(NA, head(x, -1)))
## ---- 5. Assemble the analysis data set --------------------------------------
dat <- merge(spt,
apt[, c("pat_id", "gender", "age_baseline", "smoking",
"n_teeth_baseline", "diagnosis", "subsample_5years",
"percent_from4mm_baseline", "percent_from4mm_follow_up")],
by = "pat_id", all.x = TRUE)
# Reference categories. Gender is referenced to male so that the estimate is the
# effect of being female, which is the quantity Table 2 reports; see the header.
dat$gender <- relevel(factor(dat$gender), ref = "male")
dat$smoking <- relevel(factor(dat$smoking), ref = "non-smoker")
dat$diagnosis <- relevel(factor(dat$diagnosis), ref = "1")
model_vars <- c("pat_id", "ppd_change", "gender", "age_baseline", "smoking",
"n_teeth_baseline", "diagnosis", "true_interval_days",
"percent_from4mm_baseline", "percent_from4mm_follow_up",
"bop_percent", "percent_from4mm_spt", "subsample_5years")
dat <- dat[complete.cases(dat[, model_vars]), model_vars]
dat <- dat[dat$true_interval_days > 0, ]
## ---- 6. Published values, for comparison ------------------------------------
# Effect sizes and F values as printed in Table 2, except for the two
# categorical F values, which are taken from the statistician's working document
# of 3 December 2017; see the header.
published <- data.frame(
term = c("gender", "age_baseline", "smoking", "n_teeth_baseline",
"diagnosis", "true_interval_days", "percent_from4mm_baseline",
"percent_from4mm_follow_up", "bop_percent", "percent_from4mm_spt"),
F_value = c(0.2222, 0.0961, 0.0123, 0.0674, 0.6330, 37.5900,
0.1878, 3.7042, 283.1533, 1460.3165),
stringsAsFactors = FALSE
)
## ---- 7. Fit -----------------------------------------------------------------
# The outcome is rank transformed: the residuals of the untransformed model were
# not normally distributed (section 2.7). Estimates are therefore on the rank
# scale and are not interpretable in percentage points; their sign and relative
# magnitude carry the interpretation, and their absolute size depends on the
# number of observations.
#
# Term order matters for the sequential tests. Bleeding on probing is entered
# before the residual probing depth at the same visit, which is the order that
# reproduces the published F column.
fit_model <- function(d, label, file) {
d$ppd_change_rank <- rank(d$ppd_change)
cat("\n", strrep("=", 78), "\n", label, "\n", strrep("=", 78), "\n", sep = "")
cat(sprintf("%d SPT visits from %d patients\n",
nrow(d), length(unique(d$pat_id))))
fit <- lmer(
ppd_change_rank ~ gender + age_baseline + smoking + n_teeth_baseline +
diagnosis + true_interval_days +
percent_from4mm_baseline + percent_from4mm_follow_up +
bop_percent + percent_from4mm_spt +
(1 | pat_id),
data = d, REML = TRUE
)
cat("\n--- Estimates and 95% confidence intervals ---\n")
cat("These correspond to the effect sizes and intervals of Table 2.\n\n")
ct <- summary(fit)$coefficients
ci <- confint(fit, method = "Wald")
print(round(cbind(Estimate = ct[, "Estimate"],
Std.Error = ct[, "Std. Error"],
ci[rownames(ct), , drop = FALSE]), 4))
a1 <- anova(fit, type = 1)
a3 <- anova(fit, type = 3)
m <- match(published$term, rownames(a1))
cat("\n--- Type I, sequential: corresponds to the F column of Table 2 ---\n")
cat("Depends on the order in which terms were entered. Not a test of an\n")
cat("independent contribution.\n\n")
print(data.frame(term = published$term,
published_F = published$F_value,
reconstructed = round(a1[m, "F value"], 4),
ratio = round(a1[m, "F value"] / published$F_value, 3),
row.names = NULL))
cat("\n--- Type III, marginal: use these ---\n")
cat("Each term adjusted for all others. These correspond to the confidence\n")
cat("intervals above.\n\n")
print(round(a3, 4))
out <- data.frame(variable = rownames(ct),
estimate = ct[, "Estimate"],
std_error = ct[, "Std. Error"],
df = ct[, "df"],
t_value = ct[, "t value"],
p_value = ct[, "Pr(>|t|)"],
row.names = NULL)
write.csv(out, file, row.names = FALSE)
cat(sprintf("\nWritten: %s\n", file))
invisible(fit)
}
## ---- 8. The model of Table 2 -------------------------------------------------
fit_model(dat[dat$subsample_5years == 1, ],
"Table 2: five-year subsample, disease severity as recorded",
"mixed_model_results_subsample.csv")
## ---- 9. The same model in the full cohort ------------------------------------
# Not the published fit. Disease severity is derived from the baseline probing
# depth profile for the 438 patients for whom it was not recorded; that
# derivation agrees with every one of the 445 recorded values and reproduces the
# totals of Table 1 exactly. Coefficients are on a larger rank scale and are not
# comparable in size with those above.
fit_model(dat,
"All patients: not the published fit, shown for comparison",
"mixed_model_results_all_patients.csv")
sessionInfo()