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