Heller UNT457 Pooled analysis from Ph1-3 studies
#libraries
rm(list=ls())
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(PKNCA)## Warning: package 'PKNCA' was built under R version 4.4.2##
## Attaching package: 'PKNCA'## The following object is masked from 'package:stats':
##
## filterlibrary(xpose4)## Warning: package 'xpose4' was built under R version 4.4.2## Loading required package: latticelibrary(knitr)
library(tidyr)## Warning: package 'tidyr' was built under R version 4.4.2library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows#theme
my_theme<-function(x){theme_bw()+
theme(text = element_text(size=20))+
theme(axis.line.y = element_line(size = 2.0))+
theme(axis.line.x = element_line(size = 2.0))+
theme(axis.ticks = element_line(size = 1.5,colour="black"))+
theme(axis.ticks.length= unit(0.45, "cm"))+
theme(axis.title.y =element_text(vjust=1.2))+
theme(axis.title.x =element_text(vjust=-0.2))+
theme(axis.text=element_text(colour="black"))+
theme(panel.background = element_rect(fill ="white"))}#data import
unt457ph123<-read.csv("/Heller/PHAR7384/wk8and9/unt457ph123.csv",stringsAsFactors = F)
unt457ph123$DOSEC<-factor(unt457ph123$DOSE,labels = c("10 mg","25 mg","50 mg","100 mg"))
unique(unt457ph123$DOSE)## [1] 10 25 50 100unt457ph123sum<-unt457ph123%>%filter(EVID==0)%>%group_by(STUDY,DOSEC,DAY,NTIME,TIME)%>%summarise(cmean=mean(DV),stdev=sd(DV))## `summarise()` has grouped output by 'STUDY', 'DOSEC', 'DAY', 'NTIME'. You can
## override using the `.groups` argument.#new theme
my_theme <- function() {
theme(
plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
strip.text = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.margin = unit(c(0.5, 0.5, 0.5,0.5), "cm")
)
}#All data population plot#mean plot all data# Plot for study 1 Rising Single Dose
cmax_data_study1 <- unt457ph123sum %>%
filter(STUDY == 1) %>%
group_by(DOSEC, STUDY, DAY) %>%
summarize(Cmax = max(cmean), Cmax_time = NTIME[which.max(cmean)], .groups = 'drop')
ggplot(data = unt457ph123sum %>% filter(STUDY == 1), aes(NTIME, cmean, color = DOSEC, group = DOSEC)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
scale_x_continuous(limits = c(0, 24), breaks = c(0, 1, 2, 4, 6, 8, 12, 24)) +
geom_errorbar(aes(ymin = cmean - stdev, ymax = cmean + stdev), width = .2) +
geom_point(data = cmax_data_study1, aes(x = Cmax_time, y = Cmax), color = 'red', size = 3, shape = 8) + # Highlight Cmax
geom_text(data = cmax_data_study1, aes(x = Cmax_time, y = Cmax, label = round(Cmax, 2)), vjust = -1, color = 'black', size = 3) + # Add Cmax values
theme_bw() +
my_theme() +
labs(x = "Time after dose (hour)", y = "Plasma concentration (ng/ml)") +
facet_wrap(vars(STUDY, DOSEC, DAY))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Plot for study 2 Rising Multiple Dose
cmax_data_study2 <- unt457ph123sum %>%
filter(STUDY == 2) %>%
group_by(DOSEC, STUDY, DAY) %>%
summarize(Cmax = max(cmean), Cmax_time = NTIME[which.max(cmean)], .groups = 'drop')
ggplot(data = unt457ph123sum %>% filter(STUDY == 2), aes(NTIME, cmean, color = DOSEC, group = DOSEC)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
scale_x_continuous(limits = c(0, 24), breaks = c(0, 1, 2, 4, 6, 8, 12, 24)) +
geom_errorbar(aes(ymin = cmean - stdev, ymax = cmean + stdev), width = .2) +
geom_point(data = cmax_data_study2, aes(x = Cmax_time, y = Cmax), color = 'red', size = 3, shape = 8) + # Highlight Cmax
geom_text(data = cmax_data_study2, aes(x = Cmax_time, y = Cmax, label = round(Cmax, 2)), vjust = -1, color = 'black', size = 3) + # Add Cmax values
theme_bw() +
my_theme() +
labs(x = "Time after dose (hour)", y = "Plasma concentration (ng/ml)") +
facet_wrap(vars(STUDY, DOSEC, DAY))
# Plot for study 3 Phase 2
cmax_data_study3 <- unt457ph123sum %>%
filter(STUDY == 3) %>%
group_by(DOSEC, STUDY, DAY) %>%
summarize(Cmax = max(cmean), Cmax_time = NTIME[which.max(cmean)], .groups = 'drop')
ggplot(data = unt457ph123sum %>% filter(STUDY == 3), aes(NTIME, cmean, color = DOSEC, group = DOSEC)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
scale_x_continuous(limits = c(0, 24), breaks = c(0, 1, 2, 4, 6, 8, 12, 24)) +
geom_errorbar(aes(ymin = cmean - stdev, ymax = cmean + stdev), width = .2) +
geom_point(data = cmax_data_study3, aes(x = Cmax_time, y = Cmax), color = 'red', size = 3, shape = 8) + # Highlight Cmax
geom_text(data = cmax_data_study3, aes(x = Cmax_time, y = Cmax, label = round(Cmax, 2)), vjust = -1, color = 'black', size = 3) + # Add Cmax values
theme_bw() +
my_theme() +
labs(x = "Time after dose (hour)", y = "Plasma concentration (ng/ml)") +
facet_wrap(vars(STUDY, DOSEC, DAY))
# Plot for study 4
cmax_data_study4 <- unt457ph123sum %>%
filter(STUDY == 4) %>%
group_by(DOSEC, STUDY, DAY) %>%
summarize(Cmax = max(cmean), Cmax_time = NTIME[which.max(cmean)], .groups = 'drop')
ggplot(data = unt457ph123sum %>% filter(STUDY == 4), aes(NTIME, cmean, color = DOSEC, group = DOSEC)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
scale_x_continuous(limits = c(0, 24), breaks = c(0, 1, 2, 4, 6, 8, 12, 24)) +
geom_errorbar(aes(ymin = cmean - stdev, ymax = cmean + stdev), width = .2) +
geom_point(data = cmax_data_study4, aes(x = Cmax_time, y = Cmax), color = 'red', size = 3, shape = 8) + # Highlight Cmax
geom_text(data = cmax_data_study4, aes(x = Cmax_time, y = Cmax, label = round(Cmax, 2)), vjust = -1, color = 'black', size = 3) + # Add Cmax values
theme_bw() +
my_theme() +
labs(x = "Time after dose (hour)", y = "Plasma concentration (ng/ml)") +
facet_wrap(vars(STUDY, DOSEC, DAY))
#Exploratory plots
#All data population plot#Covariate Summary
unt457phase123COV<-unt457ph123%>%filter(!duplicated(ID))# this selects the first row of each individual
unt457phase123COV%>%group_by(STUDY)%>%summarise(WTmean=mean(WT),WTSD=sd(WT),n=n_distinct(ID),CRCLmean=mean(CRCL),CRCLsd=sd(CRCL),Male=sum(SEX),FEMALE=n_distinct(ID)-Male)## # A tibble: 4 x 8
## STUDY WTmean WTSD n CRCLmean CRCLsd Male FEMALE
## <int> <dbl> <dbl> <int> <dbl> <dbl> <int> <int>
## 1 1 71.0 7.66 48 128. 10.3 22 26
## 2 2 72.2 7.32 48 125. 9.50 26 22
## 3 3 72.7 7.98 100 125. 9.73 52 48
## 4 4 72.5 7.68 300 125. 9.95 161 139#Diagnostics for Run31
run31<-xpose.data(31,dir="/Heller/PHAR7384/wk8and9")##
## Looking for NONMEM table files.
## Reading /Heller/PHAR7384/wk8and9/sdtab31
## Reading /Heller/PHAR7384/wk8and9/patab31
## Reading /Heller/PHAR7384/wk8and9/catab31
## Reading /Heller/PHAR7384/wk8and9/cotab31
## Table files read.
## Reading /Heller/PHAR7384/wk8and9/run31.phi
##
## Looking for NONMEM simulation table files.
## No simulated table files read.dv.vs.ipred(run31,type="p")
dv.vs.pred(run31,type="p")
cwres.vs.idv(run31,type="p")
cwres.vs.pred(run31,type="p")
ind.plots(run31)
ranpar.hist(run31)
cov.hist(run31)
ranpar.vs.cov(run31)
#pcvpc stratify by Day
#vpc run31.mod -lst=run31.lst -samples=1000 -dir=pcvpc1 -stratify=DAY -predcorr -seed=12345
pcvpc1.file <- "/Heller/PHAR7384/wk8and9/pcvpc1/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc1/vpctab31"
xpose.VPC(vpc.info=pcvpc1.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")
#stratify by study
#vpc run31.mod -lst=run31.lst -samples=1000 -dir=pcvpc2 -stratify=STUDY -predcorr -seed=12345 #error I left out -dir=pcvpc2 Ended up with vpc_dir1 file!!
vpc_dir1.file <- "/Heller/PHAR7384/wk8and9/vpc_dir1/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/vpc_dir1/vpctab31"
xpose.VPC(vpc.info=vpc_dir1.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")
#Revise run31 to run32 with additional code for stratified VPC output
#execute run32
#vpc run32.mod -lst=run32.lst -samples=1000 -dir=pcvpc2 -stratify=STRT -predcorr -seed=12345
pcvpc2.file <- "/Heller/PHAR7384/wk8and9/pcvpc2/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc2/vpctab32"
xpose.VPC(vpc.info=pcvpc2.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")
#same as above, only logy=“false”
pcvpc2.file <- "/Heller/PHAR7384/wk8and9/pcvpc2/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc2/vpctab32"
xpose.VPC(vpc.info=pcvpc2.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="FALSE")
#bootstrap # bootstrap run31.mod -samples=200 -seed=12345 -dir=bootstrap1
#library(tidyverse)
library(knitr)
# 1. Read the bootstrap simulation results
bsres <- read.csv("C:\\Heller\\PHAR7384\\wk8and9\\bootstrap1/raw_results_run31.csv", stringsAsFactors = FALSE)
# 2. Set global print options for 3 significant digits
options(digits = 3)
# 3. Extract target columns (21 to 31)
target_cols <- bsres[, 20:34]
# 4. Calculate summary metrics for each column
summary_mat <- t(sapply(target_cols, function(x) {
c(
median(x, na.rm = TRUE),
quantile(x, prob = 0.025, na.rm = TRUE),
quantile(x, prob = 0.975, na.rm = TRUE)
)
}))
# 5. Convert to data frame safely using numeric positions
summary_df <- data.frame(
parameter = rownames(summary_mat),
MEDIAN = summary_mat[, 1],
LL = summary_mat[, 2],
UL = summary_mat[, 3],
stringsAsFactors = FALSE
)
# 6. Apply your exact row ordering safely
#custom_row_order <- c(7, 6, 10, 11, 8, 2, 1, 4, 5, 3, 9,12,13,14)
custom_row_order <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
final_table <- summary_df[custom_row_order, ]
# 7. Print the final result
print(final_table, row.names = FALSE)## parameter MEDIAN LL UL
## ofv 1.70e+04 1.60e+04 1.82e+04
## KA 4.48e-01 3.80e-01 5.03e-01
## CL 4.06e+01 3.99e+01 4.14e+01
## V2 3.77e+02 3.06e+02 4.37e+02
## V3 7.72e+02 7.27e+02 8.26e+02
## Q 1.87e+02 1.78e+02 1.97e+02
## CRCLonCL 5.39e-01 3.21e-01 7.67e-01
## WTonV3 1.32e+00 1.02e+00 1.54e+00
## WTonQ 1.45e+00 1.09e+00 1.83e+00
## BSVKA 4.36e-02 2.57e-02 6.57e-02
## BSVCL 3.76e-02 3.34e-02 4.19e-02
## BSVV2 3.97e-02 3.89e-03 6.03e-02
## BSVV3 3.42e-02 2.47e-02 4.68e-02
## BSVQ 3.51e-02 1.76e-02 5.65e-02
## PROPORTIONAL 9.67e-03 9.16e-03 1.02e-02#Revise run24 to run33 with additional code for stratified VPC output to cover (P1 and P2) #execute run33
#vpc run33.mod -lst=run33.lst -samples=1000 -dir=pcvpc3 -stratify=STRT -predcorr -seed=12345
pcvpc3.file <- "/Heller/PHAR7384/wk8and9/pcvpc3/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc3/vpctab33"
xpose.VPC(vpc.info=pcvpc3.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")
#same as above run33 equiv to run 24 - only change is logy=“FALSE”
pcvpc3.file <- "/Heller/PHAR7384/wk8and9/pcvpc3/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc3/vpctab33"
xpose.VPC(vpc.info=pcvpc3.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="FALSE")
#now do a startified PCVPC for the former run 20 #add vPC code to run 20 - rename run34 #execute run 34 (modified Run 20) #vpc run34.mod -lst=run34.lst -samples=1000 -dir=pcvpc4 -stratify=STRT -predcorr -seed=12345
pcvpc4.file <- "/Heller/PHAR7384/wk8and9/pcvpc4/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc4/vpctab34"
xpose.VPC(vpc.info=pcvpc4.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")
#same as above but logy = FALSE
pcvpc4.file <- "/Heller/PHAR7384/wk8and9/pcvpc4/vpc_results.csv"
vpctab <- "/Heller/PHAR7384/wk8and9/pcvpc4/vpctab34"
xpose.VPC(vpc.info=pcvpc4.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="FALSE")