#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, union
library(PKNCA)
## Warning: package 'PKNCA' was built under R version 4.4.2
## 
## Attaching package: 'PKNCA'
## The following object is masked from 'package:stats':
## 
##     filter
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 analysis

wk4_tds<-read.csv("C:\\Heller\\PHAR7383\\week5\\tds\\tdsnmdt.csv",stringsAsFactors = F)
wk4_tds_sum<-wk4_tds%>%group_by(TIME,CMT)%>%summarise(cmean=mean(CONC),stdev=sd(CONC))
## `summarise()` has grouped output by 'TIME'. You can override using the
## `.groups` argument.

#Exploratory plots

#Population Plot
ggplot(data=wk4_tds,aes(TIME,CONC,group=ID))+
geom_line(size=0.5)+
geom_point(size=1)+
scale_x_continuous(limits = c(0,144),breaks = c(0,1,2,4,6,8,12,24,48,72,96,120,144))+
theme_bw()+
my_theme()+
labs(x="Time dose start (hour)",y="Plasma concentration (ng/ml)")
## 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.

#Individual Plot 
ggplot(data=wk4_tds,aes(TIME,CONC))+
geom_line(size=0.5)+
geom_point(size=1)+
scale_x_continuous(limits = c(0,144),breaks = c(0,1,2,4,6,8,12,24,48,72,96,120,144))+
theme_bw()+
my_theme()+
labs(x="Time after dose start (hour)",y="Plasma concentration (ng/ml)")+
facet_wrap(vars(ID))

#average plot with standard deviation (+/-1SD)
ggplot(data=wk4_tds_sum,aes(TIME,cmean))+
geom_line(size=0.5)+
geom_point(size=2)+
scale_x_continuous(limits = c(0,72),breaks = c(0,1,2,4,6,8,12,24,48,72))+
geom_errorbar(aes(ymin=cmean-stdev, ymax=cmean+stdev), width=.2)+
theme_bw()+
my_theme()+
labs(x="Time after dose start (hour)",y="Plasma concentration (ng/ml)")
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

#model building`

install.packages(“xpose4”)

#Run1- One comp model, zero order absorption; BSV on CL; PROP RUV

#install.packages("xpose4") 
library(xpose4)
## Warning: package 'xpose4' was built under R version 4.4.2
## Loading required package: lattice
#Diagnostic plots
run1<-xpose.data(1,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab1 
##     Reading C:\Heller\PHAR7383\week5\tds/patab1 
##     Reading C:\Heller\PHAR7383\week5\tds/catab1 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab1 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/run1.phi 
## 
## 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)

ranpar.hist(run1)

#individual parameters
run1_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab1",header=T,skip=1)%>%
  select(ID,CL,V)%>%
  filter(!duplicated(ID))
print (run1_ebe)
##    ID     CL      V
## 1   1 35.058 743.83
## 2   2 49.305 743.83
## 3   3 41.420 743.83
## 4   4 48.288 743.83
## 5   5 43.755 743.83
## 6   6 55.298 743.83
## 7   7 34.652 743.83
## 8   8 51.684 743.83
## 9   9 42.035 743.83
## 10 10 45.455 743.83
## 11 11 35.843 743.83
## 12 12 44.330 743.83

#Run2- One comp model, zero order absorption; BSV on CL,V; PROP RUV

#Diagnostic plots
run2<-xpose.data(2,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab2 
##     Reading C:\Heller\PHAR7383\week5\tds/patab2 
##     Reading C:\Heller\PHAR7383\week5\tds/catab2 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab2 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/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)

#individual parameters
run2_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab2",header=T,skip=1)%>%
  select(ID,CL,V)%>%
  filter(!duplicated(ID))
print (run2_ebe)
##    ID     CL      V
## 1   1 38.169 827.78
## 2   2 49.511 749.90
## 3   3 42.650 771.18
## 4   4 46.444 711.60
## 5   5 42.818 724.88
## 6   6 47.152 619.89
## 7   7 37.002 806.00
## 8   8 50.538 726.92
## 9   9 41.521 732.61
## 10 10 45.011 736.21
## 11 11 38.414 810.12
## 12 12 48.017 823.05

print (run2_ebe)

#Run3- Two comp model, zero order absorption; BSV on CL; PROP RUV

#Diagnostic plots
run3<-xpose.data(3,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab3 
##     Reading C:\Heller\PHAR7383\week5\tds/patab3 
##     Reading C:\Heller\PHAR7383\week5\tds/catab3 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab3 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/run3.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run3,type="p")

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

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

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

ind.plots(run3)

ranpar.hist(run3)

#individual parameters
run3_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab3",header=T,skip=1)%>%
  select(ID,CL,V1,Q,V2)%>%
  filter(!duplicated(ID))
print (run3_ebe)
##    ID      CL     V1      Q     V2
## 1   1  53.869 359.76 162.88 859.33
## 2   2  87.959 359.76 162.88 859.33
## 3   3  69.696 359.76 162.88 859.33
## 4   4  85.003 359.76 162.88 859.33
## 5   5  73.658 359.76 162.88 859.33
## 6   6 103.640 359.76 162.88 859.33
## 7   7  53.704 359.76 162.88 859.33
## 8   8  93.548 359.76 162.88 859.33
## 9   9  70.269 359.76 162.88 859.33
## 10 10  77.797 359.76 162.88 859.33
## 11 11  56.542 359.76 162.88 859.33
## 12 12  74.420 359.76 162.88 859.33

#Run4- Two comp model, zero order absorption; BSV on CL,V1; PROP RUV

#Diagnostic plots
run4<-xpose.data(4,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab4 
##     Reading C:\Heller\PHAR7383\week5\tds/patab4 
##     Reading C:\Heller\PHAR7383\week5\tds/catab4 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab4 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/run4.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run4,type="p")

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

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

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

ind.plots(run4)

ranpar.hist(run4)

#individual parameters
run4_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab4",header=T,skip=1)%>%
  select(ID,CL,V1,Q,V2)%>%
  filter(!duplicated(ID))
print (run4_ebe)
##    ID     CL     V1      Q     V2
## 1   1 55.643 463.91 147.31 829.18
## 2   2 89.483 416.42 147.31 829.18
## 3   3 71.188 427.71 147.31 829.18
## 4   4 83.066 335.78 147.31 829.18
## 5   5 71.872 336.84 147.31 829.18
## 6   6 96.313 227.00 147.31 829.18
## 7   7 53.999 403.91 147.31 829.18
## 8   8 92.728 362.03 147.31 829.18
## 9   9 69.222 357.25 147.31 829.18
## 10 10 76.522 351.25 147.31 829.18
## 11 11 58.086 446.10 147.31 829.18
## 12 12 80.429 556.03 147.31 829.18

#Run5- Two comp model, zero order absortion; BSV on CL,V1,Q; PROP RUV

#Diagnostic plots
run5<-xpose.data(5,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab5 
##     Reading C:\Heller\PHAR7383\week5\tds/patab5 
##     Reading C:\Heller\PHAR7383\week5\tds/catab5 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab5 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/run5.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run5,type="p")

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

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

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

ind.plots(run5)

ranpar.hist(run5)

#individual parameters
run5_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab5",header=T,skip=1)%>%
  select(ID,CL,V1,Q,V2)%>%
  filter(!duplicated(ID))
print (run5_ebe)
##    ID     CL     V1      Q     V2
## 1   1 55.080 450.43 167.68 836.51
## 2   2 89.599 413.05 149.96 836.51
## 3   3 73.602 436.36 124.48 836.51
## 4   4 84.300 335.94 140.61 836.51
## 5   5 71.547 330.50 155.10 836.51
## 6   6 97.993 224.82 141.60 836.51
## 7   7 56.785 434.09 103.11 836.51
## 8   8 87.629 340.46 199.58 836.51
## 9   9 68.274 347.24 165.14 836.51
## 10 10 77.306 350.88 142.85 836.51
## 11 11 59.982 462.42 118.81 836.51
## 12 12 78.289 533.42 185.58 836.51

#Run6- Two comp model, zero order absortion; BSV on CL,V1,Q,V2;PROP RUV

#Diagnostic plots
run6<-xpose.data(6,dir="C:\\Heller\\PHAR7383\\week5\\tds")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5\tds/sdtab6 
##     Reading C:\Heller\PHAR7383\week5\tds/patab6 
##     Reading C:\Heller\PHAR7383\week5\tds/catab6 
##     Reading C:\Heller\PHAR7383\week5\tds/cotab6 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5\tds/run6.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run6,type="p")

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

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

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

ind.plots(run6)

ranpar.hist(run6)

#individual parameters
run6_ebe<-read.table("C:\\Heller\\PHAR7383\\week5\\tds/patab6",header=T,skip=1)%>%
  select(ID,CL,V1,Q,V2)%>%
  filter(!duplicated(ID))
print (run6_ebe)
##    ID     CL     V1      Q      V2
## 1   1 55.440 449.12 165.22  848.07
## 2   2 88.609 421.24 143.45  810.81
## 3   3 82.225 382.41 156.16 1071.50
## 4   4 83.038 343.20 134.21  806.96
## 5   5 70.771 334.99 150.45  816.14
## 6   6 95.710 233.52 130.53  791.36
## 7   7 57.795 426.26 105.03  873.54
## 8   8 82.998 390.71 157.24  708.77
## 9   9 72.392 318.12 183.55  951.18
## 10 10 74.039 374.79 128.15  753.72
## 11 11 66.047 425.35 131.46 1037.20
## 12 12 74.902 574.92 163.17  728.25