Heller— title: “Phar7383 BE Class 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)
beclass<-read.csv("C:\\Heller\\PHAR7383\\week2\\beclass.csv",stringsAsFactors = F)
beclasssum<-beclass%>%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=beclass,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=beclasssum,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 <-beclass
## 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 | 475.704482 | 18.58979 | 473.5558 | 364.2501 | 737.6779 | 468.323044 |
| R | cmax | 61.942083 | 16.77527 | 63.1250 | 43.8530 | 80.6490 | 61.079445 |
| R | tmax | 1.208333 | 24.14864 | 1.0000 | 1.0000 | 2.0000 | 1.178257 |
| T | aucinf.obs | 471.139481 | 18.60555 | 484.3006 | 328.0943 | 699.5488 | 463.410721 |
| T | cmax | 43.300833 | 17.03613 | 41.9880 | 32.4580 | 59.2580 | 42.716541 |
| T | tmax | 2.312500 | 55.14758 | 1.7500 | 1.0000 | 6.0000 | 2.059767 |
#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 22 26589.060 <.0001
## TRT 1 22 1.005 0.3271
## PERIOD 1 22 0.687 0.4160
## SEQ 1 22 0.051 0.8226alpha<-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: 97.18–100.75%, PE: 98.95%#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 22 14059.080 <.0001
## TRT 1 22 234.751 <.0001
## PERIOD 1 22 2.763 0.1107
## SEQ 1 22 0.215 0.6471alpha<-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: 67.19– 72.80%, PE: 69.94%