GLMMs for Epilepsy Medication: How to Get Them Right (and Wrong)

This mini-project uses simulated epilepsy data to show how generalized linear mixed models (GLMMs) can be correctly specified, and how small modeling mistakes can lead to confident but wrong conclusions.

The goal is to explain GLMMs simply and clearly, using a realistic clinical question and a sequence of Q&A-style examples.

---

Clinical Question

Does a new add-on anti-seizure medication reduce seizure frequency in drug-resistant temporal lobe epilepsy (TLE) patients over 3–6 months?

We simulate a cohort that mimics a typical clinical scenario:

• 50 patients with drug-resistant TLE.
• Each already on a background medication regimen.
• A new medication is added at time 0.
• We observe monthly seizure counts:
  • Pre-drug: 3 months.
  • Post-drug: 6 months.

Patients differ in:

• Age
• Sex
• Baseline seizure rate

We then fit a series of models and ask:

What goes wrong if we use the wrong model?
What improves when we use the right structure?

---

Data Structure

The script glmm_epilepsy_project.R simulates a data frame with one row per patient × month:

• Patient_ID: patient identifier (factor)
• Month: index of monitoring month
• RelMonth: month relative to drug start
  • Negative = baseline months
  • 0, 1, 2, ... = months after starting drug
• NewDrug: baseline vs new_drug
• Time_Since_Start: 0 before the drug, then 0, 1, 2, ... after
• Age: age in years at baseline
• Sex: F or M
• BaselineRate: patient’s underlying baseline seizure rate
• Seizure_Count: monthly seizure count
• Days_Recorded: number of days in that month (30 here)

The true data-generating model is a Negative Binomial GLMM with:

• Fixed effects for:
  • New drug
  • Time since start
  • Age
  • Sex
  • Baseline seizure rate
• Random intercept and random slope for time per patient

This lets us compare the true mechanism to what each fitted model recovers.

---

Q&A: Models from Naive to Reasonable

Q1 — Can I just compare average seizure counts before and after the drug?

Model 0: Naive Poisson GLM

model_0 <- glm(
  Seizure_Count ~ NewDrug,
  family = poisson,
  offset = log(Days_Recorded),
  data   = epilepsy_df
)

Answer

You can, but you probably should not.

• This model assumes:
  • No differences between patients.
  • No time trends.
  • No age/sex/baseline effects.
  • Poisson variance (variance approximately equal to mean).
• All repeated measures are treated as independent.

The result is usually:

• Strongly biased standard errors.
• Overconfident p-values for the drug.
• A very clean but very misleading story.

This is a good example of a model that is easy to fit and easy to misinterpret.

---

Q2 — What if I at least acknowledge that patients are different?

Model 1: Poisson GLMM with random intercept for patient

model_1 <- glmer(
  Seizure_Count ~ NewDrug + (1 | Patient_ID),
  family = poisson,
  offset = log(Days_Recorded),
  data   = epilepsy_df
)

Answer

This is already a big step forward.

• We allow each patient to have their own baseline level via (1 | Patient_ID).
• We still treat the distribution as Poisson, and still ignore time and covariates.

Moving from Model 0 to Model 1 typically shows:

• The estimated effect of the drug may change.
• The standard error for the drug effect usually increases.
• The p-value becomes less extreme (often much less dramatic).

Key idea:

Once you acknowledge between-patient variability, the apparent certainty about the treatment effect often drops.

But we still have problems:

• Overdispersion is likely.
• No explicit modeling of time.
• No adjustment for age, sex, or baseline rate.

---

Q3 — The counts look more variable than Poisson. Should I care?

Model 2: Negative Binomial GLMM with random intercept

model_2 <- glmer.nb(
  Seizure_Count ~ NewDrug + Time_Since_Start + Age + Sex +
    log(BaselineRate + 1e-6) +
    (1 | Patient_ID),
  offset = log(Days_Recorded),
  data   = epilepsy_df
)

Answer

Yes, you should care.

• Real seizure counts are usually overdispersed.
• The Negative Binomial GLMM (glmer.nb) allows variance greater than the mean.

In this model we also:

• Add time since start as a fixed effect.
• Include Age, Sex, and BaselineRate.

Now the drug effect is:

• Estimated conditional on time and baseline factors.
• More robust to extra variability in the counts.

Teaching points:

• Overdispersion inflates uncertainty. Ignoring it makes results look falsely precise.
• Adjusting for baseline seizure rate and demographics changes the interpretation:
  • Effect of the drug controlling for other factors, not just simple before/after.

---

Q4 — Do all patients improve at the same rate over time?

Model 3: Negative Binomial GLMM with random intercept and random slope

model_3 <- glmer.nb(
  Seizure_Count ~ NewDrug * Time_Since_Start + Age + Sex +
    log(BaselineRate + 1e-6) +
    (1 + Time_Since_Start | Patient_ID),
  offset = log(Days_Recorded),
  data   = epilepsy_df
)

Answer

Probably not. Patients often differ in how quickly they respond.

Here we:

• Allow each patient to have:
  • Their own baseline level (random intercept).
  • Their own time trend (random slope for Time_Since_Start).
• Include an interaction NewDrug * Time_Since_Start.

This lets us separate:

• Immediate level change at drug start.
• Change in slope over time.

With this model you can:

• Plot predicted seizure trajectories for a typical patient.
• Show patient-specific curves for a subset of patients.
• Directly visualize heterogeneity in drug response.

---

Q5 — Can I treat Age, Sex, or previous medications as random effects?

Bad idea 1: Random effect for Sex

# Not recommended
Seizure_Count ~ NewDrug + (1 | Sex) + (1 | Patient_ID)

Answer

No, not in this setting.

• Sex has only two levels (F and M).
• Random effects need many levels to estimate a variance reliably.
• Here, you actually want an explicit, interpretable fixed effect:
  • Difference between males and females.

Using (1 | Sex) with only two levels is unstable and conceptually odd.

---

Bad idea 2: Random effect for detailed medication combinations

# Not recommended if most levels are unique
Seizure_Count ~ NewDrug + (1 | MedCombo) + (1 | Patient_ID)

If each patient has a unique or nearly unique combination of previous drugs:

• (1 | MedCombo) behaves like another patient-level random intercept.
• Variance estimation becomes unstable.
• It does not give a clear, interpretable summary.

Better alternatives:

• Capture medication burden in simple fixed covariates, like:
  • Number of prior AEDs.
  • Polytherapy vs monotherapy.
• Or keep this part out of the model in this toy example, and mention it as a limitation.

---

Visualizations and GIF Ideas

The script includes example plots built with ggplot2:

1. Observed mean seizures per month (baseline vs new drug)

   • Shows mean ± standard error across RelMonth.
   • Good for a big-picture view of before vs after.

2. Predicted trajectory for a typical patient (Model 3)

   • Uses fixed effects only (re.form = NA).
   • Shows how the model expects seizures to change around drug start.

3. Subject-specific trajectories

   • Small multiples for a random subset of patients.
   • Illustrates between-patient heterogeneity.

For GIFs, you can:

• Animate the predicted curve as you move from Model 0 → Model 3.
• Animate confidence intervals shrinking or widening.
• Animate subject-specific lines overlaid with model-based averages.

---

How to Run the Script

1. Install required packages:

install.packages(c('lme4', 'MASS', 'ggplot2', 'dplyr', 'tidyr'))

2. Run the script:

source('glmm_epilepsy_project.R')

3. Inspect outputs:

• Console summaries for Models 0–3.
• Plots for:
  • Observed means over time.
  • Predicted trajectory from Model 3.
  • Example patient trajectories.

---

Key Takeaways

• GLMMs are powerful but easy to misuse.
• Small modeling choices (ignoring patients, time, or overdispersion) can change your story.
• Random effects are for grouping factors with many levels, not for small categorical variables like Sex or continuous variables like Age.
• Always ask:
  • What is my data structure?
  • What is my scientific question?
  • Does my model reflect both?