##This is synthetic patient-level clinical trial data available from: https://www.kaggle.com/datasets/dillonmyrick/bells-palsy-clinical-trial, re-created based on data from a clinical trial for corticosteroids and antiviral agents as treatment for Bell’s Palsy: https://www.nejm.org/doi/full/10.1056/nejmoa072006# ##the protocol is also available as supplementary material. ##The authors conducted a double-blind, placebo-controlled, randomized, factorial trial involving patients with Bell’s Palsy who were recruited within 72 hours after the onset of symptoms. Patients were randomly assigned to receive 10 days of treatment with prednisolone, acyclovir, both agents, or placebo. The primary outcome was recovery of facial function, as rated on the House–Brackmann scale.
##Load data: done by copying the csv file into the folder for the project.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsbells <- read.csv("BellsPalsy.csv")
glimpse(bells)## Rows: 494
## Columns: 12
## $ Patient.ID <int> 1, 2, 3, 4, 5, 6…
## $ Sex <chr> "Female", "Femal…
## $ Age <int> 77, 61, 46, 46, …
## $ Baseline.Score.on.House.Brackmann.scale <int> 6, 6, 4, 3, 3, 4…
## $ Time.between.onset.of.symptoms.and.start.of.treatment <chr> "Within 24 hr", …
## $ Treatment.Group <chr> "Prednisolone–Pl…
## $ Received.Prednisolone <chr> "Yes", "Yes", "Y…
## $ Received.Acyclovir <chr> "No", "No", "No"…
## $ X3.Month.Score.on.House.Brackmann.scale <int> 2, 1, 1, 1, 1, 1…
## $ Full.Recovery.in.3.Months <chr> "No", "Yes", "Ye…
## $ X9.Month.Score.on.House.Brackmann.scale <int> 2, 1, 1, 1, 1, 1…
## $ Full.Recovery.in.9.Months <chr> "No", "Yes", "Ye…##Cleaning the code ##Creating binary outcomes
bells <- bells |>
mutate(
prednisolone = ifelse(Received.Prednisolone == "Yes", 1,0),
acyclovir = ifelse(Received.Acyclovir == "Yes", 1,0),
threemthrecovery = ifelse(Full.Recovery.in.3.Months == "Yes", 1,0),
ninemthrecovery = ifelse(Full.Recovery.in.9.Months == "Yes", 1,0)
)##Creating baseline table
bells <- bells |>
mutate(
group = case_when(
prednisolone == 0 & acyclovir == 0 ~ "Neither",
prednisolone == 1 & acyclovir == 0 ~ "Prednisolone only",
prednisolone == 0 & acyclovir == 1 ~ "Acyclovir only",
prednisolone == 1 & acyclovir == 1 ~ "Both"
)
)##Coding time to start of treatment
bells <- bells |>
mutate(
time_group = case_when(
`Time.between.onset.of.symptoms.and.start.of.treatment` == "Within 24 hr" ~ 1,
`Time.between.onset.of.symptoms.and.start.of.treatment` == ">24 to ≤48 hr" ~ 2,
`Time.between.onset.of.symptoms.and.start.of.treatment` == ">48 to ≤72 hr"~ 3,
TRUE ~ NA_real_
)
)##Converting variables to proper types
bells <- bells |>
mutate(
sex = as.factor(Sex),
baseline_score = as.numeric(Baseline.Score.on.House.Brackmann.scale)
)##installing packages: install.packages(“gtsummary”) # run once
library(gtsummary)##Create baseline table
tbl <- bells |>
select(group, sex, Age, baseline_score, time_group) |>
tbl_summary(
by = group,
statistic = list(
all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{n} ({p}%)"
),
digits = all_continuous() ~ 1 #rounds all continuous data to one decimal point
) |>
add_overall() %>% # adds TOTAL column
modify_header(label = "**Baseline Characteristic**") |>
bold_labels()
tbl| Baseline Characteristic | Overall N = 4941 |
Acyclovir only N = 1221 |
Both N = 1271 |
Neither N = 1221 |
Prednisolone only N = 1231 |
|---|---|---|---|---|---|
| sex | |||||
| Female | 238 (48%) | 56 (46%) | 62 (49%) | 65 (53%) | 55 (45%) |
| Male | 256 (52%) | 66 (54%) | 65 (51%) | 57 (47%) | 68 (55%) |
| Age | 44.9 (14.6) | 44.6 (14.9) | 44.7 (13.4) | 45.7 (15.4) | 44.5 (14.6) |
| baseline_score | |||||
| 2 | 84 (17%) | 17 (14%) | 24 (19%) | 19 (16%) | 24 (20%) |
| 3 | 135 (27%) | 28 (23%) | 37 (29%) | 35 (29%) | 35 (28%) |
| 4 | 163 (33%) | 48 (39%) | 37 (29%) | 38 (31%) | 40 (33%) |
| 5 | 79 (16%) | 20 (16%) | 23 (18%) | 19 (16%) | 17 (14%) |
| 6 | 33 (6.7%) | 9 (7.4%) | 6 (4.7%) | 11 (9.0%) | 7 (5.7%) |
| time_group | |||||
| 1 | 249 (51%) | 68 (56%) | 62 (49%) | 66 (54%) | 53 (43%) |
| 2 | 205 (42%) | 41 (34%) | 53 (42%) | 55 (45%) | 56 (46%) |
| 3 | 39 (7.9%) | 12 (9.9%) | 12 (9.4%) | 1 (0.8%) | 14 (11%) |
| Unknown | 1 | 1 | 0 | 0 | 0 |
| 1 n (%); Mean (SD) | |||||
##Next plan to read the results section of NEJM to see which test was done. ##To test for association between P and A at 3 months. ##Recovery has been coded as HB grade 1 and is therefore binary ##Logistic regression is best. ##Testing for interaction
model <- glm(threemthrecovery ~ prednisolone * acyclovir + Age + Baseline.Score.on.House.Brackmann.scale,
family = binomial,
data = bells)
##Generating odds ratio
exp(coef(model))## (Intercept) prednisolone
## 48.8659828 3.1768886
## acyclovir Age
## 0.8541283 0.9536483
## Baseline.Score.on.House.Brackmann.scale prednisolone:acyclovir
## 0.7509903 0.8010728##Creating p values, CI and publication type tables
# Fit model (if not already saved)
model <- glm(
threemthrecovery ~ prednisolone * acyclovir + Age +
Baseline.Score.on.House.Brackmann.scale,
family = binomial,
data = bells
)
# Extract results
coef_table <- summary(model)$coefficients
OR <- exp(coef(model))
CI <- exp(confint(model))## Waiting for profiling to be done...results <- data.frame(
Variable = rownames(coef_table),
OR = OR,
CI_lower = CI[,1],
CI_upper = CI[,2],
p_value = coef_table[,4]
)
print(results)## Variable
## (Intercept) (Intercept)
## prednisolone prednisolone
## acyclovir acyclovir
## Age Age
## Baseline.Score.on.House.Brackmann.scale Baseline.Score.on.House.Brackmann.scale
## prednisolone:acyclovir prednisolone:acyclovir
## OR CI_lower CI_upper
## (Intercept) 48.8659828 17.1073121 148.9371969
## prednisolone 3.1768886 1.6781423 6.2022264
## acyclovir 0.8541283 0.4874143 1.4928699
## Age 0.9536483 0.9382095 0.9685820
## Baseline.Score.on.House.Brackmann.scale 0.7509903 0.6167942 0.9109114
## prednisolone:acyclovir 0.8010728 0.3291584 1.9333840
## p_value
## (Intercept) 1.676195e-12
## prednisolone 5.024323e-04
## acyclovir 5.800976e-01
## Age 4.858712e-09
## Baseline.Score.on.House.Brackmann.scale 3.911803e-03
## prednisolone:acyclovir 6.225793e-01##Final table will need formatting to be done.##Parking this for now. Shall revisit.