Pre-operative anti-hyperglycemic agents could mask effects

[davebridges]

Gap 4: Pre-operative medication and comorbidity status are not adjusted — likely masks glucose interaction and (more concerningly) confounds the blood pressure attenuation and liver enzyme synergy

Labels: reviewer-response analysis data-extraction priority-high

Reviewer summary

Multiple reviewers have flagged non-adjustment for pre-existing conditions and medications that could confound the main outcomes:

• Reviewer A (original Gap 4): Glucose analyses don't adjust for antihyperglycemic medications or steroidogenesis inhibitors.
• Reviewer B (formerly Gap 5): Glucose analyses don't adjust for pre-existing diabetes diagnosis.
• Reviewer C (liver enzyme confounding): ALT/AST analyses don't adjust for alcohol use, viral hepatitis, autoimmune liver disease, or preoperative hepatotoxic medications.
• Reviewer D (new — BP treatment confounding): BP analyses don't adjust for antihypertensive medication use or intensity. Reviewer proposes HDD (Hypertension Daily Dose) or Therapeutic Intensity Score quantification via ATC/RxNorm mapping, and a longitudinal trajectory analysis with IPTW for time-varying treatment.

My critique of the critiques

Legitimate concerns:

• Antihyperglycemic confounding for glucose is plausible.
• Steroidogenesis inhibitors pull cortisol down directly AND have documented hepatotoxicity (especially ketoconazole) — confound glucose, liver enzymes, and potentially BP.
• Pre-existing diabetes/hypertension/liver disease (ICD) are upstream of outcomes; ICD codes are reliably available.
• Alcohol use disorder is a genuine hepatotoxic confounder.
• Antihypertensive confounding for BP is the most plausible confounder for the most biologically striking result in the paper — this is what Reviewer D has correctly identified.

Overstated or misdirected:

• 11β-HSD1 mechanism is about biology, not medication confounding.
• Formal IPTW + marginal structural models across multiple exposure windows is overkill.
• "Case-time-control" methodology is a misapplication.
• "HCV ceiling effect" argument doesn't mechanistically explain an interaction.
• Reviewer C's ALP/GGT/bilirubin phenotyping is a different paper.
• Reviewer D's HDD/Therapeutic Intensity Score quantification is high engineering cost for marginal value over the cleaner argument. Requires RxNorm/ATC mapping, per-drug DDD lookups, dose fields (often missing or free-text in EHR), and dose-duration integration. The naive-restriction subset provides a stronger test at a fraction of the cost: untreated patients can't be confounded by treatment intensity. If obese Cushing's patients with zero antihypertensive exposure still show attenuated BP elevations vs lean untreated Cushing's patients, medication intensity is not the explanation. Skip HDD/TIS.
• Reviewer D's Option 2 (longitudinal BP trajectories with time-varying treatment, IPTW, and Cushing's×obesity×time interactions) is scope creep. This is a longitudinal paper, not a sensitivity analysis. Current paper's claim is cross-sectional; cross-sectional sensitivity tests it appropriately. Defer as follow-up.

What needs to change

A. Data extraction — three medication classes + three ICD flag bundles + antihypertensive class count

ICD flags (cheap, do first):

1. Prior diabetes — ICD-10 E10., E11., E13.; ICD-9 250.
2. Prior hypertension — ICD-10 I10–I15; ICD-9 401–405
3. Hepatotoxic/hepatic comorbidities — alcohol-related F10./Z72.1, viral hepatitis B18./B19., autoimmune K75.4/K74.3/K83.01, cirrhosis K74.6

Medication flags (180-day window), three classes — antihyperglycemics, steroidogenesis inhibitors, antihypertensives.

New: antihypertensive class count (cheap — free from existing extraction). Rather than attempt full HDD/TIS, count the number of distinct antihypertensive classes (ACE/ARB, diuretic, CCB, beta blocker, alpha blocker, MRA, central) each patient is exposed to in the 180-day window. This serves as a simple intensity proxy: 0 classes = untreated, 1 = monotherapy, 2 = dual therapy, 3+ = resistant hypertension regimen. Use this both descriptively (prevalence by stratum) and as an ordinal covariate in the BP models.

library(dplyr)
library(lubridate)
library(stringr)

# ----- ICD-based comorbidity flags -----
icd_patterns <- list(
  prior_diabetes   = c("^E10", "^E11", "^E13", "^250"),
  prior_htn        = c("^I1[0-5]", "^40[1-5]"),
  alcohol_use      = c("^F10", "^Z72\\.1", "^305\\.0", "^303"),
  viral_hepatitis  = c("^B18", "^B19", "^070"),
  autoimmune_liver = c("^K75\\.4", "^K74\\.3", "^K83\\.01", "^571\\.6", "^571\\.42"),
  cirrhosis        = c("^K74\\.6", "^571\\.5")
)

flag_prior_dx <- function(diagnoses, analysis_df, pattern, flag_name) {
  diagnoses |>
    filter(str_detect(icd_code, paste(pattern, collapse = "|"))) |>
    inner_join(analysis_df |> select(DeID_PatientID, index_date),
               by = "DeID_PatientID") |>
    filter(dx_date <= index_date) |>
    distinct(DeID_PatientID) |>
    mutate("{flag_name}" := TRUE)
}

for (nm in names(icd_patterns)) {
  analysis_df <- analysis_df |>
    left_join(flag_prior_dx(diagnoses, analysis_df, icd_patterns[[nm]], nm),
              by = "DeID_PatientID")
}

analysis_df <- analysis_df |>
  mutate(across(c(prior_diabetes, prior_htn, alcohol_use,
                  viral_hepatitis, autoimmune_liver, cirrhosis),
                \(x) replace_na(x, FALSE)),
         any_chronic_liver = viral_hepatitis | autoimmune_liver | cirrhosis)

# ----- Medication flags + antihypertensive class count -----
antihyperglycemic_rx <- c(
  "metformin", "glipizide", "glyburide", "glimepiride",
  "insulin", "semaglutide", "liraglutide", "dulaglutide", "tirzepatide",
  "empagliflozin", "dapagliflozin", "canagliflozin",
  "sitagliptin", "linagliptin", "saxagliptin",
  "pioglitazone", "rosiglitazone"
)

steroidogenesis_rx <- c(
  "ketoconazole", "metyrapone", "osilodrostat", "mitotane",
  "levoketoconazole", "mifepristone", "pasireotide"
)

antihypertensive_classes <- list(
  ace_arb    = c("lisinopril", "enalapril", "ramipril", "benazepril",
                 "losartan", "valsartan", "olmesartan", "irbesartan", "candesartan"),
  diuretic   = c("hydrochlorothiazide", "chlorthalidone", "furosemide", "torsemide",
                 "spironolactone", "eplerenone", "amiloride"),
  ccb        = c("amlodipine", "nifedipine", "felodipine", "diltiazem", "verapamil"),
  beta       = c("metoprolol", "atenolol", "carvedilol", "propranolol", "bisoprolol"),
  alpha_cent = c("doxazosin", "terazosin", "clonidine", "prazosin")
)

classify_antihypertensive_class <- function(name) {
  name_lc <- tolower(name)
  for (cls in names(antihypertensive_classes)) {
    if (str_detect(name_lc, str_c(antihypertensive_classes[[cls]], collapse = "|"))) {
      return(cls)
    }
  }
  NA_character_
}

# Count distinct antihypertensive classes per patient in the 180-day window
htn_class_counts <- medications |>
  mutate(htn_class = map_chr(med_name, classify_antihypertensive_class)) |>
  filter(!is.na(htn_class)) |>
  inner_join(analysis_df |> select(DeID_PatientID, index_date),
             by = "DeID_PatientID") |>
  filter(start_date <= index_date,
         start_date >= index_date - days(180)) |>
  distinct(DeID_PatientID, htn_class) |>
  count(DeID_PatientID, name = "n_htn_classes")

analysis_df <- analysis_df |>
  left_join(htn_class_counts, by = "DeID_PatientID") |>
  mutate(n_htn_classes = replace_na(n_htn_classes, 0L),
         antihypertensive = n_htn_classes > 0)

# Binary flags for antihyperglycemics and steroidogenesis (same as before)
# ... (existing code unchanged)

B. Descriptive table first

Supplementary table of prevalence by Cushings × obesity stratum for all flags, now including mean antihypertensive class count and distribution of 0/1/2/3+ classes:

analysis_df |>
  mutate(group = paste(if_else(Cushings == 1, "Cushing's", "Control"),
                       if_else(BMI >= 30, "obese", "lean"))) |>
  group_by(group) |>
  summarise(
    n                       = n(),
    prior_diabetes_pct      = mean(prior_diabetes) * 100,
    prior_htn_pct           = mean(prior_htn) * 100,
    alcohol_use_pct         = mean(alcohol_use) * 100,
    any_chronic_liver_pct   = mean(any_chronic_liver) * 100,
    antihyperglycemic_pct   = mean(antihyperglycemic) * 100,
    steroidogenesis_pct     = mean(steroidogenesis) * 100,
    antihypertensive_pct    = mean(antihypertensive) * 100,
    htn_classes_mean        = mean(n_htn_classes),
    htn_classes_median      = median(n_htn_classes)
  )

The htn_classes_mean row answers the intensity question directly: if obese Cushing's patients are on 2.5 classes on average while lean Cushing's patients are on 0.5, a binary flag is an underspecified adjustment and the class-count covariate matters.

C. Re-fit interaction models with the relevant covariates per outcome

Glucose and HbA1c — prior diabetes + antihyperglycemic + steroidogenesis:

lm.glucose.medadj <- lm(
  value ~ GenderCode + RaceEthnicity + Cushings * Obesity +
          prior_diabetes + antihyperglycemic + steroidogenesis,
  data = glucose.data
)

ALT and AST — alcohol use + chronic liver disease + steroidogenesis:

lm.alt.medadj <- lm(
  value ~ GenderCode + RaceEthnicity + Cushings * Obesity +
          alcohol_use + any_chronic_liver + steroidogenesis,
  data = alt.data
)

MAP, SBP, DBP — prior HTN + antihypertensive class count + steroidogenesis:

lm.map.medadj <- lm(
  MAP_imputed ~ GenderCode + RaceEthnicity + Cushings * Obesity +
                prior_htn + n_htn_classes + steroidogenesis,
  data = bp.data
)

Using class count rather than binary exposure is the intensity-adjustment response to Reviewer D. Document that this was chosen over HDD/TIS as a pragmatic alternative.

Add new rows to Table 4:

• Adjusted for prior diabetes + antihyperglycemic + steroidogenesis (glucose, HbA1c)
• Adjusted for alcohol use + chronic liver disease + steroidogenesis (ALT, AST)
• Adjusted for prior HTN + antihypertensive class count + steroidogenesis (MAP, SBP, DBP)

D. Naive-restriction sensitivity (primary answer to Reviewer D)

This is the cleaner test. For BP, restrict to patients on zero antihypertensive classes in the 180-day window AND no prior HTN diagnosis:

bp_untreated <- bp.data |>
  filter(n_htn_classes == 0, !prior_htn, !steroidogenesis)

lm.map.untreated <- lm(
  MAP_imputed ~ GenderCode + RaceEthnicity + Cushings * Obesity,
  data = bp_untreated
)

If the obesity × Cushing's interaction on BP persists in this untreated subset, treatment intensity cannot be the explanation for the attenuation — biology is. If the interaction collapses, differential treatment is the explanation. Either answer is scientifically informative and directly addresses Reviewer D's concern more cleanly than HDD quantification would.

Report the untreated-subset n per group prominently — it may be small, especially for the obese Cushing's cell.

E. Optional — Bayesian parallel with posterior probabilities

(Unchanged from previous Gap 4 version — applies here too, especially useful for the untreated BP subset where n will be small.)

F. Manuscript text

Methods (new paragraph under Statistical Analyses):

"Pre-existing diabetes, hypertension, and chronic liver conditions were ascertained from ICD-9 and ICD-10 codes with diagnosis dates preceding the relevant measurement. Pre-operative medication exposures in the 180 days before each measurement were extracted from the EHR for antihyperglycemics, steroidogenesis inhibitors, and antihypertensives. For antihypertensives, we counted distinct pharmacologic classes (ACE/ARB, diuretic, calcium channel blocker, beta blocker, alpha/central agent, mineralocorticoid receptor antagonist) as a proxy for treatment intensity. Sensitivity analyses re-estimated each interaction with these diagnoses and exposures as binary or ordinal covariates and, separately, restricted to patients with no relevant diagnosis or exposure."

Discussion (BP section — substantial revision, replaces "ceiling effect" speculation with evidence):

Replace "One possibility is that patients with obesity and Cushing's disease are already medicated to control their hypertension levels. Another possibility is there is already a ceiling effect as counter-regulatory mechanisms for blood pressure are already maximally engaged."

with: "Antihypertensive use differed substantially by obesity stratum among Cushing's patients: [X]% of obese Cushing's patients were on antihypertensives (mean [Y] classes) compared to [Z]% of lean Cushing's patients ([W] classes). After adjustment for antihypertensive class count and prior hypertension diagnosis, the obesity-attenuation of Cushing's-related blood pressure elevation [persisted / was substantially reduced / reversed] (Table 4). Restricting to patients with no antihypertensive exposure and no prior hypertension diagnosis (n = [N] per group) [confirmed / did not confirm] the attenuation pattern, [supporting biological attenuation as the dominant mechanism / suggesting that differential treatment explains much of the observed attenuation / indicating that both mechanisms contribute]."

Fill in based on results.

Abstract — hedge the blood pressure claim per Reviewer D:

Replace: "Cushing's disease also had a moderating effect on blood pressure, with participants a BMI under 30 kg/m² increasing by 12.6 mmHg and participants with obesity increasing by only 7.9 mmHg."

with: "Cushing's disease was associated with smaller clinic blood pressure increases in patients with obesity (7.9 mmHg vs 12.6 mmHg in lean patients); this pattern [persisted / attenuated] after adjustment for antihypertensive treatment."

Acceptance criteria

• ICD flag extraction code reviewed for all comorbidity categories.
• Medication extraction code reviewed for all three drug classes.
• Antihypertensive class count variable constructed and verified.
• Supplementary table: prevalence of all flags + mean/median antihypertensive class count by Cushings × obesity stratum.
• Table 4: new rows for ICD + medication-adjusted and fully-naive sensitivity analyses, for all three outcome groups.
• Methods paragraph describing ICD, medication extraction, and class-count methodology.
• Limitations paragraph updated.
• BP discussion section substantively revised — the "ceiling effect" speculation should be replaced with evidence from the actual adjusted and untreated-subset analyses. This is the single most consequential prose change across all six gaps.
• Liver discussion section updated with hepatotoxic-confounder sensitivity result.
• Abstract BP claim hedged.
• Optional but recommended: Bayesian parallel models with posterior probabilities reported for the three headline interactions, especially the untreated-BP subset where n will be small.
• Skip HDD / Therapeutic Intensity Score quantification (Reviewer D Option 1 in its formal form). Document in the response letter: "We elected to use count of antihypertensive pharmacologic classes as a pragmatic intensity proxy, and additionally to restrict to patients with no antihypertensive exposure as the primary test of treatment confounding. The untreated-subset analysis provides a direct test that cannot be confounded by treatment intensity regardless of measurement granularity."
• Skip longitudinal BP trajectory analysis with time-varying IPTW (Reviewer D Option 2). Document in the response letter: "The current paper's claim is cross-sectional; the cross-sectional untreated-subset sensitivity directly tests it. Longitudinal trajectory analysis is a separate question beyond the scope of the present manuscript."

Deferred to a follow-up paper

• Hepatocellular vs cholestatic injury phenotyping via ALP, GGT, total bilirubin (Reviewer C Option 2). Create a separate issue labeled follow-up-paper.
• Longitudinal BP trajectories with time-varying antihypertensive exposure and IPTW (Reviewer D Option 2). Same.

Notes / open questions

• Do we have access to structured medication data for controls? If not, ICD flags become the primary defense and the medication analysis is restricted or dropped.
• Untreated-BP subset n. Likely to be small, especially in the obese-Cushing's cell (obese Cushing's patients overwhelmingly will be on antihypertensives). Worth estimating before committing to this as the primary test. If n < 10 per cell, fall back to the class-count adjustment and the Bayesian posterior.
• 180-day window. 90-day secondary sensitivity worth adding if primary changes the interaction.
• Sparse confounder n. Steroidogenesis inhibitors, viral hepatitis, autoimmune liver disease will all be small. Bayesian regularization helps.

References

1. Masuzaki et al. / 11β-HSD1 adipose expression (Reviewer A — biology, don't cite for adjustment).
2. Downstream cortisol amplification (Reviewer A — same).
3. Insulin therapy effects (Reviewer A).
4. HCV ALT/AST ceiling (Reviewer C — misapplied).
5. Ketoconazole hepatotoxicity monitoring (Reviewer C — relevant).
6. HDD / Therapeutic Intensity Score methodology (Reviewer D, refs 1, 2) — cite only if defending the choice to use class count instead, or omit.

(Locate full cites before response letter.)