phar7383_bql_student_Heller

Heller date: “2025-09-20” title: “PHAR7383 Missing data” output: html_document —

#libraries

library(dplyr)
library(ggplot2)
library(xpose4)
library(tidyr)
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"))}

#import data

#unt456<-read.csv("/work/ac0837/PHAR7383/wk9-missingdata/unt456.csv",stringsAsFactors = F)
#unt234<-read.csv("C:\\Heller\\PHAR7383\\week8\\nmdtunt234.csv",stringsAsFactors = F)
unt456<-read.csv("C:\\Heller\\PHAR7383\\week9\\unt456.csv",stringsAsFactors = F)
#percentage BQL 
length(unt456$DV[unt456$BQL==1])*100/length(unt456$DV[unt456$EVID==0])

## [1] 22.10526

#Exploratory data analysis

#linear plot mean versus dose
ggplot(data=unt456%>%filter(EVID==0),aes(TIME,DV))+
geom_point(size=2)+
scale_x_continuous(limits = c(0,24),breaks = c(0,2,4,8,12,24))+
theme_bw()+
my_theme()+
labs(x="Time after dose (hour)",y="Plasma concentration (ng/ml)")

## 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.

#Twocomp M1 method Diagnostic plots (not using xpose4, just to remind you. If you want to use xpose4 for these plots, you need to include required table files in model file eg. sdtab,cotab etc.)

#twocompm1<-read.table("/work/ac0837/PHAR7383/wk9-missingdata/twocompm1.res",skip=1,header = T)
#unt456<-read.csv("C:\\Heller\\PHAR7383\\week9\\unt456.csv",stringsAsFactors = F)
twocompm1<-read.table("C:\\Heller\\PHAR7383\\week9\\twocompm1.res",skip=1,header = T)  
#observed and individual predictions
ggplot(data=twocompm1%>%filter(TIME>0),aes(DV,IPRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,20))+
  scale_y_continuous(limits = c(0,20))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Individual Predictions (ng/ml)")

#observed and population predictions
ggplot(data=twocompm1%>%filter(TIME>0),aes(DV,PRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,20))+
  scale_y_continuous(limits = c(0,20))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Population Predictions (ng/ml)")

#TIME versus CRES plot
ggplot(data=twocompm1%>%filter(TIME>0),aes(TIME,CWRES))+
  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 (hour)",y="Conditional Weighted residuals (CWRES)")

#Predictions versus CWRES plot
ggplot(data=twocompm1%>%filter(TIME>0),aes(PRED,CWRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Population Predictions (ug/L)",y="Weighted residuals (CWRES)")

#Twocomp M3 method Diagnostic plots (not using xpose4, just to remind you. If you want to use xpose4 for these plots, you need to include required table files in model file eg. sdtab,cotab etc.)

#twocompm3<-read.table("/work/ac0837/PHAR7383/wk9-missingdata/twocompm3.res",skip=1,header = T)
twocompm3<-read.table("C:\\Heller\\PHAR7383\\week9\\twocompm3.res",skip=1,header = T) 

#note that BQL==0 is set for plots.  This is to only use non-BQL data for plots.  The predictons for BQL data is not concentration for probability. So, we do not plot.  More advanced plots are avialble in PSN but we will visit them in future.
#observed and individual predictions
ggplot(data=twocompm3%>%filter(TIME>0&BQL==0),aes(DV,IPRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,20))+
  scale_y_continuous(limits = c(0,20))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Individual Predictions (ng/ml)")

#observed and population predictions
ggplot(data=twocompm3%>%filter(TIME>0&BQL==0),aes(DV,PRED))+
  geom_point(shape=19,colour="blue")+
  geom_abline(intercept =0, slope =1)+
  scale_x_continuous(limits = c(0,20))+
  scale_y_continuous(limits = c(0,20))+
  theme_bw()+
  my_theme()+
labs(x="Observations (ng/ml)",y="Population Predictions (ng/ml)")

#TIME versus CRES plot
ggplot(data=twocompm3%>%filter(TIME>0&BQL==0),aes(TIME,CWRES))+
  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 (hour)",y="Conditional Weighted residuals (CWRES)")

#Predictions versus CWRES plot
ggplot(data=twocompm3%>%filter(TIME>0&BQL==0),aes(PRED,CWRES))+
  geom_point(shape=19) +
  geom_hline(yintercept = 0, colour="black") +
  geom_hline(yintercept = c(-6,6),linetype = 2)+
  theme_bw()+
  my_theme()+
  labs(x=" Population Predictions (ug/L)",y="Weighted residuals (CWRES)")

#Twocomp M1 method Diagnostic plots using xpose4, including required table files in model file eg. sdtab,cotab etc.)

run1<-xpose.data(1,dir="/Heller/PHAR7383/week9")

## 
## Looking for NONMEM table files.
##     Reading /Heller/PHAR7383/week9/sdtab1 
##     Reading /Heller/PHAR7383/week9/patab1 
##     Reading /Heller/PHAR7383/week9/catab1 
##     Reading /Heller/PHAR7383/week9/cotab1 
## Table files read.
##     Reading /Heller/PHAR7383/week9/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)

#Twocomp M3 method Diagnostic plots using xpose4, including required table files in model file eg. sdtab,cotab etc.)

run2<-xpose.data(2,dir="/Heller/PHAR7383/week9")

## 
## Looking for NONMEM table files.
##     Reading /Heller/PHAR7383/week9/sdtab2 
##     Reading /Heller/PHAR7383/week9/patab2 
##     Reading /Heller/PHAR7383/week9/catab2 
##     Reading /Heller/PHAR7383/week9/cotab2 
## Table files read.
##     Reading /Heller/PHAR7383/week9/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)