Heller— title: “Phar7383 BE Homework Ex Wk2 01_29_25” 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, unionlibrary(PKNCA)## Warning: package 'PKNCA' was built under R version 4.4.2##
## Attaching package: 'PKNCA'## The following object is masked from 'package:stats':
##
## filterlibrary(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 (student starts from here)
behome<-read.csv("C:\\Heller\\PHAR7383\\week2a\\behomework.csv",stringsAsFactors = F)
behomesum<-behome%>%group_by(TRT,TIME)%>%summarise(cmean=mean(CONC),stdev=sd(CONC))## `summarise()` has grouped output by 'TRT'. You can override using the `.groups`
## argument.#plot
#Population Plot
ggplot(data=behome,aes(TIME,CONC,group=ID ))+
geom_line(size=0.5)+
scale_x_continuous(limits = c(0,72),breaks = c(0,1,2,4,6,8,12,24,48,72))+
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.
#average plot with standard deviation (+/-1SD)
ggplot(data=behomesum,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 <-behome
## 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 | 460.458306 | 13.70656 | 464.6776 | 357.8251 | 546.0039 | 456.36751 |
| R | cmax | 47.734667 | 14.04817 | 47.7765 | 36.8360 | 58.4150 | 47.30221 |
| R | tmax | 1.666667 | 23.35497 | 1.7500 | 1.0000 | 2.0000 | 1.61887 |
| T | aucinf.obs | 470.320293 | 13.70070 | 475.5687 | 357.6764 | 587.8078 | 466.18995 |
| T | cmax | 43.019167 | 13.11335 | 42.3880 | 32.7880 | 52.8420 | 42.67583 |
| T | tmax | 2.333333 | 64.17740 | 1.5000 | 1.5000 | 6.0000 | 2.03081 |
#BE Analysis- AUC
be_auc<-beanalysis%>%filter(PPTESTCD=="aucinf.obs")
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapsemod <- 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 22475.910 <.0001
## TRT 1 10 2.363 0.1553
## PERIOD 1 10 0.179 0.6811
## SEQ 1 10 0.492 0.4991alpha<-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: 99.62–104.75%, PE:102.15%#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 9645.008 <.0001
## TRT 1 10 13.922 0.0039
## PERIOD 1 10 1.034 0.3332
## SEQ 1 10 0.090 0.7701alpha<-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: 85.82– 94.84%, PE: 90.22%