Heller-Phar7383_qual-Wk8-22JUN2025.knit

Heller title: “PHAR7383 Model Qualification” date: “2025-06-22” output: html_document — #libraries

library(dplyr)
library(ggplot2)
library(xpose4)
library(tidyr)
library(knitr)

#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"))}

#import data

unt234<-read.csv("C:\\Heller\\PHAR7383\\week8\\nmdtunt234.csv",stringsAsFactors = F)
unt234sum<-unt234%>%group_by(DOSE,TIME)%>%summarise(cmean=mean(DV))

## `summarise()` has grouped output by 'DOSE'. You can override using the
## `.groups` argument.

#Individual Population Plot
ggplot(data=unt234,aes(TIME,DV,group=CID))+
geom_line(size=0.5)+
geom_point(size=1)+
scale_x_continuous(limits = c(0,24),breaks = c(0,2,4,8,12,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(DOSE))

## 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: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i 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.

#Exploratory data analysis

unt234sum$DOSE<-factor(unt234sum$DOSE,labels = c("5 mg","10 mg","25 mg","50 mg","100 mg"))
##linear plot mean versus dose
ggplot(data=unt234sum,aes(TIME,cmean,color=DOSE))+
geom_line(size=0.5)+
geom_point(size=2)+
scale_x_continuous(limits = c(0,24),breaks = c(0,2,4,8,12,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")

#Twocomp Diagnostic plots (not using xpose4, just to remind you. If you want to use xpose4 for these plots, you need to include required table files in model file eg. sdtab,cotab etc.)

#unt234<-read.csv("C:\\Heller\\PHAR7383\\week8\\nmdtunt234.csv",stringsAsFactors = F)
twocomp<-read.table("C:\\Heller\\PHAR7383\\week8\\twocomp.res",skip=1,header = T)
#twocomp<-read.table("/work/ac0837/PHAR7383/wk8-qual/twocomp.res",skip=1,header = T)
#observed and individual predictions
ggplot(data=twocomp%>%filter(TIME>0),aes(DV,IPRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,500))+
  scale_y_continuous(limits = c(0,500))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Individual Predictions (ng/ml)")

#observed and population predictions
ggplot(data=twocomp%>%filter(TIME>0),aes(DV,PRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,500))+
  scale_y_continuous(limits = c(0,500))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Population Predictions (ng/ml)")

#TIME versus CRES plot
ggplot(data=twocomp%>%filter(TIME>0),aes(TIME,CWRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Conditional Weighted residuals (CWRES)")

#Predictions versus CWRES plot
ggplot(data=twocomp%>%filter(TIME>0),aes(PRED,CWRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Population Predictions (ug/L)",y="Weighted residuals (CWRES)")

#Individual plots
twocomp$DOSE<-factor(twocomp$DOSE,labels = c("5 mg","10 mg","25 mg","50 mg","100 mg"))
##individual plots 5 mg
ggplot(data=twocomp%>%filter(DOSE=="5 mg"),aes(TIME,DV))+
geom_point(size=1)+
geom_line(aes(TIME,IPRED),size=0.5,linetype = "dashed")+
scale_x_continuous(limits = c(0,24),breaks = c(0,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

##individual plots 10 mg
ggplot(data=twocomp%>%filter(DOSE=="10 mg"),aes(TIME,DV))+
geom_point(size=1)+
geom_line(aes(TIME,IPRED),size=0.5,linetype = "dashed")+
scale_x_continuous(limits = c(0,24),breaks = c(0,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

##individual plots 25 mg
ggplot(data=twocomp%>%filter(DOSE=="25 mg"),aes(TIME,DV))+
geom_point(size=1)+
geom_line(aes(TIME,IPRED),size=0.5,linetype = "dashed")+
scale_x_continuous(limits = c(0,24),breaks = c(0,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

##individual plots 50 mg
ggplot(data=twocomp%>%filter(DOSE=="50 mg"),aes(TIME,DV))+
geom_point(size=1)+
geom_line(aes(TIME,IPRED),size=0.5,linetype = "dashed")+
scale_x_continuous(limits = c(0,24),breaks = c(0,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

##individual plots 100 mg
ggplot(data=twocomp%>%filter(DOSE=="100 mg"),aes(TIME,DV))+
geom_point(size=1)+
geom_line(aes(TIME,IPRED),size=0.5,linetype = "dashed")+
scale_x_continuous(limits = c(0,24),breaks = c(0,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

#Using Xpose4

run2<-xpose.data(2,dir="/Heller/PHAR7383/week8")

## 
## Looking for NONMEM table files.
##     Reading /Heller/PHAR7383/week8/sdtab2 
##     Reading /Heller/PHAR7383/week8/patab2 
##     Reading /Heller/PHAR7383/week8/catab2 
##     Reading /Heller/PHAR7383/week8/cotab2 
## Table files read.
##     Reading /Heller/PHAR7383/week8/run2.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.

dv.vs.ipred(run2,type="p")

dv.vs.pred(run2,type="p")

cwres.vs.idv(run2,type="p")

cwres.vs.pred(run2,type="p")

ind.plots(run2)

ranpar.hist(run2)

library(xpose4)

#VPC plot

#VPC
## to be more clear about which files should be read in
vpc.file <- "C:\\Heller\\PHAR7383\\week8\\vpc/vpc_results.csv"
vpctab <- "C:\\Heller\\PHAR7383\\week8/vpc/vpctab"
xpose.VPC(vpc.info=vpc.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")

#pcVPC
vpc.file <- "C:\\Heller\\PHAR7383\\week8\\pcvpc/vpc_results.csv"
vpctab <- "C:\\Heller\\PHAR7383\\week8/pcvpc/vpctab"
xpose.VPC(vpc.info=vpc.file,vpctab=vpctab,PI.ci.area.smooth="TRUE",logy="TRUE")

#Bootstrap

bsres<-read.csv ("C:\\Heller\\PHAR7383\\week8\\bootstrap_dir1/raw_results_twocomp.csv",stringsAsFactors = F)
options(digits=2)
bsres%>%gather(parameter,value,21:31)%>%
  select(parameter, value)%>%
  group_by(parameter)%>%
summarise(MEDIAN=median(value),LL=quantile(value,prob=0.025),UL=quantile(value,prob=0.975))%>%slice(7,6,10,11,8,2,1,4,5,3,9)%>%kable()

parameter	MEDIAN	LL	UL
KA	0.33	0.16	0.62
CL	27.52	25.31	29.98
V2	98.95	48.03	180.05
V3	204.07	148.50	258.06
Q	84.84	28.64	153.30
BSVKA	0.33	0.17	0.57
BSVCL	0.07	0.04	0.11
BSVV2	0.11	0.00	0.32
BSVV3	0.33	0.09	0.83
BSVQ	0.25	0.00	0.79
SIGMA.1.1.	0.05	0.04	0.06