PKPD Analysis
Sara Angelucci and John Shorter
2025-10-28
- Load libraries.
- Let`s read in the data, one that is all data points, another that is only first visit. Load and remove the 2 patients from the data file.
- Make a few plots,
first number of visits.
- Next plot, gender.
- Next plot, current age, full years.
- BMI histogram.
- Beneficial effect on seizures 3 months AFTER perampanel start against weight.
- Beneficial effect on seizures 3 months AFTER perampanel start against mg/kg.
- Did the Patient become seizure free on perampanel against mg/kg.
- Plotting side effects of PER against mg/kg.
- Beneficial effect on seizures 3 months AFTER perampanel start against mg/kg with all data points.
- Did the Patient become seizure free on perampanel against mg/kg with all data points.
- Plotting side effects of PER against mg/kg for all data points.
- Table 1 (in manuscript), using file with 1 visit per patient. Making a data table on patient profiles.
- Supplemental Table 1 (in manuscript), medication list with KABLE (including Perampanel).
- Supplemental Table 2 (in manuscript), 10 most frequently coadministered medications per patient.
- Supplemental Table 3 (in manuscript), adverse effects.
- Calculations.
- Basic stats using GLM and logistic regression
- Stats for analyzing with the repeated measures.
- Stats for analyzing with the repeated measures.
- Draft of sankey diagram
- Figure 1 (in manuscript), Sankey diagram.
- stats on Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED. Adding seizure type here. Remember this is only the first visit, as it labels the beneficial effect at 3 months.
- Showing that there is no major side effect of groups.
- Showing that there is no effect of groups.
- ANOVA on the comparison with blood levels and the 3 age groups.
- Plotting just the mkmg vs blood value with rec range as hlines.
- Coloring seizure type within the 3 months after start of perampanel.
- Plotting side effects of PER against mg/kg for all data points.
- Dotplot (ng/mL)/(mg/kg).
- Figure 2 (in manuscript), plasma levels based on dosage, and the effect of perampanel on seizure burden within 3 months.
- Figure 3 (in manuscript), relationship between plasma levels of perampanel and the dose of perampanel.
- Figure 4 (in manuscript), intra-patient variability in C/D-ratio for patients, sorted by increasing age using random patient ID, and adjusted perampanel levels by co-medication type.
- New suggestion from Allan to know if having a high dose resulted in a subsequent visit with a lower dose.
Load libraries.
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'tibble' was built under R version 4.4.3## Warning: package 'tidyr' was built under R version 4.4.3## Warning: package 'readr' was built under R version 4.4.3## Warning: package 'purrr' was built under R version 4.4.3## Warning: package 'dplyr' was built under R version 4.4.3## Warning: package 'forcats' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── 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 errorslibrary(ggplot2)
library(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(table1)## Warning: package 'table1' was built under R version 4.4.3##
## Attaching package: 'table1'
##
## The following objects are masked from 'package:base':
##
## units, units<-library(lme4)## Warning: package 'lme4' was built under R version 4.4.3## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(dplyr)
library(Hmisc)## Warning: package 'Hmisc' was built under R version 4.4.3##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:table1':
##
## label, label<-, units
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(table1)
library(d3Network)## Warning: package 'd3Network' was built under R version 4.4.3library(igraph)## Warning: package 'igraph' was built under R version 4.4.3##
## Attaching package: 'igraph'
##
## The following objects are masked from 'package:lubridate':
##
## %--%, union
##
## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union
##
## The following objects are masked from 'package:purrr':
##
## compose, simplify
##
## The following object is masked from 'package:tidyr':
##
## crossing
##
## The following object is masked from 'package:tibble':
##
## as_data_frame
##
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
##
## The following object is masked from 'package:base':
##
## unionlibrary(networkD3) #New one to help with Sankey diagram. ## Warning: package 'networkD3' was built under R version 4.4.3library(corrplot)## Warning: package 'corrplot' was built under R version 4.4.3## corrplot 0.95 loadedlibrary(ggcorrplot)## Warning: package 'ggcorrplot' was built under R version 4.4.3library(likert)## Warning: package 'likert' was built under R version 4.4.3## Loading required package: xtable## Warning: package 'xtable' was built under R version 4.4.3##
## Attaching package: 'xtable'
##
## The following objects are masked from 'package:Hmisc':
##
## label, label<-
##
## The following objects are masked from 'package:table1':
##
## label, label<-
##
##
## Attaching package: 'likert'
##
## The following object is masked from 'package:dplyr':
##
## recodelibrary(ggpubr)## Warning: package 'ggpubr' was built under R version 4.4.3library(ggsankey)
library(tidyr)
library(knitr)## Warning: package 'knitr' was built under R version 4.4.3library(kableExtra) ## Warning: package 'kableExtra' was built under R version 4.4.3##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(scales) ## Warning: package 'scales' was built under R version 4.4.3##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(emmeans)## Warning: package 'emmeans' was built under R version 4.4.3## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(stats)Let`s read in the data, one that is all data points, another that is only first visit. Load and remove the 2 patients from the data file.
setwd("C:/Users/saraa/Downloads/PKPD")
#setwd("~/PKPD_PER")
#"C:/Users/johnsh/OneDrive - Roskilde Universitet/Documents/PKPD_PER"
PKPD_JRS_transposed_all <- read_excel("PKPD_JRS_transposed_JRS_SA_7new.xlsx")## New names:
## • `` -> `...48`PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all[!(PKPD_JRS_transposed_all$Name %in% c('Patient 16', 'Patient 65')), ]
PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all %>%
dplyr::filter(!is.na(`Perampanel blood test value`))
PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all %>%
dplyr::filter(!is.na(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`))
PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all %>%
dplyr::filter(`Dose of PER from list of medications column` != 0)
# subset data set
PKPD_JRS_transposed <- read_excel("PKPD_JRS_transposed-Sara_JRS_1VisitxPatient.xlsx")
PKPD_JRS_transposed <- PKPD_JRS_transposed[!(PKPD_JRS_transposed$Name %in% c('Patient 16', 'Patient 65')), ]
PKPD_JRS_transposed <- PKPD_JRS_transposed %>%
dplyr::filter(!is.na(`Perampanel blood test value`))
PKPD_JRS_transposed <- PKPD_JRS_transposed %>%
dplyr::filter(!is.na(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`))
PKPD_JRS_transposed <- PKPD_JRS_transposed %>%
dplyr::filter(`Dose of PER from list of medications column` != 0)Looks like there is a now a duplicated column with the same title. R has changed this for us.
Make a few plots, first number of visits.
ggplot(data = PKPD_JRS_transposed_all, aes(x = Visit_number)) +
geom_bar() +
labs(title = "Frequency of Visit number",
x = "Visit number",
y = "Count") +
theme_minimal()
Next plot, gender.
ggplot(data = PKPD_JRS_transposed, aes(x = `Gender (0=F, 1=M)`)) +
geom_bar() +
labs(title = "Distribution of gender",
x = "Gender (0=F, 1=M)",
y = "Count") +
theme_minimal()
Next plot, current age, full years.
ggplot(data = PKPD_JRS_transposed, aes(x = `Current age, full years`)) +
geom_bar() +
labs(title = "Current age, full years",
x = "Current age, full years",
y = "Count") +
theme_minimal()
BMI histogram.
ggplot(data = PKPD_JRS_transposed, aes(x = `BMI (kg/m^2)`)) +
geom_histogram() +
labs(title = "BMI distribution",
x = "BMI",
y = "Count") +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 26 rows containing non-finite outside the scale range
## (`stat_bin()`).
Beneficial effect on seizures 3 months AFTER perampanel start against weight.
PKPD_JRS_transposed$`Weight at plasma measure` <- as.numeric(PKPD_JRS_transposed$`Weight at plasma measure`)
PKPD_JRS_transposed$`Beneficial effect (0=no, 1=yes) FIXED` <- as.factor(PKPD_JRS_transposed$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`)
ggplot(data = PKPD_JRS_transposed, aes(x = `Weight at plasma measure`,
y = `Perampanel blood test value`,
color = `Beneficial effect (0=no, 1=yes) FIXED`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Weight",
x = "Patient Weight (kg)",
y = "Plasma Level of PER",
color = "Beneficial Effect of Perampanel") +
theme_minimal(base_size = 15) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using formula = 'y ~ x'
Beneficial effect on seizures 3 months AFTER perampanel start against mg/kg.
PKPD_JRS_transposed$`PER dose` <- as.numeric(PKPD_JRS_transposed$`Dose of PER from list of medications column`)
PKPD_JRS_transposed$mgkg = PKPD_JRS_transposed$`PER dose`/PKPD_JRS_transposed$`Weight at plasma measure`
ggplot(data = PKPD_JRS_transposed, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `Beneficial effect (0=no, 1=yes) FIXED`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Beneficial effect 3 months AFTER perampanel start") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Did the Patient become seizure free on perampanel against mg/kg.
PKPD_JRS_transposed$`PER dose` <- as.numeric(PKPD_JRS_transposed$`Dose of PER from list of medications column`)
PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` <- as.factor(PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`)
PKPD_JRS_transposed$mgkg = PKPD_JRS_transposed$`PER dose`/PKPD_JRS_transposed$`Weight at plasma measure`
ggplot(data = PKPD_JRS_transposed, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `Did the Patient become seizure free on perampanel? (0=no, 1=yes)`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Did the Patient become seizure free on perampanel") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Plotting side effects of PER against mg/kg.
PKPD_JRS_transposed$`side effects FIXED` <- as.factor(PKPD_JRS_transposed$`side effects (0 = none, 1=something)`)
ggplot(data = PKPD_JRS_transposed, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `side effects FIXED`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose with side effects",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Side effects from Perampanel") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Beneficial effect on seizures 3 months AFTER perampanel start against mg/kg with all data points.
PKPD_JRS_transposed_all$`Weight at plasma measure` <- as.numeric(PKPD_JRS_transposed_all$`Weight at plasma measure`)
PKPD_JRS_transposed_all$`Beneficial effect (0=no, 1=yes) FIXED` <- as.factor(PKPD_JRS_transposed_all$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`)
PKPD_JRS_transposed_all$`PER dose` <- as.numeric(PKPD_JRS_transposed_all$`Dose of PER from list of medications column`)
PKPD_JRS_transposed_all$mgkg = PKPD_JRS_transposed_all$`PER dose`/PKPD_JRS_transposed_all$`Weight at plasma measure`
ggplot(data = PKPD_JRS_transposed_all, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `Beneficial effect (0=no, 1=yes) FIXED`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Beneficial effect 3 months AFTER perampanel start") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Did the Patient become seizure free on perampanel against mg/kg with all data points.
PKPD_JRS_transposed_all$`PER dose` <- as.numeric(PKPD_JRS_transposed_all$`Dose of PER from list of medications column`)
PKPD_JRS_transposed_all$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` <- as.factor(PKPD_JRS_transposed_all$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`)
PKPD_JRS_transposed_all$mgkg = PKPD_JRS_transposed_all$`PER dose`/PKPD_JRS_transposed_all$`Weight at plasma measure`
ggplot(data = PKPD_JRS_transposed_all, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `Did the Patient become seizure free on perampanel? (0=no, 1=yes)`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Did the Patient become seizure free on perampanel") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Plotting side effects of PER against mg/kg for all data points.
PKPD_JRS_transposed_all$`side effects FIXED` <- as.factor(PKPD_JRS_transposed_all$`side effects (0 = none, 1=something)`)
ggplot(data = PKPD_JRS_transposed_all, aes(x = `mgkg`,
y = `Perampanel blood test value`,
color = `side effects FIXED`)) +
geom_point(size = 3, alpha = 0.7) + # Adjust point size and transparency
#geom_jitter(size = 3, alpha = 0.7, width = 0.3, height = 0, seed = 666) + # Add jitter to points
geom_smooth(method = "lm", se = TRUE, linetype = "dashed", size = 1) + # Add dashed regression lines
labs(title = "Perampanel Plasma Level vs Dose with side effects",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)",
color = "Side effects from Perampanel") +
theme_minimal(base_size = 14) + # Use a larger base size for better readability
scale_color_manual(values = c("0" = "red3", "1" = "blue4"), # Custom color scheme
labels = c("No", "Yes")) + # Custom labels for the legend
theme(
plot.title = element_text(hjust = 0.5, face = "bold"), # Center and bold the title
axis.title = element_text(face = "bold"), # Bold the axis titles
legend.position = "top", # Position the legend at the top
legend.title = element_text(face = "bold"), # Bold the legend title
legend.background = element_rect(fill = "lightgray", color = NA), # Light gray background for legend
legend.key = element_rect(fill = "white", color = NA) # White background for legend keys
)## `geom_smooth()` using formula = 'y ~ x'
Table 1 (in manuscript), using file with 1 visit per patient. Making a data table on patient profiles.
# Convert to factors
PKPD_JRS_transposed$`Gender (0=F, 1=M)` <- as.factor(PKPD_JRS_transposed$`Gender (0=F, 1=M)`)
levels(PKPD_JRS_transposed$`Gender (0=F, 1=M)`) <- c("Female", "Male")
PKPD_JRS_transposed$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` <- as.factor(PKPD_JRS_transposed$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`)
levels(PKPD_JRS_transposed$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) <- c("No", "Yes")
PKPD_JRS_transposed$`side effects (0 = none, 1=something)` <- as.factor(PKPD_JRS_transposed$`side effects (0 = none, 1=something)`)
levels(PKPD_JRS_transposed$`side effects (0 = none, 1=something)`) <- c("No", "Yes")
PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` <- as.factor(PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`)
levels(PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`) <- c("No", "Yes")
PKPD_JRS_transposed$`Is the Patient currently on perampanel treatment? (Y/N). If "No" then go to row 21 (0=No, 1=Yes)` <- as.factor(PKPD_JRS_transposed$`Is the Patient currently on perampanel treatment? (Y/N). If "No" then go to row 21 (0=No, 1=Yes)`)
levels(PKPD_JRS_transposed$`Is the Patient currently on perampanel treatment? (Y/N). If "No" then go to row 21 (0=No, 1=Yes)`) <- c("No", "Yes")
# Apply descriptive labels
table1::label(PKPD_JRS_transposed$`Current age, full years`) <- "Current age (years)"
table1::label(PKPD_JRS_transposed$`Age at start of perampanel (full years)`) <- "Age at start of perampanel (years)"
table1::label(PKPD_JRS_transposed$`Gender (0=F, 1=M)`) <- "Sex"
table1::label(PKPD_JRS_transposed$`Weight at plasma measure`) <- "Weight (kg)"
table1::label(PKPD_JRS_transposed$`side effects (0 = none, 1=something)`) <- "Adverse effects"
table1::label(PKPD_JRS_transposed$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) <- "Reduced seizure burden 3 months after perampanel start"
table1::label(PKPD_JRS_transposed$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`) <- "Patient became seizure free on perampanel"
table1::label(PKPD_JRS_transposed$`Is the Patient currently on perampanel treatment? (Y/N). If "No" then go to row 21 (0=No, 1=Yes)`) <- "Patient is currently on perampanel"
table1::label(PKPD_JRS_transposed$`Total duration of perampanel treatment (in total years) (remains to be calculated)`) <- "Total duration of perampanel treatment (years)"
table1::label(PKPD_JRS_transposed$`Seizure type (focal, generalized, both focal+generalized, unknown origin)`) <- "Seizure type"
# Function to render missing values
render_missing <- function(x, ...) c("Missing" = sum(is.na(x)))
# Custom continuous rendering function:
# Show both median [Q1–Q3] and median [min–max] for age and duration,
# only median [Q1–Q3] for others.
my.render.cont <- function(x, ...) {
lbl <- attr(x, "label")
q1 <- quantile(x, 0.25, na.rm = TRUE)
q3 <- quantile(x, 0.75, na.rm = TRUE)
mn <- min(x, na.rm = TRUE)
mx <- max(x, na.rm = TRUE)
med <- median(x, na.rm = TRUE)
if (lbl %in% c("Age at start of perampanel (years)",
"Total duration of perampanel treatment (years)")) {
c("",
"Median [Q1–Q3]" = sprintf("%.1f [%.1f–%.1f]", med, q1, q3),
"Median [min–max]" = sprintf("%.1f [%.1f–%.1f]", med, mn, mx))
} else {
c("",
"Median [Q1–Q3]" = sprintf("%.1f [%.1f–%.1f]", med, q1, q3))
}
}
# Calculate total number of patients
total_patients <- nrow(PKPD_JRS_transposed[!is.na(PKPD_JRS_transposed$Name), ])
# Create age group variable
PKPD_JRS_transposed$`Age group` <- cut(
PKPD_JRS_transposed$`Age at start of perampanel (full years)`,
breaks = c(-Inf, 5, 12, Inf),
labels = c("<6", "6–12", ">12"),
right = TRUE
)
table1::label(PKPD_JRS_transposed$`Age group`) <- "Age group (age at start of perampanel)"
# Generate summary table
table_output <- table1(
~ `Age at start of perampanel (full years)` +
`Age group` +
`Gender (0=F, 1=M)` +
`Weight at plasma measure` +
`Seizure type (focal, generalized, both focal+generalized, unknown origin)` +
`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` +
`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` +
`side effects (0 = none, 1=something)` +
`Is the Patient currently on perampanel treatment? (Y/N). If "No" then go to row 21 (0=No, 1=Yes)` +
`Total duration of perampanel treatment (in total years) (remains to be calculated)`,
data = PKPD_JRS_transposed,
render.missing = render_missing,
render.continuous = my.render.cont,
overall = "Perampanel",
topclass = "Rtable1-basic",
rowlabelhead = "Patient characteristics"
)
table_output| Patient characteristics | Perampanel (N=68) |
|---|---|
| Age at start of perampanel (years) | |
| Median [Q1–Q3] | 10.5 [6.0–14.0] |
| Median [min–max] | 10.5 [1.0–18.0] |
| Age group (age at start of perampanel) | |
| <6 | 15 (22.1%) |
| 6–12 | 29 (42.6%) |
| >12 | 24 (35.3%) |
| Sex | |
| Female | 36 (52.9%) |
| Male | 32 (47.1%) |
| Weight (kg) | |
| Median [Q1–Q3] | 41.3 [25.0–54.1] |
| Seizure type | |
| Focal | 35 (51.5%) |
| Focal and generalized | 18 (26.5%) |
| Generalized | 15 (22.1%) |
| Reduced seizure burden 3 months after perampanel start | |
| No | 13 (19.1%) |
| Yes | 55 (80.9%) |
| Patient became seizure free on perampanel | |
| No | 54 (79.4%) |
| Yes | 14 (20.6%) |
| Adverse effects | |
| No | 31 (45.6%) |
| Yes | 37 (54.4%) |
| Patient is currently on perampanel | |
| No | 26 (38.2%) |
| Yes | 42 (61.8%) |
| Total duration of perampanel treatment (years) | |
| Median [Q1–Q3] | 2.8 [1.0–5.0] |
| Median [min–max] | 2.8 [0.1–12.0] |
Supplemental Table 1 (in manuscript), medication list with KABLE (including Perampanel).
medication_df <- PKPD_JRS_transposed_all %>%
select(Name, Medications_correct_to_john_alphabetic) %>%
separate_rows(Medications_correct_to_john_alphabetic, sep = ",") %>%
mutate(Medications_correct_to_john_alphabetic = trimws(Medications_correct_to_john_alphabetic))
#Total number of patients
total_patients <- n_distinct(medication_df$Name)
#Medication summary by drug (including perampanel)
medication_summary_incl <- medication_df %>%
group_by(Medications_correct_to_john_alphabetic) %>%
summarise(
`Number of Patients` = n_distinct(Name),
.groups = 'drop'
) %>%
arrange(desc(`Number of Patients`))
#percentage of patients
medication_summary_incl <- medication_summary_incl %>%
mutate(
`Percentage of Patients` = percent(`Number of Patients` / total_patients)
)
#Number of medications per patient (including perampanel)
medications_per_patient_incl <- medication_df %>%
group_by(Name) %>%
summarise(
Medications_Per_Patient = n_distinct(Medications_correct_to_john_alphabetic)
)
#statistics for medications per patient (including perampanel)
overall_median_incl <- median(medications_per_patient_incl$Medications_Per_Patient, na.rm = TRUE)
overall_iqr_incl <- quantile(medications_per_patient_incl$Medications_Per_Patient, c(0.25, 0.75), na.rm = TRUE)
# Add only the median [IQR] row (incl. perampanel)
medication_summary_incl <- medication_summary_incl %>%
add_row(
Medications_correct_to_john_alphabetic = "Median of medications administered [IQR] (incl. perampanel)",
`Number of Patients` = NA,
`Percentage of Patients` = paste0(overall_median_incl, " [", round(overall_iqr_incl[1], 1), " - ", round(overall_iqr_incl[2], 1), "]")
)
#Final medication summary
medication_summary <- medication_summary_incl %>%
distinct(Medications_correct_to_john_alphabetic, .keep_all = TRUE)
#Convert 'Number of Patients' to character to handle NA properly for summary row
medication_summary$`Number of Patients` <- as.character(medication_summary$`Number of Patients`)
medication_summary$`Number of Patients`[is.na(medication_summary$`Number of Patients`)] <- ""
# Move the statistical row to the bottom
statistical_rows <- medication_summary %>%
dplyr::filter(grepl("Median", Medications_correct_to_john_alphabetic))
medication_summary_main <- medication_summary %>%
dplyr::filter(!grepl("Median", Medications_correct_to_john_alphabetic))
#Combine main summary with the statistical row at the bottom
medication_summary_final <- bind_rows(medication_summary_main, statistical_rows)
#row to bold
rows_to_bold <- which(medication_summary_final$Medications_correct_to_john_alphabetic ==
"Median of medications administered [IQR] (incl. perampanel)")
#table
medication_summary_final %>%
kbl(
caption = "<b>Medications administered to patients</b>",
col.names = c("Medications", "Number of patients", "Percentage of patients"),
align = "lcc",
escape = FALSE
) %>%
kable_classic(full_width = FALSE, html_font = "<Calibri>") %>%
row_spec(0, bold = TRUE) %>%
row_spec(0, extra_css = "padding-top:7px;padding-bottom:7px;") %>%
row_spec(rows_to_bold, bold = TRUE)| Medications | Number of patients | Percentage of patients |
|---|---|---|
| Perampanel | 68 | 100.0% |
| Clobazam | 21 | 30.9% |
| Zonisamide | 18 | 26.5% |
| Eslicarbazepine | 17 | 25.0% |
| Lacosamide | 15 | 22.1% |
| Levetiracetam | 15 | 22.1% |
| Cannabidiol | 14 | 20.6% |
| Valproate | 12 | 17.6% |
| Brivaracetam | 9 | 13.2% |
| Topiramate | 9 | 13.2% |
| Lamotrigine | 6 | 8.8% |
| Rufinamide | 6 | 8.8% |
| Carbamazepine | 5 | 7.4% |
| Ethosuximide | 5 | 7.4% |
| Vigabatrin | 4 | 5.9% |
| Gabapentin | 3 | 4.4% |
| Oxcarbazepine | 3 | 4.4% |
| Clonazepam | 2 | 2.9% |
| Sultiame | 2 | 2.9% |
| Acetazolamide | 1 | 1.5% |
| Potassium bromide | 1 | 1.5% |
| Pregabalin | 1 | 1.5% |
| Risperidone | 1 | 1.5% |
| Median of medications administered [IQR] (incl. perampanel) | 3 [2 - 4] |
medication_summary_final## # A tibble: 24 × 3
## Medications_correct_to_john_alp…¹ `Number of Patients` Percentage of Patien…²
## <chr> <chr> <chr>
## 1 Perampanel 68 100.0%
## 2 Clobazam 21 30.9%
## 3 Zonisamide 18 26.5%
## 4 Eslicarbazepine 17 25.0%
## 5 Lacosamide 15 22.1%
## 6 Levetiracetam 15 22.1%
## 7 Cannabidiol 14 20.6%
## 8 Valproate 12 17.6%
## 9 Brivaracetam 9 13.2%
## 10 Topiramate 9 13.2%
## # ℹ 14 more rows
## # ℹ abbreviated names: ¹Medications_correct_to_john_alphabetic,
## # ²`Percentage of Patients`Supplemental Table 2 (in manuscript), 10 most frequently coadministered medications per patient.
#Combination per patient per visit
combination_df <- PKPD_JRS_transposed_all %>%
select(
Name,
Visit_number,
beneficial_raw = `Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`,
side_effect_raw = `new side effects (0 = none, 1=something)`,
Medications_correct_to_john_alphabetic
) %>%
separate_rows(Medications_correct_to_john_alphabetic, sep = ",") %>%
mutate(Medications_correct_to_john_alphabetic = trimws(Medications_correct_to_john_alphabetic)) %>%
group_by(Name, Visit_number) %>%
summarise(
Combination = paste(
c("Perampanel", sort(setdiff(unique(Medications_correct_to_john_alphabetic), "Perampanel"))),
collapse = ", "
),
reduced_seizure_burden = max(beneficial_raw, na.rm = TRUE),
adverse_effects = max(side_effect_raw, na.rm = TRUE),
.groups = 'drop'
)
#Patient-level combinations (collapse across visits)
patient_combination <- combination_df %>%
group_by(Name, Combination) %>%
summarise(
reduced_seizure_burden = ifelse(all(is.na(reduced_seizure_burden)), NA, max(reduced_seizure_burden, na.rm = TRUE)),
adverse_effects = ifelse(all(is.na(adverse_effects)), NA, max(adverse_effects, na.rm = TRUE)),
.groups = 'drop'
)
N_total_patients <- patient_combination %>%
distinct(Name) %>%
nrow() # all perampanel patients (mono + combo)
N_combo_patients <- patient_combination %>%
filter(Combination != "Perampanel") %>%
distinct(Name) %>%
nrow() # only patients with combos
# summary table (only combos)
summary_table <- patient_combination %>%
filter(Combination != "Perampanel") %>%
group_by(Combination) %>%
summarise(
`Number of patients` = n(),
`Reduced seizure burden` = sum(reduced_seizure_burden, na.rm = TRUE),
`Adverse effects` = sum(adverse_effects, na.rm = TRUE),
.groups = 'drop'
) %>%
mutate(
`% of patients with reduced seizure burden` = paste0(round((`Reduced seizure burden` / `Number of patients`) * 100, 1), "%"),
`% of patients with adverse effects` = paste0(round((`Adverse effects` / `Number of patients`) * 100, 1), "%")
) %>%
arrange(desc(`Reduced seizure burden`)) %>%
slice_head(n = 10)
#table
summary_table %>%
kbl(
caption = paste0(
"<b>First 10 most frequent co-medications administered per patient (with higher reduced seizure burden first)</b> ",
"(N = ", N_combo_patients, " patients out of ", N_total_patients, " total)"
),
align = "lccccc",
escape = FALSE
) %>%
kable_classic(full_width = FALSE, html_font = "Calibri") %>%
column_spec(1:ncol(summary_table), extra_css = "white-space: nowrap; background-color:white;") %>%
row_spec(0, bold = TRUE, extra_css = "padding-top:7px; padding-bottom:7px; background-color:white;") | Combination | Number of patients | Reduced seizure burden | Adverse effects | % of patients with reduced seizure burden | % of patients with adverse effects |
|---|---|---|---|---|---|
| Perampanel, Eslicarbazepine | 10 | 9 | 5 | 90% | 50% |
| Perampanel, Lacosamide | 5 | 5 | 2 | 100% | 40% |
| Perampanel, Zonisamide | 8 | 5 | 1 | 62.5% | 12.5% |
| Perampanel, Valproate | 3 | 3 | 1 | 100% | 33.3% |
| Perampanel, Brivaracetam, Clobazam | 2 | 2 | 1 | 100% | 50% |
| Perampanel, Clobazam, Eslicarbazepine | 2 | 2 | 2 | 100% | 100% |
| Perampanel, Lamotrigine, Zonisamide | 2 | 2 | 2 | 100% | 100% |
| Perampanel, Levetiracetam | 4 | 2 | 3 | 50% | 75% |
| Perampanel, Levetiracetam, Topiramate | 2 | 2 | 2 | 100% | 100% |
| Perampanel, Oxcarbazepine | 2 | 2 | 1 | 100% | 50% |
Supplemental Table 3 (in manuscript), adverse effects.
#Define symptom groups
side_effect_groups <- list(
"None" = c("none"),
"Central nervous system" = c("fatigue", "dizziness", "reduced cognitive performance", "increased seizures","memory difficulties"),
"Psychobehavioral" = c("aggressivity", "irritability",
"depression / sadness", "mood swings", "restlessness"),
"Gastrointestinal" = c("weight gain / increased appetite", "abdominal pain", "loose stool", "nausea"),
"Other" = c("restless sleep")
)
#Map symptoms to group
symptom_to_group <- unlist(lapply(names(side_effect_groups), function(group) {
setNames(rep(group, length(side_effect_groups[[group]])), side_effect_groups[[group]])
}))
#medication_df from patient data
medication_df <- PKPD_JRS_transposed_all %>%
select(Name, `side effects per patient`) %>%
filter(!is.na(`side effects per patient`)) %>%
separate_rows(`side effects per patient`, sep = ",") %>%
mutate(symptom = tolower(trimws(`side effects per patient`))) %>%
filter(symptom %in% names(symptom_to_group)) %>%
mutate(
Category = symptom_to_group[symptom],
Symptom = str_to_sentence(symptom)
)
# Total number of patients
total_patients <- n_distinct(PKPD_JRS_transposed_all$Name)
#Count symptoms per category
symptom_summary <- medication_df %>%
group_by(Category, Symptom) %>%
summarise(
Patients = n_distinct(Name),
.groups = "drop"
) %>%
mutate(Percentage = percent(Patients / total_patients)) %>%
arrange(match(Category, names(side_effect_groups)), Category, desc(Patients))
#Add bold category only on first row of each group
symptom_summary <- symptom_summary %>%
group_by(Category) %>%
mutate(
Category_display = ifelse(row_number() == 1, paste0("<strong>", Category, "</strong>"), "")
) %>%
ungroup()
#Calculate average side effect groups per patient
average_side_effects <- medication_df %>%
group_by(Name) %>%
summarise(UniqueGroups = n_distinct(Category)) %>%
summarise(Average = mean(UniqueGroups)) %>%
pull(Average)
#final table
final_table <- symptom_summary %>%
select(Category = Category_display, Symptom, Patients, Percentage) %>%
mutate(across(everything(), as.character)) %>%
add_row(
Category = "<strong>Average number of adverse effect groups per patient</strong>",
Symptom = "",
Patients = "",
Percentage = as.character(round(average_side_effects, 2))
)
#Render styled kable table
manual_rows_to_space <- c(1, 6, 11, 15, 16)
final_table %>%
kbl(
format = "html",
col.names = c("Category", "Adverse effects", "Number of patients", "Percentage of patients"),
caption = paste0("<b>Patients experiencing specific adverse effects</b> (N = ", total_patients, ")"),
align = "lccc",
escape = FALSE
) %>%
kable_classic(full_width = FALSE, html_font = "Calibri") %>%
row_spec(0, bold = TRUE, extra_css = "padding-top:7px; padding-bottom:7px;") %>%
row_spec(manual_rows_to_space, extra_css = "padding-bottom:12px;")| Category | Adverse effects | Number of patients | Percentage of patients |
|---|---|---|---|
| None | None | 26 | 38.2% |
| Central nervous system | Fatigue | 22 | 32.4% |
| Dizziness | 13 | 19.1% | |
| Increased seizures | 2 | 2.9% | |
| Reduced cognitive performance | 2 | 2.9% | |
| Memory difficulties | 1 | 1.5% | |
| Psychobehavioral | Aggressivity | 9 | 13.2% |
| Irritability | 4 | 5.9% | |
| Depression / sadness | 2 | 2.9% | |
| Mood swings | 1 | 1.5% | |
| Restlessness | 1 | 1.5% | |
| Gastrointestinal | Weight gain / increased appetite | 4 | 5.9% |
| Abdominal pain | 1 | 1.5% | |
| Loose stool | 1 | 1.5% | |
| Nausea | 1 | 1.5% | |
| Other | Restless sleep | 2 | 2.9% |
| Average number of adverse effect groups per patient | 1.15 |
Calculations.
# Reference range
ref_low <- 250
ref_high <- 2850
df <- PKPD_JRS_transposed_all %>%
mutate(
Perampanel_val = as.numeric(`Perampanel blood test value`),
Date_Perampanel = as.Date(`Date of Perampanel`, format = "%d/%m/%Y")
) %>%
filter(!is.na(Perampanel_val))
# Count measurements below/above reference range (all measurements)
range_counts <- df %>%
summarise(
below_n = sum(Perampanel_val < ref_low, na.rm = TRUE),
below_patients = n_distinct(Name[Perampanel_val < ref_low]),
above_n = sum(Perampanel_val > ref_high, na.rm = TRUE),
above_patients = n_distinct(Name[Perampanel_val > ref_high])
)
results <- list(
measurements_out_of_range = range_counts
)
results## $measurements_out_of_range
## # A tibble: 1 × 4
## below_n below_patients above_n above_patients
## <int> <int> <int> <int>
## 1 5 5 38 14Basic stats using GLM and logistic regression
stats on Did the Patient become seizure free on perampanel?
# Fit the logistic regression model
logistic_model <- glm(`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` ~ `Perampanel blood test value` + `Weight at plasma measure` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)` ,
data = PKPD_JRS_transposed, family = binomial)
# Display the summary of the model
summary(logistic_model)##
## Call:
## glm(formula = `Did the Patient become seizure free on perampanel? (0=no, 1=yes)` ~
## `Perampanel blood test value` + `Weight at plasma measure` +
## `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)`,
## family = binomial, data = PKPD_JRS_transposed)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -0.9337471 1.0105512 -0.924
## `Perampanel blood test value` -0.0003498 0.0003638 -0.961
## `Weight at plasma measure` -0.0154774 0.0297100 -0.521
## `Gender (0=F, 1=M)`Male 0.9261620 0.6587541 1.406
## `Age at start of perampanel (full years)` 0.0182113 0.1270660 0.143
## Pr(>|z|)
## (Intercept) 0.355
## `Perampanel blood test value` 0.336
## `Weight at plasma measure` 0.602
## `Gender (0=F, 1=M)`Male 0.160
## `Age at start of perampanel (full years)` 0.886
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 69.149 on 67 degrees of freedom
## Residual deviance: 65.212 on 63 degrees of freedom
## AIC: 75.212
##
## Number of Fisher Scoring iterations: 5# Calculate the odds ratios and 95% confidence intervals
exp_coef <- exp(coef(logistic_model)) # Odds ratios
conf_int <- exp(confint(logistic_model)) # Confidence intervals## Waiting for profiling to be done...# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 0.3930780 0.05036256 2.824748
## `Perampanel blood test value` 0.9996503 0.99878880 1.000245
## `Weight at plasma measure` 0.9846418 0.92588750 1.042691
## `Gender (0=F, 1=M)`Male 2.5248005 0.71930041 9.959287
## `Age at start of perampanel (full years)` 1.0183781 0.78801526 1.308405stats on Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED.
# Fit the logistic regression model
logistic_model <- glm(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~ `Perampanel blood test value` + `Weight at plasma measure` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)`,
data = PKPD_JRS_transposed, family = binomial)
# Display the summary of the model
summary(logistic_model)##
## Call:
## glm(formula = `Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~
## `Perampanel blood test value` + `Weight at plasma measure` +
## `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)`,
## family = binomial, data = PKPD_JRS_transposed)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 2.1331743 0.9975234 2.138
## `Perampanel blood test value` -0.0001018 0.0002682 -0.379
## `Weight at plasma measure` -0.0383054 0.0267919 -1.430
## `Gender (0=F, 1=M)`Male 0.1704953 0.6442807 0.265
## `Age at start of perampanel (full years)` 0.1060630 0.1237407 0.857
## Pr(>|z|)
## (Intercept) 0.0325 *
## `Perampanel blood test value` 0.7044
## `Weight at plasma measure` 0.1528
## `Gender (0=F, 1=M)`Male 0.7913
## `Age at start of perampanel (full years)` 0.3914
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.358 on 67 degrees of freedom
## Residual deviance: 63.555 on 63 degrees of freedom
## AIC: 73.555
##
## Number of Fisher Scoring iterations: 4# Calculate the odds ratios and 95% confidence intervals
exp_coef <- exp(coef(logistic_model)) # Odds ratios
conf_int <- exp(confint(logistic_model)) # Confidence intervals## Waiting for profiling to be done...# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 8.4416206 1.3111318 69.575238
## `Perampanel blood test value` 0.9998983 0.9993969 1.000490
## `Weight at plasma measure` 0.9624190 0.9105742 1.013315
## `Gender (0=F, 1=M)`Male 1.1858921 0.3343165 4.357660
## `Age at start of perampanel (full years)` 1.1118919 0.8791701 1.440640stats on side effects FIXED.
# Fit the logistic regression model
logistic_model <- glm(`side effects FIXED` ~ `Perampanel blood test value` + `Weight at plasma measure` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)`,
data = PKPD_JRS_transposed, family = binomial)
# Display the summary of the model
summary(logistic_model)##
## Call:
## glm(formula = `side effects FIXED` ~ `Perampanel blood test value` +
## `Weight at plasma measure` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)`,
## family = binomial, data = PKPD_JRS_transposed)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 0.5471353 0.7907390 0.692
## `Perampanel blood test value` -0.0002475 0.0002382 -1.039
## `Weight at plasma measure` -0.0305599 0.0228309 -1.339
## `Gender (0=F, 1=M)`Male -0.1823183 0.5181287 -0.352
## `Age at start of perampanel (full years)` 0.1365029 0.0995862 1.371
## Pr(>|z|)
## (Intercept) 0.489
## `Perampanel blood test value` 0.299
## `Weight at plasma measure` 0.181
## `Gender (0=F, 1=M)`Male 0.725
## `Age at start of perampanel (full years)` 0.170
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 93.738 on 67 degrees of freedom
## Residual deviance: 88.973 on 63 degrees of freedom
## AIC: 98.973
##
## Number of Fisher Scoring iterations: 4# Calculate the odds ratios and 95% confidence intervals
exp_coef <- exp(coef(logistic_model)) # Odds ratios
conf_int <- exp(confint(logistic_model)) # Confidence intervals## Waiting for profiling to be done...# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 1.7282948 0.3717614 8.555309
## `Perampanel blood test value` 0.9997525 0.9992375 1.000200
## `Weight at plasma measure` 0.9699023 0.9253943 1.013376
## `Gender (0=F, 1=M)`Male 0.8333361 0.3000421 2.318814
## `Age at start of perampanel (full years)` 1.1462582 0.9475720 1.406836Stats for analyzing with the repeated measures.
Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED.
# Ensure 'Name' and 'Visit_number' are factors
PKPD_JRS_transposed_all$Name <- as.factor(PKPD_JRS_transposed_all$Name)
PKPD_JRS_transposed_all$Visit_number <- as.factor(PKPD_JRS_transposed_all$Visit_number)
visit_counts <- PKPD_JRS_transposed_all %>%
group_by(Visit_number) %>%
tally()
# # Filter out 'Visit_number's greater than 10 and count occurrences
PKPD_JRS_filtered <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Remove visits greater than 10
group_by(Visit_number) %>%
tally()
# Further filter to remove 'Visit_number's with only one observation
PKPD_JRS_cleaned <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Keep only visits <= 10
filter(Visit_number %in% visit_counts$Visit_number[visit_counts$n > 1])
### Hmmm this is tricky to remove, because it deletes most rows.
#PKPD_JRS_cleaned_duration <- PKPD_JRS_cleaned %>%
# filter(as.numeric(`Total duration of perampanel treatment (in total years) (remains to be calculated)`) <= 5)
# Fit the mixed-effects logistic regression model
# Random effects for Name and Visit_number (nested within Name)
mixed_logistic_model <- glmer(`Beneficial effect (0=no, 1=yes) FIXED` ~ `Perampanel blood test value` + `Weight at plasma measure` + (1 | Visit_number),
data = PKPD_JRS_cleaned, family = binomial)## Warning: Some predictor variables are on very different scales: consider
## rescaling## boundary (singular) fit: see help('isSingular')# Display the summary of the model
summary(mixed_logistic_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## `Beneficial effect (0=no, 1=yes) FIXED` ~ `Perampanel blood test value` +
## `Weight at plasma measure` + (1 | Visit_number)
## Data: PKPD_JRS_cleaned
##
## AIC BIC logLik -2*log(L) df.resid
## 279.2 293.9 -135.6 271.2 284
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1766 0.3397 0.4226 0.5152 0.7383
##
## Random effects:
## Groups Name Variance Std.Dev.
## Visit_number (Intercept) 0 0
## Number of obs: 288, groups: Visit_number, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7378408 0.4298095 6.370 1.89e-10 ***
## `Perampanel blood test value` -0.0001759 0.0001183 -1.486 0.13715
## `Weight at plasma measure` -0.0219878 0.0067409 -3.262 0.00111 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `Pbtv`
## `Prmpnbtvl` -0.527
## `Wghtapmsr` -0.809 0.073
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')# Calculate odds ratios and 95% confidence intervals
exp_coef <- exp(fixef(mixed_logistic_model)) # Odds ratios for fixed effects
conf_int <- exp(confint(mixed_logistic_model, parm = "beta_", method = "Wald")) # Confidence intervals
# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 15.4535817 6.6554075 35.8825794
## `Perampanel blood test value` 0.9998241 0.9995922 1.0000560
## `Weight at plasma measure` 0.9782522 0.9654126 0.9912625Did the Patient become seizure free on perampanel? (0=no, 1=yes) including only visits 1 to 10.
# Ensure 'Name' and 'Visit_number' are factors
PKPD_JRS_transposed_all$Name <- as.factor(PKPD_JRS_transposed_all$Name)
PKPD_JRS_transposed_all$Visit_number <- as.factor(PKPD_JRS_transposed_all$Visit_number)
visit_counts <- PKPD_JRS_transposed_all %>%
group_by(Visit_number) %>%
tally()
# # Filter out 'Visit_number's greater than 10 and count occurrences
PKPD_JRS_filtered <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Remove visits greater than 10
group_by(Visit_number) %>%
tally()
# Further filter to remove 'Visit_number's with only one observation
PKPD_JRS_cleaned <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Keep only visits <= 10
filter(Visit_number %in% visit_counts$Visit_number[visit_counts$n > 1])
# Fit the mixed-effects logistic regression model
# Random effects for Name and Visit_number (nested within Name)
mixed_logistic_model <- glmer(`Did the Patient become seizure free on perampanel? (0=no, 1=yes)` ~ `Perampanel blood test value` + `Weight at plasma measure` + (1 | Visit_number),
data = PKPD_JRS_cleaned, family = binomial)## Warning: Some predictor variables are on very different scales: consider
## rescaling## boundary (singular) fit: see help('isSingular')# Display the summary of the model
summary(mixed_logistic_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: `Did the Patient become seizure free on perampanel? (0=no, 1=yes)` ~
## `Perampanel blood test value` + `Weight at plasma measure` +
## (1 | Visit_number)
## Data: PKPD_JRS_cleaned
##
## AIC BIC logLik -2*log(L) df.resid
## 272.3 286.9 -132.1 264.3 284
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7565 -0.5271 -0.4344 -0.2741 3.0296
##
## Random effects:
## Groups Name Variance Std.Dev.
## Visit_number (Intercept) 0 0
## Number of obs: 288, groups: Visit_number, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2544400 0.4211542 -0.604 0.5457
## `Perampanel blood test value` -0.0003991 0.0001782 -2.240 0.0251 *
## `Weight at plasma measure` -0.0160640 0.0078982 -2.034 0.0420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `Pbtv`
## `Prmpnbtvl` -0.563
## `Wghtapmsr` -0.751 0.021
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')# Calculate odds ratios and 95% confidence intervals
exp_coef <- exp(fixef(mixed_logistic_model)) # Odds ratios for fixed effects
conf_int <- exp(confint(mixed_logistic_model, parm = "beta_", method = "Wald")) # Confidence intervals
# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 0.7753506 0.3396339 1.7700486
## `Perampanel blood test value` 0.9996010 0.9992520 0.9999501
## `Weight at plasma measure` 0.9840644 0.9689482 0.9994163Stats for analyzing with the repeated measures.
Did the Patient become seizure free on perampanel? including only visits 1 to 10.
# Ensure 'Name' and 'Visit_number' are factors
PKPD_JRS_transposed_all$Name <- as.factor(PKPD_JRS_transposed_all$Name)
PKPD_JRS_transposed_all$Visit_number <- as.factor(PKPD_JRS_transposed_all$Visit_number)
visit_counts <- PKPD_JRS_transposed_all %>%
group_by(Visit_number) %>%
tally()
# # Filter out 'Visit_number's greater than 10 and count occurrences
PKPD_JRS_filtered <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Remove visits greater than 10
group_by(Visit_number) %>%
tally()
# Further filter to remove 'Visit_number's with only one observation
PKPD_JRS_cleaned <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 10) %>% # Keep only visits <= 10
filter(Visit_number %in% visit_counts$Visit_number[visit_counts$n > 1])
# Fit the mixed-effects logistic regression model
# Random effects for Name and Visit_number (nested within Name)
mixed_logistic_model <- glmer(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~ `Perampanel blood test value` + `Weight at plasma measure` + (1 | Visit_number),
data = PKPD_JRS_cleaned, family = binomial)## Warning: Some predictor variables are on very different scales: consider
## rescaling## boundary (singular) fit: see help('isSingular')# Display the summary of the model
summary(mixed_logistic_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## `Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~
## `Perampanel blood test value` + `Weight at plasma measure` +
## (1 | Visit_number)
## Data: PKPD_JRS_cleaned
##
## AIC BIC logLik -2*log(L) df.resid
## 279.2 293.9 -135.6 271.2 284
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1766 0.3397 0.4226 0.5152 0.7383
##
## Random effects:
## Groups Name Variance Std.Dev.
## Visit_number (Intercept) 0 0
## Number of obs: 288, groups: Visit_number, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7378408 0.4298095 6.370 1.89e-10 ***
## `Perampanel blood test value` -0.0001759 0.0001183 -1.486 0.13715
## `Weight at plasma measure` -0.0219878 0.0067409 -3.262 0.00111 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `Pbtv`
## `Prmpnbtvl` -0.527
## `Wghtapmsr` -0.809 0.073
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')# Calculate odds ratios and 95% confidence intervals
exp_coef <- exp(fixef(mixed_logistic_model)) # Odds ratios for fixed effects
conf_int <- exp(confint(mixed_logistic_model, parm = "beta_", method = "Wald")) # Confidence intervals
# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 15.4535817 6.6554075 35.8825794
## `Perampanel blood test value` 0.9998241 0.9995922 1.0000560
## `Weight at plasma measure` 0.9782522 0.9654126 0.9912625Draft of sankey diagram
Table to get the numbers for the sankey diagram, can be changed to match what to plot.
PKPD_JRS_transposed1 <- PKPD_JRS_transposed_all %>%
filter(as.numeric(Visit_number) <= 1)
table1(~ as.factor(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) + `Did the Patient become seizure free on perampanel? (0=no, 1=yes)` + `Seizure type` |
`side effects FIXED`,
data = PKPD_JRS_transposed1,
render.missing = render_missing)| 0 (N=31) |
1 (N=37) |
Overall (N=68) |
|
|---|---|---|---|
| as.factor(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) | |||
| 0 | 6 (19.4%) | 7 (18.9%) | 13 (19.1%) |
| 1 | 25 (80.6%) | 30 (81.1%) | 55 (80.9%) |
| Did the Patient become seizure free on perampanel? (0=no, 1=yes) | |||
| 0 | 27 (87.1%) | 27 (73.0%) | 54 (79.4%) |
| 1 | 4 (12.9%) | 10 (27.0%) | 14 (20.6%) |
| Seizure type | |||
| Focal | 14 (45.2%) | 21 (56.8%) | 35 (51.5%) |
| Focal and generalized | 8 (25.8%) | 10 (27.0%) | 18 (26.5%) |
| Generalized | 9 (29.0%) | 6 (16.2%) | 15 (22.1%) |
Table to get the numbers for the sankey diagram, can be changed to match what to plot.
table1(~ as.factor(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) |
`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`,
data = PKPD_JRS_transposed1,
render.missing = render_missing)| 0 (N=54) |
1 (N=14) |
Overall (N=68) |
|
|---|---|---|---|
| as.factor(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`) | |||
| 0 | 13 (24.1%) | 0 (0%) | 13 (19.1%) |
| 1 | 41 (75.9%) | 14 (100%) | 55 (80.9%) |
Figure 1 (in manuscript), Sankey diagram.
PKPD_JRS_transposed1$`Beneficial effects within 3 months` <- PKPD_JRS_transposed1$`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`
PKPD_JRS_transposed1$`Beneficial effects within 3 months` <- ifelse(PKPD_JRS_transposed1$`Beneficial effects within 3 months` == 1, "Yes", "No")
PKPD_JRS_transposed1$`Seizure free on perampanel` <- PKPD_JRS_transposed1$`Did the Patient become seizure free on perampanel? (0=no, 1=yes)`
PKPD_JRS_transposed1$`Seizure free on perampanel` <- ifelse(PKPD_JRS_transposed1$`Seizure free on perampanel` == 1, "Yes", "No")
# Renamed here
PKPD_JRS_transposed1$`Adverse effects` <- PKPD_JRS_transposed1$`side effects FIXED`
PKPD_JRS_transposed1$`Adverse effects` <- ifelse(PKPD_JRS_transposed1$`Adverse effects` == 1, "Yes", "No")
df <- PKPD_JRS_transposed1 %>%
make_long(`Seizure type`, `Beneficial effects within 3 months`, `Seizure free on perampanel`, `Adverse effects`)
ggplot(df, aes(x = x, next_x = next_x, node = node, next_node = next_node, fill = factor(node), label = node)) +
geom_sankey(flow.alpha = .6, node.color = "gray30") +
geom_sankey_label(size = 5, color = "white", fill = "gray40") +
scale_fill_viridis_d(drop = FALSE) +
theme_sankey(base_size = 18) +
labs(x = NULL) +
scale_x_discrete(position = "top") +
theme(legend.position = "none",
plot.title = element_text(hjust = .5)) +
ggtitle("")
stats on Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED. Adding seizure type here. Remember this is only the first visit, as it labels the beneficial effect at 3 months.
PKPD_JRS_transposed_duration <- PKPD_JRS_transposed1 %>%
filter(as.numeric(`Total duration of perampanel treatment (in total years) (remains to be calculated)`) <= 5)
# Set the baseline level for Seizure type
PKPD_JRS_transposed1$`Seizure type` <- as.factor(PKPD_JRS_transposed1$`Seizure type`)
PKPD_JRS_transposed1$`Seizure type` <- relevel(PKPD_JRS_transposed1$`Seizure type`, ref = "Focal and generalized")
# Fit the logistic regression model
logistic_model <- glm(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~ `Seizure type` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)` + `C:(D/kg)`,
data = PKPD_JRS_transposed1, family = binomial)
# Display the summary of the model
summary(logistic_model)##
## Call:
## glm(formula = `Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~
## `Seizure type` + `Gender (0=F, 1=M)` + `Age at start of perampanel (full years)` +
## `C:(D/kg)`, family = binomial, data = PKPD_JRS_transposed1)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 9.537e-01 1.035e+00 0.922
## `Seizure type`Focal 1.772e+00 8.289e-01 2.137
## `Seizure type`Generalized 8.980e-01 8.476e-01 1.060
## `Gender (0=F, 1=M)` 2.619e-01 7.012e-01 0.374
## `Age at start of perampanel (full years)` 1.163e-02 7.542e-02 0.154
## `C:(D/kg)` -6.648e-05 3.511e-05 -1.893
## Pr(>|z|)
## (Intercept) 0.3566
## `Seizure type`Focal 0.0326 *
## `Seizure type`Generalized 0.2894
## `Gender (0=F, 1=M)` 0.7088
## `Age at start of perampanel (full years)` 0.8775
## `C:(D/kg)` 0.0583 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.358 on 67 degrees of freedom
## Residual deviance: 57.989 on 62 degrees of freedom
## AIC: 69.989
##
## Number of Fisher Scoring iterations: 5# Calculate the odds ratios and 95% confidence intervals
exp_coef <- exp(coef(logistic_model)) # Odds ratios
conf_int <- exp(confint(logistic_model)) # Confidence intervals## Waiting for profiling to be done...# Combine odds ratios and confidence intervals into a data frame
OR_results <- data.frame(Odds_Ratio = exp_coef, Lower_CI = conf_int[,1], Upper_CI = conf_int[,2])
# Display the results
OR_results## Odds_Ratio Lower_CI Upper_CI
## (Intercept) 2.5953777 0.3460583 21.4479634
## `Seizure type`Focal 5.8814420 1.2469595 35.2229000
## `Seizure type`Generalized 2.4546489 0.4892588 14.7076989
## `Gender (0=F, 1=M)` 1.2994393 0.3323559 5.4816237
## `Age at start of perampanel (full years)` 1.0116932 0.8700178 1.1760050
## `C:(D/kg)` 0.9999335 0.9998493 0.9999947Showing that there is no major side effect of groups.
PKPD_JRS_transposed_all$`New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` <- as.factor(PKPD_JRS_transposed_all$`New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`)
PKPD_JRS_transposed_all$`New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` <- relevel(PKPD_JRS_transposed_all$`New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`, ref = "B")
# Refit the full logistic regression model if not already in environment
logistic_model <- glm(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` ,
data = PKPD_JRS_transposed_all,
family = binomial)
summary(logistic_model)##
## Call:
## glm(formula = `Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED` ~
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`,
## family = binomial, data = PKPD_JRS_transposed_all)
##
## Coefficients:
## Estimate
## (Intercept) 1.13140
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 1.68700
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.05237
## Std. Error
## (Intercept) 0.36367
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.63022
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.40160
## z value
## (Intercept) 3.111
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 2.677
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.130
## Pr(>|z|)
## (Intercept) 0.00186
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.00743
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.89625
##
## (Intercept) **
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A **
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 299.19 on 303 degrees of freedom
## Residual deviance: 285.43 on 301 degrees of freedom
## AIC: 291.43
##
## Number of Fisher Scoring iterations: 5# Get estimated marginal means (EMMs) for seizure type
emm <- emmeans(logistic_model, ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`, type = "response")
# Pairwise comparisons with odds ratios and adjusted p-values
pairs(emm, adjust = "tukey")## contrast odds.ratio SE df null z.ratio p.value
## B / A 0.185 0.117 Inf 1 -2.677 0.0203
## B / (Non-interacting) 0.949 0.381 Inf 1 -0.130 0.9907
## A / (Non-interacting) 5.128 2.780 Inf 1 3.015 0.0072
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scale# Refit the full logistic regression model if not already in environment
logistic_model <- glm(`new side effects (0 = none, 1=something)` ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` ,
data = PKPD_JRS_transposed_all,
family = binomial)
summary(logistic_model)##
## Call:
## glm(formula = `new side effects (0 = none, 1=something)` ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`,
## family = binomial, data = PKPD_JRS_transposed_all)
##
## Coefficients:
## Estimate
## (Intercept) -1.7636
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.4463
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.1542
## Std. Error
## (Intercept) 0.4419
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.5289
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.4824
## z value
## (Intercept) -3.991
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.844
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.320
## Pr(>|z|)
## (Intercept) 6.57e-05
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A 0.399
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting 0.749
##
## (Intercept) ***
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`A
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`Non-interacting
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 281.32 on 303 degrees of freedom
## Residual deviance: 280.37 on 301 degrees of freedom
## AIC: 286.37
##
## Number of Fisher Scoring iterations: 4# Get estimated marginal means (EMMs) for seizure type
emm <- emmeans(logistic_model, ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`, type = "response")
# Pairwise comparisons with odds ratios and adjusted p-values
pairs(emm, adjust = "tukey")## contrast odds.ratio SE df null z.ratio p.value
## B / A 0.640 0.339 Inf 1 -0.844 0.6758
## B / (Non-interacting) 0.857 0.414 Inf 1 -0.320 0.9453
## A / (Non-interacting) 1.339 0.468 Inf 1 0.836 0.6804
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scaleShowing that there is no effect of groups.
# Filter data to exclude values above 100000
filtered_data <- PKPD_JRS_transposed_all %>%
filter(`C:(D/kg)` <= 100000)
# One-way ANOVA: does perampanel levels differ by New grouping?
anova_model <- aov(`C:(D/kg)` ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`, data = filtered_data)
# View the summary table
summary(anova_model)## Df
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` 2
## Residuals 300
## Sum Sq
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` 2.498e+09
## Residuals 3.366e+10
## Mean Sq
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` 1.249e+09
## Residuals 1.122e+08
## F value
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` 11.13
## Residuals
## Pr(>F)
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` 2.17e-05
## Residuals
##
## `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers` ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Post-hoc pairwise comparisons with Tukey adjustment
library(emmeans)
emm <- emmeans(anova_model, ~ `New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers`)
pairs(emm, adjust = "tukey")## contrast estimate SE df t.ratio p.value
## B - A 6978 2080 300 3.358 0.0025
## B - (Non-interacting) 245 1820 300 0.135 0.9901
## A - (Non-interacting) -6733 1470 300 -4.573 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesANOVA on the comparison with blood levels and the 3 age groups.
library(dplyr)
library(emmeans)
# 1) Filter outliers and define AgeGroup (<6, 6–12, 13–17)
filtered_data <- PKPD_JRS_transposed_all %>%
filter(`C:(D/kg)` <= 100000) %>%
mutate(
AgeGroup = cut(
`Age at start of perampanel (full years)`,
breaks = c(-Inf, 6, 13, 18), # <6, 6–12, 13–17
labels = c("<6", "6–12", "13–17"),
right = FALSE
)
) %>%
filter(!is.na(AgeGroup))
# 2) One-way ANOVA
anova_mod <- aov(`C:(D/kg)` ~ AgeGroup, data = filtered_data)
summary(anova_mod)## Df Sum Sq Mean Sq F value Pr(>F)
## AgeGroup 2 1.028e+09 514140576 4.255 0.0151 *
## Residuals 275 3.323e+10 120820735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 3) Post-hoc pairwise comparisons (Tukey-adjusted)
emm <- emmeans(anova_mod, ~ AgeGroup)
pairs(emm, adjust = "tukey") # pairwise differences## contrast estimate SE df t.ratio p.value
## <6 - 6–12 -3440 1820 275 -1.888 0.1441
## <6 - 13–17 -5510 1890 275 -2.913 0.0108
## 6–12 - 13–17 -2069 1470 275 -1.408 0.3380
##
## P value adjustment: tukey method for comparing a family of 3 estimatesemm # group means (on the original scale)## AgeGroup emmean SE df lower.CL upper.CL
## <6 7013 1540 275 3983 10043
## 6–12 10453 975 275 8533 12374
## 13–17 12523 1100 275 10359 14687
##
## Confidence level used: 0.95# If you prefer base-R Tukey:
# TukeyHSD(anova_mod)Plotting just the mkmg vs blood value with rec range as hlines.
PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all %>%
mutate(value_range = ifelse(`Perampanel blood test value` >= 250 & `Perampanel blood test value` <= 2850,
"Within range", "Outside range"))
A <- ggplot(PKPD_JRS_transposed_all, aes(x = mgkg, y = `Perampanel blood test value`, color = value_range)) +
geom_point(size = 3, alpha = 0.85) +
scale_color_brewer(palette = "Set1", name = "Recommended level") +
labs(
title = " ",
x = "Perampanel Dose (mg/kg)",
y = "Plasma Perampanel (ng/mL)"
) +
geom_hline(yintercept = 250, linetype = "dashed", color = "black") +
geom_hline(yintercept = 2850, linetype = "dashed", color = "black") +
scale_y_continuous(breaks = seq(0, 6000, by = 1000), limits = c(0, 6000)) + # <-- here
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title = element_text(size = 12),
legend.position = "top",
#legend.title = element_text(face = "bold"),
legend.background = element_rect(fill = "white", color = NA),
legend.key = element_rect(fill = "white")
)Coloring seizure type within the 3 months after start of perampanel.
# Create a column with the variable names (without intercept)
OR_results$Variables <- rownames(OR_results)
#OR_results <- OR_results[-1, ] # Remove intercept
# Relabel Gender → Sex
OR_results$Variables <- gsub("Gender", "Sex", OR_results$Variables)
# Create the forest plot
B <- ggplot(OR_results, aes(x = Variables, y = Odds_Ratio, ymin = Lower_CI, ymax = Upper_CI)) +
geom_pointrange() + # Points with error bars (CIs)
geom_hline(yintercept = 1, linetype = "dashed") + # Null effect line
coord_flip() + # Flip coordinates
theme_minimal() + # Clean theme
labs(
title = " ",
x = "Variables",
y = "Odds Ratio (OR)"
) +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12)
)C <- ggplot(PKPD_JRS_transposed_all,
aes(x = as.numeric(`Dose of PER from list of medications column`),
y = `Perampanel blood test value`,
color = factor(`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`,
levels = c(0, 1),
labels = c("No", "Yes")))) +
geom_point(size = 3, alpha = 0.85) +
scale_color_brewer(palette = "Set1", name = "Reduced seizure burden") +
geom_hline(yintercept = 300, linetype = "dashed", color = "black") +
geom_hline(yintercept = 2800, linetype = "dashed", color = "black") +
labs(
title = " ",
x = "Perampanel Dose (mg)",
y = "Plasma Perampanel (ng/mL)"
) +
scale_x_continuous(
breaks = seq(0, max(as.numeric(PKPD_JRS_transposed_all$`Dose of PER from list of medications column`), na.rm = TRUE), by = 1)
) +
scale_y_continuous(
breaks = seq(0, max(PKPD_JRS_transposed_all$`Perampanel blood test value`, na.rm = TRUE), by = 1000)
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title = element_text(size = 12),
legend.position = "top",
#legend.title = element_text(face = "bold"),
legend.background = element_rect(fill = "white", color = NA),
legend.key = element_rect(fill = "white")
)Plotting side effects of PER against mg/kg for all data points.
PKPD_JRS_transposed_all$`side effects FIXED` <- factor(
PKPD_JRS_transposed_all$`side effects (0 = none, 1=something)`,
levels = c(1, 0),
labels = c("Any", "None")
)
D <- ggplot(PKPD_JRS_transposed_all,
aes(x = as.numeric(`Dose of PER from list of medications column`),
y = `Perampanel blood test value`,
color = `side effects FIXED`)) +
geom_point(size = 3, alpha = 0.85) +
scale_color_brewer(
palette = "Set1",
name = "Adverse effects"
) +
geom_hline(yintercept = 300, linetype = "dashed", color = "black") +
geom_hline(yintercept = 2800, linetype = "dashed", color = "black") +
labs(
title = " ",
x = "Perampanel Dose (mg)",
y = "Plasma Perampanel (ng/mL)"
) +
scale_x_continuous(
breaks = seq(0, max(as.numeric(PKPD_JRS_transposed_all$`Dose of PER from list of medications column`), na.rm = TRUE), by = 1)) +
scale_y_continuous(
breaks = seq(0, max(PKPD_JRS_transposed_all$`Perampanel blood test value`, na.rm = TRUE), by = 1000)
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title = element_text(size = 12),
legend.position = "top",
#legend.title = element_text(face = "bold"),
legend.background = element_rect(fill = "white", color = NA),
legend.key = element_rect(fill = "white")
)Dotplot (ng/mL)/(mg/kg).
PKPD_JRS_transposed_all <- PKPD_JRS_transposed_all %>%
mutate(
`Perampanel blood test value` = as.numeric(`Perampanel blood test value`),
Dose_mg = as.numeric(`Dose of PER from list of medications column`),
Weight_kg = as.numeric(`Weight at plasma measure`),
Age_years = as.numeric(`Age at start of perampanel (full years)`),
Benefit = factor(
`Beneficial effect on seizures 3 months AFTER perampanel start (0=no, 1=yes) FIXED`,
levels = c(0, 1), labels = c("No", "Yes")
)
) %>%
filter(
!is.na(`Perampanel blood test value`),
!is.na(Dose_mg), !is.na(Weight_kg),
Dose_mg > 0, Weight_kg > 0
) %>%
mutate(
conc_ng_ml = `Perampanel blood test value` * 349.3 / 1000,
conc_umol_L = `Perampanel blood test value` / 1000,
Dose_per_kg = Dose_mg / Weight_kg,
`C/D ratio (ng/mL*mg)` = conc_ng_ml * Dose_mg,
`C:(D/kg) ratio (umol/L/mg/kg)` = conc_umol_L / Dose_per_kg,
AgeGroup = cut(Age_years, breaks = c(-Inf, 6, 12, Inf),
labels = c("<6", "6–12", ">12"), right = FALSE),
PatientID = as.numeric(gsub("[^0-9]", "", Name))
) %>%
arrange(Age_years) %>%
mutate(PatientID = factor(PatientID, levels = unique(PatientID)))
# Plot with AgeGroup as color, Benefit as shape & alpha
E = ggplot(PKPD_JRS_transposed_all,
aes(x = PatientID, y = `C:(D/kg) ratio (umol/L/mg/kg)`)) +
geom_point(aes(color = AgeGroup, shape = Benefit, alpha = Benefit), size = 3) +
scale_color_brewer(palette = "Dark2", name = "Age group") +
scale_shape_manual(values = c("No" = 17, "Yes" = 16), name = "Reduced \nseizure burden") +
scale_alpha_manual(values = c("No" = 0.50, "Yes" = 1.00), guide = "none") +
labs(
title = " ",
x = "Patient ID (sorted by increasing age)",
y = expression("C:(D/kg) ratio ("*mu*"mol/L per mg/kg)")
) +
scale_x_discrete(labels = function(x) seq_along(x)) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
legend.position = "right",
panel.grid.major = element_line(size = 0.3),
panel.grid.minor = element_blank()
)## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.grp_col <- "New grouping A+C (inducers Carbamazepine, Phenytoin, Phenobarbital, Oxcarbazepine, Eslicarbazepine) vs B (valproate) vs non inducers"
# Data + simpler group labels (remove NA/"NA" group first)
filtered_data <- PKPD_JRS_transposed_all %>%
filter(
`C:(D/kg)` <= 100000,
!is.na(.data[[grp_col]]),
.data[[grp_col]] != "NA"
) %>%
mutate(
Group = factor(.data[[grp_col]],
levels = c("A", "B", "Non-interacting"),
labels = c("Inducers", "Valproate", "Non-inducers"))
) %>%
droplevels()
# Stats
fit <- aov(`C:(D/kg)` ~ Group, data = filtered_data)
tuk <- TukeyHSD(fit)$Group
tuk_df <- data.frame(comparison = rownames(tuk),
p = tuk[, "p adj"],
stringsAsFactors = FALSE)
p_to_stars <- function(p) ifelse(p < .001, "***",
ifelse(p < .01, "**",
ifelse(p < .05, "*", "ns")))
comparisons <- list(
c("Inducers","Valproate"),
c("Valproate","Non-inducers"),
c("Inducers","Non-inducers")
)
get_p_for <- function(a,b){
nm <- c(paste(b,a,sep="-"), paste(a,b,sep="-"))
p <- tuk_df$p[match(nm, tuk_df$comparison)]
p[!is.na(p)][1]
}
stars <- sapply(comparisons, \(x) p_to_stars(get_p_for(x[1], x[2])))
# Build the manual annotations data frame
ymax <- max(filtered_data$`C:(D/kg)`, na.rm = TRUE)
sig_df <- data.frame(
xmin = c("Inducers","Valproate","Inducers"),
xmax = c("Valproate","Non-inducers","Non-inducers"),
annotations = stars,
y_position = c(1.02, 1.12, 1.22) * ymax
)
# Plot (Y axis divided by 1000)
G = ggplot(filtered_data, aes(x = Group, y = `C:(D/kg)`, fill = Group)) +
geom_boxplot(outlier.shape = NA, width = 0.6, alpha = 0.8) +
geom_jitter(width = 0.18, size = 2, alpha = 0.6, color = "black") +
scale_fill_brewer(palette = "Set1") +
scale_y_continuous(labels = function(x) x / 1000) + # <-- scale labels
labs(
title = " ",
x = NULL, y = "C:(D/kg) ratio (μmol/L/mg/kg)"
) +
geom_signif(
data = sig_df,
aes(xmin = xmin, xmax = xmax, annotations = annotations, y_position = y_position),
manual = TRUE,
inherit.aes = FALSE,
tip_length = 0.01,
textsize = 4,
vjust = -0.1
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(size = 12)
)## Warning in geom_signif(data = sig_df, aes(xmin = xmin, xmax = xmax, annotations
## = annotations, : Ignoring unknown aesthetics: xmin, xmax, annotations, and
## y_position# see the plot
G
Figure 2 (in manuscript), plasma levels based on dosage, and the effect of perampanel on seizure burden within 3 months.
# Top: Dose vs Level, Bottom: Odds Ratios
ggarrange(A, B, labels = c("A", "B"), ncol = 1, nrow = 2, heights = c(1, 1))## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
Figure 3 (in manuscript), relationship between plasma levels of perampanel and the dose of perampanel.
# Top: Seizure-free plot, Bottom: Side effects plot
ggarrange(C, D, labels = c("A", "B"), ncol = 1, nrow = 2, heights = c(1, 1))
Figure 4 (in manuscript), intra-patient variability in C/D-ratio for patients, sorted by increasing age using random patient ID, and adjusted perampanel levels by co-medication type.
# Top: Seizure-free plot, Bottom: Side effects plot
ggarrange(E, G, labels = c("A", "B"), ncol = 1, nrow = 2, heights = c(1, 1))
New suggestion from Allan to know if having a high dose resulted in a subsequent visit with a lower dose.
suppressPackageStartupMessages({
library(dplyr)
library(readr)
})
df_raw <- PKPD_JRS_transposed_all
dose_col <- "Dose of PER from list of medications column"
level_col <- "Perampanel blood test value"
name_col <- "Name"
visit_col <- "Visit_number"
## Helper: parse numbers robustly from numeric/character/factor
safe_num <- function(x) {
if (is.numeric(x)) return(x)
readr::parse_number(as.character(x))
}
## (Optional) sanity check: see current types
str(df_raw[, c(name_col, visit_col, level_col, dose_col)])## tibble [304 × 4] (S3: tbl_df/tbl/data.frame)
## $ Name : Factor w/ 68 levels "Patient 1","Patient 10",..: 3 3 3 6 6 51 51 51 51 3 ...
## $ Visit_number : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 1 2 1 2 3 4 4 ...
## $ Perampanel blood test value : num [1:304] 930 988 1130 991 927 ...
## $ Dose of PER from list of medications column: chr [1:304] "2" "3" "4" "6" ...df <- df_raw %>%
mutate(
## Coerce types robustly
!!visit_col := safe_num(.data[[visit_col]]),
!!dose_col := safe_num(.data[[dose_col]]),
!!level_col := safe_num(.data[[level_col]])
)
df_next <- df %>%
arrange(.data[[name_col]], .data[[visit_col]]) %>%
group_by(.data[[name_col]]) %>%
mutate(
next_visit_number = dplyr::lead(.data[[visit_col]]),
next_dose = dplyr::lead(.data[[dose_col]]),
next_level = dplyr::lead(.data[[level_col]])
) %>%
ungroup()
threshold <- 2800
pairs_tbl <- df_next %>%
filter(.data[[level_col]] > threshold, !is.na(next_visit_number)) %>%
select(
{{name_col}},
current_visit = {{visit_col}},
current_level = {{level_col}},
current_dose = all_of(dose_col),
next_visit_number, next_level, next_dose
) %>%
mutate(
dose_change = next_dose - current_dose,
decreased = case_when(
is.na(current_dose) | is.na(next_dose) ~ NA,
dose_change < 0 ~ TRUE,
dose_change > 0 ~ FALSE,
TRUE ~ NA
)
) %>%
filter(!is.na(current_dose), !is.na(next_dose))
print(pairs_tbl, n = 50)## # A tibble: 32 × 9
## Name current_visit current_level current_dose next_visit_number next_level
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Patien… 6 3005 6 7 2859
## 2 Patien… 7 2859 6 8 1693
## 3 Patien… 3 4031 6 4 4455
## 4 Patien… 4 4455 6 6 2190
## 5 Patien… 3 2953 4 4 1155
## 6 Patien… 1 3509 8 2 1226
## 7 Patien… 3 2972 8 4 2074
## 8 Patien… 2 3715 12 3 2611
## 9 Patien… 4 3901 10 5 3304
## 10 Patien… 5 3304 10 6 3919
## 11 Patien… 6 3919 10 7 5426
## 12 Patien… 7 5426 10 8 3298
## 13 Patien… 8 3298 6 9 2097
## 14 Patien… 11 3020 6 12 2467
## 15 Patien… 2 5111 6 3 3237
## 16 Patien… 3 3237 6 4 5220
## 17 Patien… 1 3199 6 2 2471
## 18 Patien… 1 4109 8 2 6436
## 19 Patien… 2 6436 8 3 2134
## 20 Patien… 1 6727 10 2 5864
## 21 Patien… 2 5864 10 3 4187
## 22 Patien… 1 4955 6 2 5790
## 23 Patien… 2 5790 6 3 4386
## 24 Patien… 2 3054 8 3 1928
## 25 Patien… 2 4206 4 3 4769
## 26 Patien… 3 4769 4 4 4746
## 27 Patien… 4 4746 4 5 5474
## 28 Patien… 5 5474 4 6 4088
## 29 Patien… 6 4088 4 7 1389
## 30 Patien… 9 4674 6 10 2613
## 31 Patien… 11 3308 5 12 3477
## 32 Patien… 4 2864 10 5 1605
## # ℹ 3 more variables: next_dose <dbl>, dose_change <dbl>, decreased <lgl>desc <- pairs_tbl %>%
summarise(
n_pairs = n(),
n_decreased = sum(dose_change < 0, na.rm = TRUE),
n_increased = sum(dose_change > 0, na.rm = TRUE),
n_tied = sum(dose_change == 0, na.rm = TRUE),
median_current = median(current_dose, na.rm = TRUE),
median_next = median(next_dose, na.rm = TRUE),
median_change = median(dose_change, na.rm = TRUE),
iqr_change_low = quantile(dose_change, 0.25, na.rm = TRUE),
iqr_change_high = quantile(dose_change, 0.75, na.rm = TRUE)
)
cat("\n=== Descriptive summary ===\n"); print(desc)##
## === Descriptive summary ===## # A tibble: 1 × 9
## n_pairs n_decreased n_increased n_tied median_current median_next
## <int> <int> <int> <int> <dbl> <dbl>
## 1 32 8 1 23 6 6
## # ℹ 3 more variables: median_change <dbl>, iqr_change_low <dbl>,
## # iqr_change_high <dbl>## Tests
wilcox_res <- tryCatch(
wilcox.test(pairs_tbl$next_dose, pairs_tbl$current_dose,
paired = TRUE, alternative = "two.sided",
exact = FALSE, conf.int = TRUE),
error = function(e) e
)
cat("\n=== Wilcoxon signed-rank (next vs current) ===\n"); print(wilcox_res)##
## === Wilcoxon signed-rank (next vs current) ===##
## Wilcoxon signed rank test with continuity correction
##
## data: pairs_tbl$next_dose and pairs_tbl$current_dose
## V = 4.5, p-value = 0.0322
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## -3.000063e+00 -5.120195e-05
## sample estimates:
## (pseudo)median
## -1.999985n_dec <- sum(pairs_tbl$dose_change < 0, na.rm = TRUE)
n_inc <- sum(pairs_tbl$dose_change > 0, na.rm = TRUE)
n_eff <- n_dec + n_inc
sign_test <- if (n_eff > 0) binom.test(n_dec, n_eff, p = 0.5, alternative = "two.sided")
cat("\n=== Sign test (decrease vs increase; ties removed) ===\n"); print(sign_test)##
## === Sign test (decrease vs increase; ties removed) ===##
## Exact binomial test
##
## data: n_dec and n_eff
## number of successes = 8, number of trials = 9, p-value = 0.03906
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.5175035 0.9971909
## sample estimates:
## probability of success
## 0.8888889per_patient <- pairs_tbl %>%
mutate(direction = case_when(
dose_change < 0 ~ "decreased",
dose_change > 0 ~ "increased",
TRUE ~ "no change"
)) %>%
select({{name_col}}, current_visit, next_visit_number, current_level,
current_dose, next_dose, dose_change, direction)
cat("\n=== Per-patient pairs (triggered by level > 2800) ===\n"); print(per_patient, n = 200)##
## === Per-patient pairs (triggered by level > 2800) ===## # A tibble: 32 × 8
## Name current_visit next_visit_number current_level current_dose next_dose
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Patient… 6 7 3005 6 6
## 2 Patient… 7 8 2859 6 4
## 3 Patient… 3 4 4031 6 6
## 4 Patient… 4 6 4455 6 4
## 5 Patient… 3 4 2953 4 3.5
## 6 Patient… 1 2 3509 8 8
## 7 Patient… 3 4 2972 8 8
## 8 Patient… 2 3 3715 12 10
## 9 Patient… 4 5 3901 10 10
## 10 Patient… 5 6 3304 10 10
## 11 Patient… 6 7 3919 10 10
## 12 Patient… 7 8 5426 10 6
## 13 Patient… 8 9 3298 6 6
## 14 Patient… 11 12 3020 6 6
## 15 Patient… 2 3 5111 6 6
## 16 Patient… 3 4 3237 6 8
## 17 Patient… 1 2 3199 6 6
## 18 Patient… 1 2 4109 8 8
## 19 Patient… 2 3 6436 8 8
## 20 Patient… 1 2 6727 10 10
## 21 Patient… 2 3 5864 10 8
## 22 Patient… 1 2 4955 6 6
## 23 Patient… 2 3 5790 6 6
## 24 Patient… 2 3 3054 8 4
## 25 Patient… 2 3 4206 4 4
## 26 Patient… 3 4 4769 4 4
## 27 Patient… 4 5 4746 4 4
## 28 Patient… 5 6 5474 4 4
## 29 Patient… 6 7 4088 4 4
## 30 Patient… 9 10 4674 6 4
## 31 Patient… 11 12 3308 5 5
## 32 Patient… 4 5 2864 10 10
## # ℹ 2 more variables: dose_change <dbl>, direction <chr>