Nonlinear PK

#libraries

library(PKNCA)
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"))}

#data import and preparation for NCA

nonlinearpk<-read.csv("C:\\Heller\\PHAR7383\\week13\\nonlinearpk.csv",stringsAsFactors = F)
#for NCA and plotting
pkncadt<-nonlinearpk
pkncadt$DOSEC<-factor(pkncadt$DOSE,labels = c("10 mg","25 mg","50 mg","100 mg"))
pkncadtsum<-pkncadt%>%group_by(DOSEC,DOSE,DAY,NTIME,TIME)%>%summarise(cmean=mean(DV),stdev=sd(DV))

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

#multiple dosing, oral dosing, plots

#Population Plot
ggplot(data=pkncadt,aes(NTIME,DV,color=DOSEC,group=ID))+
geom_line(size=0.5)+
scale_x_continuous(limits = c(0,24),breaks = c(0,1,2,4,6,8,12,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(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: 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.

#average plot with standard deviation (+/-1SD)
ggplot(data=pkncadtsum,aes(NTIME,cmean,color=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)+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(DAY))

#Noncompartmental analysis- Day 1 Data (code modified from: https://cran.r-project.org/web/packages/PKNCA/vignettes/Introduction-and-Usage.html)

## Load the PK concentration data
pkday1 <-pkncadt%>%filter(DAY==1)
## Generate the dosing data
doseday1 <- pkday1%>%filter(NTIME==0)

## Create a concentration object specifying the concentration, time, and
## subject columns.  (Note that any number of grouping levels is
## supported; you are not restricted to just grouping by subject.)
conc_obj <-
  PKNCAconc(
    pkday1,
    DV~NTIME|DOSEC+DOSE+ID
  )
## Create a dosing object specifying the dose, time, and subject
## columns.  (Note that the grouping factors should be the same as or a
## subset of the grouping factors for concentration, and the grouping
## columns must have the same names between concentration and dose
## objects.)
dose_obj <-
  PKNCAdose(
    doseday1,
    DOSE~NTIME|DOSEC+ID
  )
## Combine the concentration and dosing information both to
## automatically define the intervals for NCA calculation and provide
## doses for calculations requiring dose.
data_obj <- PKNCAdata(conc_obj, dose_obj)

## Calculate the NCA parameters
results_obj <- pk.nca(data_obj)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1 points)
## Too few points for half-life calculation (min.hl.points=3 with only 1 points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1
## points)

Processing of Day1 results

#day1 results
pkday1_res<-results_obj$result
pkday1_res_sum<-pkday1_res%>%group_by(DOSEC,PPTESTCD)%>%filter(PPTESTCD=="auclast"|PPTESTCD=="cmax"|PPTESTCD=="half.life")%>%
  summarise(Average=mean(PPORRES),CV=sd(PPORRES)*100/mean(PPORRES))

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

#plotting AUClast
ggplot(pkday1_res%>%filter(PPTESTCD=="auclast"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="AUClast (ng*hr/ml)")

#plotting AUClast/DOSE
ggplot(pkday1_res%>%filter(PPTESTCD=="auclast"), aes(x=DOSEC, y=PPORRES/DOSE)) +
  geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="AUClast/DOSE")

#plotting Cmax
ggplot(pkday1_res%>%filter(PPTESTCD=="cmax"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="Cmax (ng/ml)")

#plotting cmax/DOSE
ggplot(pkday1_res%>%filter(PPTESTCD=="cmax"), aes(x=DOSEC, y=PPORRES/DOSE)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="Cmax/Dose")

#plotting half-life versus dose
ggplot(pkday1_res%>%filter(PPTESTCD=="half.life"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
scale_y_continuous(limits = c(0,100))+
labs(x="Dose (mg)",y="half life (hour)")

## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

#Noncompartmental analysis- Day 14 Data (code modified from: https://cran.r-project.org/web/packages/PKNCA/vignettes/Introduction-and-Usage.html)

## Load the PK concentration data
pkday14 <-pkncadt%>%filter(DAY==14)
## Generate the dosing data
doseday14 <- pkday14%>%filter(NTIME==0)


## Create a concentration object specifying the concentration, time, and
## subject columns.  (Note that any number of grouping levels is
## supported; you are not restricted to just grouping by subject.)
conc_obj <-
  PKNCAconc(
    pkday14,
    DV~NTIME|DOSEC+DOSE+ID
  )
## Create a dosing object specifying the dose, time, and subject
## columns.  (Note that the grouping factors should be the same as or a
## subset of the grouping factors for concentration, and the grouping
## columns must have the same names between concentration and dose
## objects.)
dose_obj <-
  PKNCAdose(
    doseday14,
    DOSE~NTIME|DOSEC+ID
  )
## Combine the concentration and dosing information both to
## automatically define the intervals for NCA calculation and provide
## doses for calculations requiring dose.
data_obj <- PKNCAdata(conc_obj, dose_obj)

## Calculate the NCA parameters
results_obj <- pk.nca(data_obj)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 1
## points)

## Warning: Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)
## Too few points for half-life calculation (min.hl.points=3 with only 2 points)

##extracting parameters from the list
pkday14_res<-results_obj$result
pkday14_res_sum<-pkday1_res%>%group_by(DOSEC,PPTESTCD)%>%filter(PPTESTCD=="auclast"|PPTESTCD=="cmax"|PPTESTCD=="tmax"|PPTESTCD=="half.life")%>%
  summarise(Average=mean(PPORRES),CV=sd(PPORRES)*100/mean(PPORRES))

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

Processing of Day14 results

#plotting AUClast
ggplot(pkday14_res%>%filter(PPTESTCD=="auclast"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="AUClast (ng*hr/ml)")

#plotting AUClast/DOSE
ggplot(pkday14_res%>%filter(PPTESTCD=="auclast"), aes(x=DOSEC, y=PPORRES/DOSE)) +
  geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="AUClast/DOSE")

#plotting Cmax
ggplot(pkday14_res%>%filter(PPTESTCD=="cmax"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="Cmax (ng/ml)")

#plotting cmax/DOSE
ggplot(pkday14_res%>%filter(PPTESTCD=="cmax"), aes(x=DOSEC, y=PPORRES/DOSE)) +
geom_boxplot()+
my_theme()+
labs(x="Dose (mg)",y="Cmax/Dose")

#plotting half-life versus dose. Be careful interpreting the half-life because the sampling schedule may not be enough to capture a long half-life
ggplot(pkday14_res%>%filter(PPTESTCD=="half.life"), aes(x=DOSEC, y=PPORRES)) +
geom_boxplot()+
my_theme()+
scale_y_continuous(limits = c(0,100))+
labs(x="Dose (mg)",y="half life (hour)")

## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

#NONMEM run processing. Linear one comp model is fit first and then nonlinear model. #linear pk

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

## 
## Looking for NONMEM table files.
##     Reading /Heller/PHAR7383/week13/sdtab2 
##     Reading /Heller/PHAR7383/week13/patab2 
##     Reading /Heller/PHAR7383/week13/catab2 
##     Reading /Heller/PHAR7383/week13/cotab2 
## Table files read.
## 
## 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,subset = "DAY==1")

ind.plots(run2,subset = "DAY==14")

ind.plots(run2)

ranpar.hist(run2)

#nonlinear pk; this may not be the best model yet and needs improvement to find optimal model by removing some parameters that are not needed. But, please interpret plots and objective function value to make model selection.

run1<-xpose.data(1,dir="/Heller/PHAR7383/week13")

## 
## Looking for NONMEM table files.
##     Reading /Heller/PHAR7383/week13/sdtab1 
##     Reading /Heller/PHAR7383/week13/patab1 
##     Reading /Heller/PHAR7383/week13/catab1 
##     Reading /Heller/PHAR7383/week13/cotab1 
## Table files read.
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.

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

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

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

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

ind.plots(run1,subset = "DAY==1")

ind.plots(run1,subset = "DAY==14")

ind.plots(run1)

ranpar.hist(run1)