Heller Phar7383 Wk4 02_15_25 title: “Individual approach oral case” output: html_document

---

#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\\week4\\oralcase\\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)")

#Noncompartmental analysis

## Load the PK concentration data
pkdt <-wk4_oral
## Generate the dosing data
dose <- pkdt%>%filter(TIME==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(
    pkdt,
    CONC~TIME|SEQ+TRT+PERIOD+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(
    dose,
    DOSE~TIME|SEQ+TRT+PERIOD+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)

Processing of NCA results

#extracting data frame with parameters
pknca_res<-results_obj$result
#creating subset of data for be analysis
beanalysis<-pknca_res%>%group_by(ID,SEQ,TRT,PERIOD,PPTESTCD)%>%filter(PPTESTCD=="aucinf.obs"|PPTESTCD=="cmax"|PPTESTCD=="tmax")%>%
  select(-start,-end,-exclude)
#summary of pk parameters
beanalysissum<-beanalysis%>%group_by(TRT,PPTESTCD)%>%filter(PPTESTCD=="aucinf.obs"|PPTESTCD=="cmax"|PPTESTCD=="tmax")%>%
  summarise(Average=mean(PPORRES),CV=sd(PPORRES)*100/mean(PPORRES),Median=median(PPORRES),mini=min(PPORRES),maxi=max(PPORRES),gmean=exp(mean(log(PPORRES))))

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

beanalysissum%>%kable

TRT	PPTESTCD	Average	CV	Median	mini	maxi	gmean
R	aucinf.obs	459.901763	13.08981	462.4901	362.8237	546.7256	456.191384
R	cmax	53.688167	14.46584	54.1360	41.5700	65.9450	53.173527
R	tmax	1.666667	23.35497	1.7500	1.0000	2.0000	1.618870
T	aucinf.obs	482.975337	13.50622	482.8950	375.6105	610.6449	478.922860
T	cmax	65.037083	13.98238	63.9500	46.8580	80.9440	64.427959
T	tmax	1.500000	20.10076	1.5000	1.0000	2.0000	1.470841

#BE Analysis- AUC

be_auc<-beanalysis%>%filter(PPTESTCD=="aucinf.obs")
library(nlme)

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

mod <- lme(log(PPORRES) ~ TRT+PERIOD+SEQ,random=~1|ID, data = be_auc)
pe  <- 100*exp(summary(mod)$tTable["TRTT", "Value"])
ci  <- 100*exp(intervals(mod, which = "fixed",
                         level = 1-2*0.05)[[1]]["TRTT", c(1, 3)])
print(anova(mod))

##             numDF denDF   F-value p-value
## (Intercept)     1    10 24750.431  <.0001
## TRT             1    10    14.008  0.0038
## PERIOD          1    10     3.543  0.0892
## SEQ             1    10     0.577  0.4648

alpha<-0.05
cat(paste0("\n",100*(1-2*alpha), "% CI:",
           sprintf("%6.2f%s%6.2f%%", ci[1], "\u2013", ci[2]),
           ", PE:", sprintf("%6.2f%%", pe), "\n")) # rounded acc. to GLs

## 
## 90% CI:102.54–107.48%, PE:104.98%

#BE Analysis- Cmax

be_cmax<-beanalysis%>%filter(PPTESTCD=="cmax")
library(nlme)
mod <- lme(log(PPORRES) ~ TRT+PERIOD+SEQ,random=~1|ID, data = be_cmax)
pe  <- 100*exp(summary(mod)$tTable["TRTT", "Value"])
ci  <- 100*exp(intervals(mod, which = "fixed",
                         level = 1-2*0.05)[[1]]["TRTT", c(1, 3)])
print(anova(mod))

##             numDF denDF  F-value p-value
## (Intercept)     1    10 9391.033  <.0001
## TRT             1    10   69.901  <.0001
## PERIOD          1    10    2.877  0.1207
## SEQ             1    10    0.011  0.9183

alpha<-0.05
cat(paste0("\n",100*(1-2*alpha), "% CI:",
           sprintf("%6.2f%s%6.2f%%", ci[1], "\u2013", ci[2]),
           ", PE:", sprintf("%6.2f%%", pe), "\n")) # rounded according to guidelines

## 
## 90% CI:116.23–126.31%, PE:121.17%

#preparation of NONMEM dataset

#nmdtoral<-wk4_oral%>%
  #filter(TRT=="T")%>%
  #select(-SEQ,-TRT,-PERIOD)%>%
  #mutate(AMT=ifelse(TIME==0,DOSE,0))%>%
  #mutate(EVID=ifelse(TIME==0,1,0))%>%
 # rename(CID=ID)
#write.csv(nmdtoral,"/work/ac0837/PHAR7383/wk4-indapproach/Oralcase/nmdtoral.csv",row.names = F,quote=F)

wk4_oral<-read.csv(“C:\Heller\PHAR7383\week4\oralcase\oralindapproach.csv”,stringsAsFactors = F)

#model building #batch processing

onecompres<-list()
for (i in 1:12){
onecompres[[i]]<-read.table(paste ("/Heller/PHAR7383/week4/oralcase/id",           i,"onecmp.res",sep=""),skip=1,header=T)}
onecompresdf<-do.call("rbind",onecompres)%>%
  mutate(model="onecomp")

twocompres<-list()
for (i in 1:12){
twocompres[[i]]<-read.table(paste("/Heller/PHAR7383/week4/oralcase/id",i,"twocmp.res",sep=""),skip=1,header=T)}
twocompresdf<-do.call("rbind",twocompres)%>%
  mutate(model="twocomp")

modelcomb<-rbind(onecompresdf,twocompresdf)

#Plotting ID 1

#Time versus oberved and predicted
ggplot(data=modelcomb%>%filter(ID==1),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==1),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==1),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 3

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==3),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==3),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==3),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 4

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==4),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==4),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==4),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 5

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==5),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==5),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==5),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 7

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==7),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==7),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==7),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 10

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==10),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==10),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==10),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 11

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==11),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==11),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==11),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Plotting ID 12

#Time versus observed and predicted
ggplot(data=modelcomb%>%filter(ID==12),aes(TIME,CONC))+
  geom_point(size=3) +
  geom_line(aes(TIME,IPRED),linetype="dashed")+
  theme_bw()+
  my_theme()+
  labs(x="Time (hour)",y="Concentration (ug/L)")+
  facet_wrap(vars(model))

#TIME versus WRES plot
ggplot(data=modelcomb%>%filter(ID==12),aes(TIME,WRES))+
  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 (min)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))

#Predictions versus WRES plot
ggplot(data=modelcomb%>%filter(ID==12),aes(IPRED,WRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Predictions (ug/L)",y="Weighted residuals (WRES)")+
  facet_wrap(vars(model))