""" Supportive Periodontal Therapy Clinical Examination Data Linear mixed-effects model Python equivalent of 02_linear_mixed_effects_model.R. See the header of that script for the provenance of the model, the analysis population, the sequential character of the F column of Table 2, and the two decimal errors in the printed table. Everything said there applies here. Both implementations were written against the same specification and agree to four decimal places in the estimates, standard errors and F values. The mixed model is fitted here without a modelling library. A random-intercept model has a closed-form profiled REML criterion in the single variance ratio sigma_u^2 / sigma_e^2, which is minimised by golden-section search; the generalised least squares solution and its covariance follow directly. The sequential (type I) tests use the contrast matrices obtained from the normalised triangular factor of X'X, which is the construction lmerTest uses, so the F values are identical to those of anova(fit, type = 1) in R. The marginal (type III) tests are the joint Wald tests of each term, which for a model without interactions is what anova(fit, type = 3) returns. Degrees of freedom are the one thing not reproduced here: lmerTest obtains the denominator degrees of freedom by Satterthwaite's method, which is not implemented in this script, so no p values are printed. The F values themselves do not depend on it. Use the R script when p values are needed. Requires: pandas, numpy. Run from the folder holding the data files: python 02_linear_mixed_effects_model.py """ import os import sys import numpy as np import pandas as pd if len(sys.argv) > 1: os.chdir(sys.argv[1]) # ---- 1. A random-intercept REML fitter -------------------------------------- def fit_random_intercept(X, y, groups): """REML fit of y = X b + u_group + e. Returns (b, cov_b, sigma2_e, ratio).""" _, gidx = np.unique(groups, return_inverse=True) n_g = np.bincount(gidx) N, p = X.shape # Group sums, which are all the sufficient statistics a random intercept needs. Sx = np.zeros((len(n_g), p)) Sy = np.zeros(len(n_g)) np.add.at(Sx, gidx, X) np.add.at(Sy, gidx, y) XtX, Xty, yty = X.T @ X, X.T @ y, y @ y def criterion(ratio): # V = sigma2_e * (I + ratio * ZZ'), block diagonal, so its inverse is # I - c_g * J_g within each group. c = ratio / (1.0 + n_g * ratio) A = XtX - (Sx.T * c) @ Sx b = Xty - Sx.T @ (c * Sy) beta = np.linalg.solve(A, b) rss = yty - np.sum(c * Sy ** 2) - b @ beta _, logdet_A = np.linalg.slogdet(A) deviance = (N - p) * np.log(rss) + np.sum(np.log(1.0 + n_g * ratio)) + logdet_A return deviance, beta, rss, A # Golden-section search on log10 of the variance ratio. phi = (np.sqrt(5.0) - 1.0) / 2.0 lo, hi = -6.0, 4.0 x1, x2 = hi - phi * (hi - lo), lo + phi * (hi - lo) f1, f2 = criterion(10 ** x1)[0], criterion(10 ** x2)[0] while hi - lo > 1e-10: if f1 < f2: hi, x2, f2 = x2, x1, f1 x1 = hi - phi * (hi - lo) f1 = criterion(10 ** x1)[0] else: lo, x1, f1 = x1, x2, f2 x2 = lo + phi * (hi - lo) f2 = criterion(10 ** x2)[0] ratio = 10 ** ((lo + hi) / 2) _, beta, rss, A = criterion(ratio) sigma2_e = rss / (N - p) return beta, sigma2_e * np.linalg.inv(A), sigma2_e, ratio def sequential_contrasts(X): """Normalised triangular factor of X'X: the type I contrast matrix.""" R = np.linalg.cholesky(X.T @ X).T return R / np.diag(R)[:, None] def wald_F(beta, cov, L): Lb = L @ beta return float(Lb @ np.linalg.solve(L @ cov @ L.T, Lb) / L.shape[0]) def average_rank(x): """R's rank(), which averages ties.""" x = np.asarray(x, dtype=float) order = np.argsort(x, kind="mergesort") ordered = x[order] out = np.empty(len(x)) i = 0 while i < len(x): j = i while j + 1 < len(x) and ordered[j + 1] == ordered[i]: j += 1 out[order[i:j + 1]] = (i + j) / 2.0 + 1 i = j + 1 return out # ---- 2. Read the published data files --------------------------------------- def read_data(filename): if not os.path.exists(filename): sys.exit(f"'{filename}' not found in the working directory ({os.getcwd()}).\n" " Run this script from the folder holding the data files, or " "pass that folder as the first argument.") return pd.read_csv(filename) apt = read_data("01_initial_periodontal_therapy_data.csv") spt = read_data("02_supportive_periodontal_therapy_data.csv") spt = spt.sort_values(["pat_id", "spt_id"]).reset_index(drop=True) # ---- 3. Covariates ----------------------------------------------------------- # Probing depth was recorded at six sites per tooth, bleeding on probing at four. apt["percent_from4mm_baseline"] = ( apt[["n_4mm_baseline", "n_5mm_baseline", "n_6mm_baseline", "n_from7mm_baseline"]].sum(axis=1) * 100 / (apt["n_teeth_baseline"] * 6)) apt["percent_from4mm_follow_up"] = ( apt[["n_4mm_follow_up", "n_5mm_follow_up", "n_6mm_follow_up", "n_from7mm_follow_up"]].sum(axis=1) * 100 / (apt["n_teeth_follow_up"] * 6)) spt["percent_from4mm_spt"] = ( spt[["n_4mm_spt", "n_5mm_spt", "n_6mm_spt", "n_from7mm_spt"]] .sum(axis=1, min_count=4) * 100 / (spt["n_teeth_spt"] * 6)) spt["bop_percent"] = spt["n_bop_pos_spt"] * 100 / (spt["n_teeth_spt"] * 4) # ---- 4. Outcome -------------------------------------------------------------- spt["ppd_change"] = (spt["percent_from4mm_spt"] - spt.groupby("pat_id")["percent_from4mm_spt"].shift(1)) dat = spt.merge( apt[["pat_id", "gender", "age_baseline", "smoking", "n_teeth_baseline", "diagnosis", "subsample_5years", "percent_from4mm_baseline", "percent_from4mm_follow_up"]], on="pat_id", how="left") MODEL_VARS = ["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[MODEL_VARS].dropna() dat = dat[dat["true_interval_days"] > 0].reset_index(drop=True) # ---- 5. Published values, for comparison ------------------------------------ # From Table 2, except the two categorical F values, which are taken from the # statistician's working document of 3 December 2017; see the R script's header. PUBLISHED_F = { "gender": 0.2222, "age_baseline": 0.0961, "smoking": 0.0123, "n_teeth_baseline": 0.0674, "diagnosis": 0.6330, "true_interval_days": 37.5900, "percent_from4mm_baseline": 0.1878, "percent_from4mm_follow_up": 3.7042, "bop_percent": 283.1533, "percent_from4mm_spt": 1460.3165, } PUBLISHED_ESTIMATE = { "gender[female]": 152.3959, "age_baseline": -0.5756, "smoking[former smoker]": -122.6808, "smoking[smoker]": -512.0860, "n_teeth_baseline": 24.9403, "diagnosis[2]": -403.3057, "diagnosis[3]": -888.5687, "true_interval_days": 0.4302, "percent_from4mm_baseline": -5.5292, "percent_from4mm_follow_up": -46.6178, "bop_percent": 0.2982, "percent_from4mm_spt": 145.1780, } # ---- 6. Fit ------------------------------------------------------------------ # The outcome is rank transformed; estimates are on the rank scale and their # absolute size depends on the number of observations entering the fit. # # 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. def fit_and_report(d, label, outfile): y = average_rank(d["ppd_change"].to_numpy()) terms = [("(Intercept)", np.ones((len(d), 1)), ["(Intercept)"])] def factor(name, levels): cols = np.column_stack([(d[name] == lv).astype(float) for lv in levels]) return (name, cols, [f"{name}[{lv}]" for lv in levels]) def numeric(name): return (name, d[[name]].to_numpy(float), [name]) terms += [factor("gender", ["female"]), numeric("age_baseline"), factor("smoking", ["former smoker", "smoker"]), numeric("n_teeth_baseline"), factor("diagnosis", [2, 3]), numeric("true_interval_days"), numeric("percent_from4mm_baseline"), numeric("percent_from4mm_follow_up"), numeric("bop_percent"), numeric("percent_from4mm_spt")] X = np.hstack([block for _, block, _ in terms]) names, index, k = [], {}, 0 for name, block, labels in terms: index[name] = list(range(k, k + block.shape[1])) names += labels k += block.shape[1] beta, cov, sigma2_e, ratio = fit_random_intercept(X, y, d["pat_id"].to_numpy()) se = np.sqrt(np.diag(cov)) L1 = sequential_contrasts(X) print("\n" + "=" * 78) print(label) print("=" * 78) print(f"{len(d)} SPT visits from {d['pat_id'].nunique()} patients") print(f"between-patient / residual variance ratio = {ratio:.4f}") print("\n--- Estimates and 95% confidence intervals ---") print("These correspond to the effect sizes and intervals of Table 2.\n") print(f"{'':<32}{'Estimate':>11}{'Std.Error':>11}{'2.5 %':>11}" f"{'97.5 %':>11} {'published':>11}") for i, name in enumerate(names): pub = PUBLISHED_ESTIMATE.get(name) print(f"{name:<32}{beta[i]:>11.4f}{se[i]:>11.4f}" f"{beta[i] - 1.96 * se[i]:>11.4f}{beta[i] + 1.96 * se[i]:>11.4f} " f"{('%.4f' % pub) if pub is not None else '':>11}") print("\n--- Type I, sequential: corresponds to the F column of Table 2 ---") print("Depends on the order in which terms were entered. Not a test of an") print("independent contribution.\n") print(f"{'term':<32}{'published':>11}{'reconstructed':>15}{'ratio':>9}") type3 = {} for name, _, _ in terms: if name == "(Intercept)": continue rows = index[name] f1 = wald_F(beta, cov, L1[rows, :]) marginal = np.zeros((len(rows), X.shape[1])) marginal[range(len(rows)), rows] = 1.0 type3[name] = wald_F(beta, cov, marginal) pub = PUBLISHED_F[name] print(f"{name:<32}{pub:>11.4f}{f1:>15.4f}{f1 / pub:>9.3f}") print("\n--- Type III, marginal: use these ---") print("Each term adjusted for all others. These correspond to the confidence") print("intervals above.\n") print(f"{'term':<32}{'F value':>11}") for name in type3: print(f"{name:<32}{type3[name]:>11.4f}") pd.DataFrame({"variable": names, "estimate": beta, "std_error": se, "z_value": beta / se}).to_csv(outfile, index=False) print(f"\nWritten: {outfile}") # ---- 7. The model of Table 2 ------------------------------------------------- fit_and_report(dat[dat["subsample_5years"] == 1].reset_index(drop=True), "Table 2: five-year subsample, disease severity as recorded", "mixed_model_results_subsample.csv") # ---- 8. 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. Coefficients # are on a larger rank scale and are not comparable in size with those above. fit_and_report(dat, "All patients: not the published fit, shown for comparison", "mixed_model_results_all_patients.csv") print(f"\npandas {pd.__version__}, numpy {np.__version__}, " f"Python {sys.version.split()[0]}")