#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_oral<-read.csv("C:\\Heller\\PHAR7383\\week5\\oralindapproach.csv",stringsAsFactors = F)
wk4_oral_sum<-wk4_oral%>%group_by(TRT,TIME)%>%summarise(cmean=mean(CONC),stdev=sd(CONC))
## `summarise()` has grouped output by 'TRT'. You can override using the `.groups`
## argument.

#Exploratory plots

#Population Plot
ggplot(data=wk4_oral,aes(TIME,CONC,group=ID))+
geom_line(size=0.5)+
geom_point(size=1)+
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(TRT))
## 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 Reference
ggplot(data=wk4_oral%>%filter(TRT=="R"),aes(TIME,CONC))+
geom_line(size=0.5)+
geom_point(size=1)+
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(ID))

#Individual Plot Test
ggplot(data=wk4_oral%>%filter(TRT=="T"),aes(TIME,CONC))+
geom_line(size=0.5)+
geom_point(size=1)+
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(ID))

#average plot with standard deviation (+/-1SD)
ggplot(data=wk4_oral_sum,aes(TIME,cmean,color=TRT))+
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 (hour)",y="Plasma concentration (ng/ml)")

#model building`

#Run1- One comp model, first order absorption; BSV on KA; ADD+PROP RUV NO INTER

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")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab1 
##     Reading C:\Heller\PHAR7383\week5/patab1 
##     Reading C:\Heller\PHAR7383\week5/catab1 
##     Reading C:\Heller\PHAR7383\week5/cotab1 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/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/patab1",header=T,skip=1)%>%
  select(ID,KA,CL,V)%>%
  filter(!duplicated(ID))
print (run1_ebe)
##    ID     KA     CL      V
## 1   1 1.8518 100.82 749.73
## 2   2 2.3932 100.82 749.73
## 3   3 1.9159 100.82 749.73
## 4   4 2.6133 100.82 749.73
## 5   5 2.5849 100.82 749.73
## 6   6 1.8517 100.82 749.73
## 7   7 1.9008 100.82 749.73
## 8   8 2.3511 100.82 749.73
## 9   9 2.4453 100.82 749.73
## 10 10 2.4132 100.82 749.73
## 11 11 3.4102 100.82 749.73
## 12 12 2.8686 100.82 749.73

#Run2- One comp model, first order absorption; BSV on KA,CL; ADD+PROP RUV

#Diagnostic plots
run2<-xpose.data(2,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab2 
##     Reading C:\Heller\PHAR7383\week5/patab2 
##     Reading C:\Heller\PHAR7383\week5/catab2 
##     Reading C:\Heller\PHAR7383\week5/cotab2 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/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/patab2",header=T,skip=1)%>%
  select(ID,KA,CL,V)%>%
  filter(!duplicated(ID))
print (run2_ebe)
##    ID     KA      CL      V
## 1   1 1.3196 120.420 736.52
## 2   2 2.1966 107.410 736.52
## 3   3 1.5981  97.003 736.52
## 4   4 2.6114  97.267 736.52
## 5   5 2.4658 106.760 736.52
## 6   6 1.7815  85.673 736.52
## 7   7 1.5682 101.800 736.52
## 8   8 2.3172  89.791 736.52
## 9   9 2.2387 120.270 736.52
## 10 10 2.2906  96.903 736.52
## 11 11 3.5341 111.490 736.52
## 12 12 2.9021  99.561 736.52

#Run3- One comp model, first order absorption; BSV on KA,CL,V; ADD+PROP RUV

#Diagnostic plots
run3<-xpose.data(3,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab3 
##     Reading C:\Heller\PHAR7383\week5/patab3 
##     Reading C:\Heller\PHAR7383\week5/catab3 
##     Reading C:\Heller\PHAR7383\week5/cotab3 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/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/patab3",header=T,skip=1)%>%
  select(ID,KA,CL,V)%>%
  filter(!duplicated(ID))
print (run3_ebe)
##    ID     KA      CL      V
## 1   1 1.7154 123.480 817.55
## 2   2 2.2589 107.590 743.68
## 3   3 2.0537 101.510 849.09
## 4   4 2.4834  96.275 715.14
## 5   5 2.3226 105.010 697.73
## 6   6 1.8869  84.639 726.02
## 7   7 1.8673 103.040 785.03
## 8   8 2.3493  89.467 741.60
## 9   9 2.3174 121.930 755.79
## 10 10 2.4049  97.931 767.49
## 11 11 2.9403 107.630 656.29
## 12 12 2.6563  97.932 701.07

#Run5- Two comp model, first order absortion; BSV on KA;ADD PROP RUV NO INTER

#Diagnostic plots
run5<-xpose.data(5,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab5 
##     Reading C:\Heller\PHAR7383\week5/patab5 
##     Reading C:\Heller\PHAR7383\week5/catab5 
##     Reading C:\Heller\PHAR7383\week5/cotab5 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/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/patab5",header=T,skip=1)%>%
  select(ID,KA,CL,V2,Q,V3)%>%
  filter(!duplicated(ID))
print (run5_ebe)
##    ID      KA     CL     V2      Q     V3
## 1   1 0.58699 94.625 322.18 177.66 476.01
## 2   2 0.83948 94.625 322.18 177.66 476.01
## 3   3 0.64244 94.625 322.18 177.66 476.01
## 4   4 0.95404 94.625 322.18 177.66 476.01
## 5   5 0.94918 94.625 322.18 177.66 476.01
## 6   6 0.70318 94.625 322.18 177.66 476.01
## 7   7 0.65105 94.625 322.18 177.66 476.01
## 8   8 0.85331 94.625 322.18 177.66 476.01
## 9   9 0.85067 94.625 322.18 177.66 476.01
## 10 10 0.85336 94.625 322.18 177.66 476.01
## 11 11 1.38060 94.625 322.18 177.66 476.01
## 12 12 1.13580 94.625 322.18 177.66 476.01

#Run6- Two comp model, first order absortion; BSV on KA,CL;ADD PROP RUV WITH INTER

run6<-xpose.data(6,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab6 
##     Reading C:\Heller\PHAR7383\week5/patab6 
##     Reading C:\Heller\PHAR7383\week5/catab6 
##     Reading C:\Heller\PHAR7383\week5/cotab6 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/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/patab6",header=T,skip=1)%>%
  select(ID,KA,CL,V2,Q,V3)%>%
  filter(!duplicated(ID))
print (run6_ebe)
##    ID      KA      CL     V2      Q     V3
## 1   1 0.57068 122.810 331.97 179.83 434.63
## 2   2 0.85995 101.750 331.97 179.83 434.63
## 3   3 0.65056 104.590 331.97 179.83 434.63
## 4   4 1.00470  90.741 331.97 179.83 434.63
## 5   5 0.98283  96.980 331.97 179.83 434.63
## 6   6 0.73191  74.228 331.97 179.83 434.63
## 7   7 0.65183  97.432 331.97 179.83 434.63
## 8   8 0.88453  81.389 331.97 179.83 434.63
## 9   9 0.89101 123.960 331.97 179.83 434.63
## 10 10 0.87151  93.811 331.97 179.83 434.63
## 11 11 1.42330 101.630 331.97 179.83 434.63
## 12 12 1.19310  96.491 331.97 179.83 434.63

#Run7- Two comp model, first order absortion; BSV on KA,CL,V2;ADD PROP RUV WITH INTER

run7<-xpose.data(7,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab7 
##     Reading C:\Heller\PHAR7383\week5/patab7 
##     Reading C:\Heller\PHAR7383\week5/catab7 
##     Reading C:\Heller\PHAR7383\week5/cotab7 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/run7.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run7,type="p")

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

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

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

ind.plots(run7)

ranpar.hist(run7)

#individual parameters
run7_ebe<-read.table("C:\\Heller\\PHAR7383\\week5/patab7",header=T,skip=1)%>%
  select(ID,KA,CL,V2,Q,V3)%>%
  filter(!duplicated(ID))
print (run7_ebe)
##    ID      KA      CL     V2      Q     V3
## 1   1 0.72667 121.090 430.58 177.31 414.89
## 2   2 0.97145 101.650 373.48 177.31 414.89
## 3   3 0.84484 102.950 432.21 177.31 414.89
## 4   4 0.99769  91.506 327.30 177.31 414.89
## 5   5 1.00820  97.474 339.23 177.31 414.89
## 6   6 0.77062  74.555 345.44 177.31 414.89
## 7   7 0.76774  96.725 391.49 177.31 414.89
## 8   8 1.00350  81.296 373.09 177.31 414.89
## 9   9 0.97300 124.080 363.68 177.31 414.89
## 10 10 1.02430  93.496 386.13 177.31 414.89
## 11 11 1.31210 102.490 309.13 177.31 414.89
## 12 12 1.17470  97.130 327.74 177.31 414.89

#Run8- Two comp model, first order absortion; BSV on KA,CL,V2,Q;ADD PROP RUV WITH INTER

run8<-xpose.data(8,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab8 
##     Reading C:\Heller\PHAR7383\week5/patab8 
##     Reading C:\Heller\PHAR7383\week5/catab8 
##     Reading C:\Heller\PHAR7383\week5/cotab8 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/run8.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run8,type="p")

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

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

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

ind.plots(run8)

ranpar.hist(run8)

#individual parameters
run8_ebe<-read.table("C:\\Heller\\PHAR7383\\week5/patab8",header=T,skip=1)%>%
  select(ID,KA,CL,V2,Q,V3)%>%
  filter(!duplicated(ID))
print (run8_ebe)
##    ID      KA      CL     V2      Q     V3
## 1   1 0.72627 121.100 430.19 177.52 414.85
## 2   2 0.97244 101.690 373.30 178.21 414.85
## 3   3 0.84458 102.970 431.79 177.77 414.85
## 4   4 0.99759  91.452 327.95 176.03 414.85
## 5   5 1.00840  97.438 339.70 176.51 414.85
## 6   6 0.77211  74.605 345.24 178.76 414.85
## 7   7 0.76799  96.749 391.27 177.90 414.85
## 8   8 1.00530  81.408 372.23 180.14 414.85
## 9   9 0.97162 123.970 364.39 175.15 414.85
## 10 10 1.02620  93.606 385.40 179.71 414.85
## 11 11 1.31370 102.460 309.67 176.83 414.85
## 12 12 1.17390  97.023 328.87 174.54 414.85

#Run9- Two comp model, first order absortion; BSV on KA,CL,V2,Q,V3;ADD PROP RUV WITH INTER

run9<-xpose.data(9,dir="C:\\Heller\\PHAR7383\\week5")
## 
## Looking for NONMEM table files.
##     Reading C:\Heller\PHAR7383\week5/sdtab9 
##     Reading C:\Heller\PHAR7383\week5/patab9 
##     Reading C:\Heller\PHAR7383\week5/cotab9 
## Table files read.
##     Reading C:\Heller\PHAR7383\week5/run9.phi 
## 
## Looking for NONMEM simulation table files.
## No simulated table files read.
dv.vs.ipred(run9,type="p")

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

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

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

ind.plots(run9)

ranpar.hist(run9)

#individual parameters
run9_ebe<-read.table("C:\\Heller\\PHAR7383\\week5/patab9",header=T,skip=1)%>%
  select(ID,KA,CL,V2,Q,V3)%>%
  filter(!duplicated(ID))
print (run9_ebe)
##    ID      KA      CL     V2      Q     V3
## 1   1 0.72652 120.820 430.47 177.02 418.66
## 2   2 0.97341 101.590 373.86 177.67 416.14
## 3   3 0.84193 102.610 430.52 177.30 423.51
## 4   4 1.00180  91.643 329.37 175.70 409.43
## 5   5 1.01110  97.496 340.74 176.11 412.28
## 6   6 0.77950  75.163 349.12 178.10 398.69
## 7   7 0.77022  96.835 392.66 177.37 411.45
## 8   8 1.00850  81.474 373.63 179.40 412.35
## 9   9 0.97138 123.600 364.45 174.85 420.77
## 10 10 1.02650  93.493 385.76 179.07 417.03
## 11 11 1.31500 102.360 310.12 176.42 416.58
## 12 12 1.17180  96.719 328.23 174.35 423.34