R/NONMEM_plot.R
In jgrevel/BAST1-R-Library: BAST R function library.
Defines functions
NONMEM_plot
Documented in
NONMEM_plot
#' Plots data from NONMEM files
#'
#' The tool provides a good overview of the available data for each patient and helps to identify severe outliers or programming errors. This can be useful when
#' doing QC work upon building a new NONMEM data file. The tool was created to visualise any NONMEM data file that contains continuous type DV. The observations
#' can be PK or PD type data, and dosing records don't necessarily have to be present. To identify potential outliers, observations for each individual can be
#' superimposed on observations for all patients.
#'
#' @author Katarzyna Nurzynska
#'
#' @param data data frame containing data in NONMEM format.
#' @param save_path path to the folder where plots will be saved.
#' @param ID_col name of the column in NONMEM file which contains patient identifiers.
#' @param ind_var_col name of the column in NONMEM file which contains the variable for use as the x-axis in plots, usually 'TIME' or 'TAD'.
#' @param DV_col name of the column in NONMEM file which contains observations for plotting on y-axis. If dose corrected DV is to be used then set DV_col='DV_DOSE_CORRECTED'.
#' @param x_unit label for units on x-axis.
#' @param y_unit label for units on y-axis.
#' @param dose_marker logical. If TRUE, add dose marker at the top of the plot. If the AMT or DOSE column is present in the data file this marker will be coloured
#' according to the dose size. A separate plot Legend.png will be generated with a key to colouring of different dose sizes.
#' @param panel logical. If FALSE, generate single plots for each ID, or if TRUE, panels of 8 plots.
#' @param keep_doses_after_last_obs logical. If TRUE, the x-axis of each plot will be scaled such that doses after the last observation are visualised. If FALSE, the
#' x-axis will only reach as far as the last observation.
#' @param calc_dose logical. If TRUE, add DOSE column (showing size of most recent non-zero AMT) to the data. Required if DOSE is not already in the data file
#' and you set either dose_marker=TRUE or DV_col=DV_DOSE_CORRECTED.
#' @param calc_tad logical. If TRUE, add TAD column to the data. Required if TAD is not already in the data file and ind_var_col='TAD'.
#' @param calc_occ logical. If TRUE, add OCC column to the data. Required if OCC is not already in the data file and you want to plot data from different occasions using
#' occ=1 or occ=2 etc.
#' @param occ numeric. Plot data for a particular occasion. Default is FALSE which leads to plots including data from all occasions. If the argument is equal to e.g. 2
#' only data from second occasion is plotted.
#' @param LLOQFLAG_col name of the LLOQFLAG column in NONMEM data file, which indicates if an observation is BLQ (1) or has been quantified (0). BLQ observations are
#' plotted as green circles. Set to FALSE if a LLOQFLAG column is not in the data file.
#' @param LLOQ_col name of the LLOQ column in NONMEM data file, which contains the values of the limit of quantification. If the argument is a number (rather than text
#' such as 'LLOQ') then the lower limit of quantification for every sample will be set to that number, and the BLQ samples will be plotted as green circles at that
#' concentration.
#' @param title_label_col vector of character strings, containing names of the columns in NONMEM data file which will be used for supplying some additional information
#' to plot titles. If title_label_col='' then no additional labels are added. For multiple dose studies, the dose shown in the title will correspond to the selected
#' occasion. If all occasions are plotted, the value of the dose in the title will correspond to the first occasion.
#' @param colour_shading numeric. Sets the relationship between the dose amounts and the colour shades. Important when the range of doses is wide. Try experimenting
#' with different values to see effect on Legend.png and the dose markers. It should be a positive number. If set to 0 then it will instead be calculated as
#' \eqn{sqrt(max(AMT)/min(AMT))/2}.
#' @param colour_range vector of colour names. Sets the colour transitions for dose markers and background data. A minimum of 2 colours is required and plot colours will
#' be interpolated between the members of colour_range, depending on the size of DOSE.
#' @param backround_data logical. If TRUE, background data (data from all other individuals) are added to the plot.
#' @param plot_IPRED logical. If TRUE, then IPRED data are added to the plots. Obviously the data file will have to be a NONMEM $TABLE file such that IPRED data are
#' available.
#' @param all_same_scale_x logical. If TRUE, x-axis range should be equal for all individuals. If FALSE, x-axis range should be different for each individual.
#' @param all_same_scale_y logical. If TRUE, y-axis range should be equal for all individuals. If FALSE, y-axis range should be different for each individual.
#' @param file_text character. Free text that will be added to each filename. Useful when creating lots of different visualisations from the same data.
#' @param ind_var_conversion numeric. Factor to be used for converting the data in the x-axis. For example set to 1/24 to convert TIME from hours to days.
#' @param extend_dose_marker logical. If TRUE, the dose markers will be extended as grey dashed lines from top to bottom of each plot. This can help determine if
#' observations close in time to a dose actually lie before or after the dosing event.
#'
#' @return The tool will plot observations and dosing records vs time, saving one plot for each individual or a panel of eight plots for subject-rich data files.
#' @details User should consult the latest Best Practice guide at
#'
#' \code{\\\\BAST_SERVER\\BAST_Shared_Files\\QA Documents\\BAST Computer Systems Validation\\final_documents\\Procedure_for_visualizing_NONMEM_files\\}
#'
#' for more in-depth documentation on the use of this function.
#' @note Only records with EVID=0 are plotted.
#' @note Function will accept standard arguments for the default R plot function, such as log='y'.
#' @note Columns 'ID', 'EVID', 'DV' and 'TIME' are expected to be present. Columns 'TAD', 'AMT', DOSE' and 'OCC' (occasion) are optional.
#' @note It is assumed that DV and IPRED are on a linear scale. If not, then the user should correct this before using the NONMEM_plot function.
#' @note Before running the function, the user is expected to filter data file to records of interest e.g. placebo or biomarker type, through filtering based
#' on columns such as CMT or DVID.
#' @export
############################################################################################################################################################################################################################################################################################################################################
#FUNCTION: NONMEM_visualising_tool
#PURPOSE: The tool provides a good overview of the available data for each patient and helps to identify severe outliers or programming errors.
# This can be useful when doing QC work upon building a new NONMEM data file.
# The tool was created to visualise any NONMEM data file that contains continuous type DV.
# The observations could be PK or PD type data, and dosing records don't necessarily have to be present.
# To identify potential outliers, observations for each individual can be superimposed on observations for all patients.
# The tool will plot observations and dosing records vs time, saving one plot for each individual or a panel of eight plots for subject-rich data files.
############################################################################################################################################################################################################################################################################################################################################
#
# Only records with EVID=0 are plotted
# Function will accept standard arguments for the default R plot command, such as log='y'
# Columns 'ID', 'EVID', 'DV' and 'TIME' are expected to be present. Columns 'TAD', 'AMT', DOSE' and 'OCC' (occasion) are optional.
# 'TAD', 'DOSE' and 'OCC' can be added to the data file before plotting by using the following arguments: calc_tad=TRUE, calc_dose=TRUE, calc_occ=TRUE.
# It is assumed that DV and IPRED are on a linear scale. If not, then the user should correct this before using the NONMEM_plot function
# Before running the function, the user is expected to filter data file to records of interest e.g. placebo or biomarker type, through filtering based on columns such as CMT or DVID
#
# Minimum function inputs:
# NAME DEFAULT VALUE DESCRIPTION
# data NONE name of the data frame containing the NONMEM-formatted data
# save_path NONE path to the folder where plots will be saved
#
# Optional function inputs:
# NAME DEFAULT VALUE DESCRIPTION
# ID_col 'ID' Name of the column in NONMEM file which contains patient identifiers.
# ind_var_col 'TIME' Name of the column in NONMEM file which contains the variable for use as the x-axis in plots. 'TAD' is the other main alternative for ind_var_col.
# DV_col 'DV' Name of the column in NONMEM file which contains observations for plotting on y-axis. If dose corrected DV is to be used then set DV_col='DV_DOSE_CORRECTED'.
# x_unit 'hours' Label for units on x-axis.
# y_unit 'ug/L' Label for units on y-axis.
# dose_marker TRUE Add dose marker at the top of the plot. If the AMT or DOSE column is present in the data file this marker will be coloured according to the dose size. A separate plot Legend.png will be generated with a key to colourin g of different dose sizes.
# panel TRUE Generate single plots for each ID when the argument is set to FALSE, or panels of 8 plots when the argument is set to TRUE.
# keep_doses_after_last_obs TRUE When set to TRUE the x-axis of each plot will be scaled such that doses after the last observation are visualised. When set to FALSE the x-xis will only reach as far as the last observation
# calc_dose TRUE Add DOSE column (showing size of most recent non-zero AMT) to the data when the argument is TRUE. Required if DOSE is not already in the data file and you set either dose_marker=TRUE or DV_col=DV_DOSE_CORRECTED.
# calc_tad TRUE Add TAD column if the argument is equal to TRUE. Required if TAD is not already in the data file and ind_var_col='TAD'
# calc_occ TRUE Add OCC column if the argument is equal to TRUE. Required if OCC is not already in the data file and you want to plot data from different occasions using occ=1 or occ=2 etc.
# occ FALSE Plot data for a particular occasion. Default is FALSE which leads to plots including data from all occasions. If the argument is equal to e.g. 2 only data from second occasion is plotted.
# LLOQFLAG_col FALSE Name of the column in NONMEM file which indicates if an observation is BLQ (1) or has been quantified (0). BLQ observations are plotted as green circles. Set to FALSE if a LLOQFLAG column is not in the data file.
# LLOQ_col FALSE Name of the column in NONMEM file which contains the values of the limit of quantification. If the argument is a number (rather than text such as 'LLOQ') then the lower limit of quantification for every sample will be set to that number, and the BLQ samples will be plotted as green circles at that concentration.
# title_label_col c('DOSE','OCC') Names of the columns in NONMEM file which will be used for supplying some additional information to plot titles. If title_label_col='' then no additional labels are added. For multiple dose studies the dose shown in the title will correspond to the selected occasion. If all occasions are plotted the value of the dose in the title will correspond to the first occasion.
# colour_shading 0 Sets the relationship between the dose amounts and the colour shades. Important when the range of doses is wide. Try experimenting with different values to see effect on Legend.png and the dose markers. It should be a positive number. If set to 0 then it will instead be calculated as sqrt(max(AMT)/min(AMT))/2.
# colour_range c('mistyrose','pink','red','red4','black') Sets the colour transitions for dose markers and background data. A minimum of 2 colours is required and plot colours will be interpolated between the members of colour_range, depending on the size of DOSE.
# backround_data FALSE If TRUE background data (data from all other individuals) are added to the plot.
# plot_IPRED FALSE If set to TRUE then IPRED data are added to the plots. Obviously the data file will have to be a NONMEM $TABLE file such that IPRED data are available.
# all_same_scale_x FALSE Determines if x-axis range should be same for all individuals (TRUE) or different for each individual (FALSE).
# all_same_scale_y TRUE Determines if y-axis range should be the same for all individuals (TRUE) or different for each individual (FALSE).
# file_text '' Free text that will be added to each filename. Useful when creating lots of different visualisations from the same data.
# ind_var_conversion 1 Factor to be used for converting the data in the x-axis. For example set to 1/24 to convert TIME from hours to days.
# extend_dose_marker FALSE When set to TRUE the dose markers will be extended as grey dashed lines from top to bottom of each plot. This can help determine if observations close in to time to a dose actually lie before or after the dosing event.
##############################################################
##############################################################
NONMEM_plot=function(data,save_path,ID_col='ID',ind_var_col='TIME',DV_col='DV',x_unit='hours',y_unit='ug/L',dose_marker=TRUE,panel=TRUE,
keep_doses_after_last_obs=TRUE,calc_dose=TRUE,calc_occ=TRUE,calc_tad=TRUE,occ=FALSE,LLOQFLAG_col=FALSE,LLOQ_col=FALSE,
title_label_col=c('DOSE','OCC'),colour_shading=0,
colour_range=c('mistyrose','pink','red','red4','black'),background_data=FALSE,plot_IPRED=FALSE,all_same_scale_x=FALSE,
all_same_scale_y=TRUE,file_text='',ind_var_conversion=1,extend_dose_marker=FALSE,...){
##############
LOG=FALSE # check if x axis is logged
if(!is.null(list(...)$log))
{
LOG=grepl('x',list(...)$log)
}
plot0=function(x,y,xlimi=NULL,...) ## A function added by Aaron for doing nice looking things with log scale on x axis
{
LOG=FALSE
x0=0
if(!is.null(list(...)$log))
{
LOG=grepl('x',list(...)$log)
}
snapx=function(x)
{
usr=par()$usr
opar <- par('xlog','ylog')
par(xlog=FALSE)
par(ylog=FALSE)
w=0.5*strwidth('E')
h=0.5*strheight('W')
polygon(x+c(w,w,-w,-w),usr[3]+c(h,-h,-h,h),col='white',lty=0,xpd=NA)
lines(x+c(w,w,NA,-w,-w),usr[3]+c(h,-h,NA,-h,h),xpd=NA)
par(opar)
}
if(LOG) # & any(x<=0) # this part has been excluded from use here so results have consistency
{
if(!is.null(xlimi))
{
r=log(xlimi,10)
}else{
r=range(log(x[x>0],10),na.rm=TRUE)
}
x0=r[1]-0.2*diff(r)
xr=x
xr[xr<=0]=10**x0
plot(xr,y,xaxt='n',xlim=10**c(x0,r[2]),...)
ax=axisTicks(r+diff(r)*0.05*c(-1,1),log=TRUE)
axis(1,at=c(10**x0,ax),labels=c(0,ax))
snapx(0.5*r[1]+0.5*x0)
}else{
if(is.null(xlimi)) xlimi=range(x,na.rm=TRUE)
plot(x,y,xlim=xlimi,...)
}
return(invisible(10**x0))
}
points0=function(x,y,x0,...)
{
if(par()$xlog)
{
x[x<=0]=x0
}
points(x,y,...)
}
lines0=function(x,y,x0,...)
{
if(par()$xlog)
{
x[x<=0]=x0
}
lines(x,y,...)
}
##############
exp_data=data
# set axes labels
if(ind_var_col=='TIME'){
x_label=paste0('Time (',x_unit,')')
}
if(ind_var_col=='TAD'){
x_label=paste0('Time after dose (',x_unit,')')
}
if(ind_var_col!='TIME' & ind_var_col!='TAD'){
x_label=paste0(ind_var_col,' (',x_unit,')')
}
y_label=paste0(DV_col,' (',y_unit,')')
# rename some variables
all_subjects_colours=colour_range
filter=file_text
NONMEM_table=plot_IPRED
if(LLOQFLAG_col!=FALSE & LLOQ_col==FALSE){
stop("LLOQFLAG_col has been set to a column name but LLOQ_col has been set to FALSE. When LLOQFLAG_col has a name then LLOQ_col must also have a name")
}
if(LLOQFLAG_col==FALSE){
LLOQ_col=FALSE
}
# if LLOQ_col is a number then it is not a column name but it is an LLOQ value to be applied to all samples where LLOQFLAG==1
LIMIT=NA
if(is.numeric(LLOQ_col)==TRUE){
LIMIT=LLOQ_col
}
# if user has specified a particular values of occ to plot, then it makes no sence to have set keep_doses_after_last_obs to TRUE, so set it to FALSE
if(occ!=FALSE){
keep_doses_after_last_obs=FALSE
}
if(!DV_col%in%c('DV','DV_DOSE_CORRECTED')){
stop('The only permitted values for DV_col are: DV or DV_DOSE_CORRECTED')
}
plot_counter=0
legend_counter=0
#-----------------------------------------------------------------------------------------------------------------
if(calc_dose==TRUE & !"AMT"%in%names(exp_data)){
stop("You have requested calculation of DOSE column but the is no AMT column in the data file. DOSE cannot be calculated")
}
# check if dose column is required
if((dose_marker==TRUE | "DOSE"%in%title_label_col | DV_col=="DV_DOSE_CORRECTED") & !"DOSE"%in%names(exp_data)){
# dose column is required and is not already in the data
if(calc_dose==FALSE){
stop("The function arguments require that a DOSE column available. DOSE is not in the data, so you much set calc_dose to TRUE to continue")
}
if(calc_dose==TRUE & !"AMT"%in%names(exp_data)){
stop("You have requested calculation of DOSE column but the is no AMT column in the data file. DOSE cannot be calculated")
}
}
#---------------------------------------------------- Option to add DOSE column to the database ------------------------------
if (calc_dose==TRUE){
exp_data=calc_dose(exp_data,ID_col=ID_col)
}
#------------------------------------------------------ DOSE column ----------------------------------
if("DOSE"%in%names(exp_data) & DV_col=="DV"){
colours_column=TRUE
plot_dose_marker=TRUE
}
if("DOSE"%in%names(exp_data) & DV_col=="DV_DOSE_CORRECTED"){
colours_column=TRUE
plot_dose_marker=TRUE
}
if(!"DOSE"%in%names(exp_data) & DV_col=="DV"){ #in this case function cannot add colours to the plot based on the dose size
colours_column=FALSE
COLOUR_occ='gray'
COLOUR_all='gray'
plot_dose_marker=FALSE
}
#---------------------------------------------------- Option to add TAD column --------------------------
if(calc_tad==TRUE & "EVID"%in%names(exp_data) & ID_col%in%names(exp_data) & "TIME"%in%names(exp_data)){
exp_data=calc_tad(exp_data,ID_col=ID_col)
}
if(ind_var_col=="TAD" & !"TAD"%in%names(exp_data)){
stop("Function arguments require that a TAD column is in the data file. TAD is not in the data file. Try setting calc_tad to TRUE")
}
exp_data[,ind_var_col]=exp_data[,ind_var_col]*ind_var_conversion # Convert independent variable. Often this will be conversion of TIME from hours to days
#---------------------------------------------------- Option to add OCC column to the databse -------------------------------
#Each occasion consisits of data records related to administrated dose and all concentrations which were measured afterwards
#Occasion is equal to zero before administration of the first dose
if(occ!=FALSE & calc_occ==FALSE & !"OCC"%in%names(exp_data)){
stop("Function arguments require OCC column to be present in the data file. Try setting calc_occ to TRUE")
}
if(calc_occ==TRUE & !"AMT"%in%names(exp_data)){
stop("Calculation of OCC column has been requested but AMT column is not in the data file. Cannot continue")
}
if(calc_occ==TRUE & "AMT"%in%names(exp_data)){
exp_data=calc_occ(exp_data,ID_col=ID_col,calc_method='by_obs')
}
#---------------------------------------------------- IPRED --------------------------
if (NONMEM_table==TRUE){
IPRED_col="IPRED"
if (DV_col=='DV_DOSE_CORRECTED'){
exp_data$IPRED_DOSE_CORRECTED=exp_data$IPRED/exp_data$DOSE #Calculate dose-corrected DV
IPRED_col="IPRED_DOSE_CORRECTED"
}
}
#------------------------------------------------ Add column with different colours for different doses ---------------------
if(colours_column==TRUE){
dose_unique=unique(exp_data[,"DOSE"])
dose_unique_sorted=sort(dose_unique)
if(colour_shading<0){
colour_shading=0
}
if(colour_shading==0){
colour_shading=sqrt(max(dose_unique)/min(dose_unique))/2
}
#Add column with with different colours for different doses to be later shown on the plot as e.g. background data
colfunc=colorRampPalette(all_subjects_colours, bias=colour_shading) #Data in the backround to be coloured as shades of pink to identify increase in dose
all_dose_colours=colfunc(length(dose_unique_sorted))
exp_data[,"ALL_COLOURS"]='9999'
for(col in 1:length(all_dose_colours)){
exp_data[which(exp_data[,"DOSE"]==dose_unique_sorted[col]),'ALL_COLOURS']=all_dose_colours[col]
}
}
if("OCC"%in%names(exp_data)){
exp_data[,"OCC_ORIG"]=exp_data[,"OCC"] #the original occasion column with OCC=0 remaining unchanged
exp_data[exp_data$OCC==0,"OCC"]=1 #occasion column used for plotting: records before the first dose are plotted as occasion 1
}
#---------------------------------------------------- Calculate dose-corrected DV -------------------------------------------
#----------------------------------------------- Handling data below the limit of quantification ---------------------------
#For data records below the limit of quantification DV_col is equal to limit of quantification; if column with the value of LOQ does not exist then the limit is equal to LIMIT
if(LLOQFLAG_col!=FALSE & DV_col=="DV_DOSE_CORRECTED"){
exp_data$DV_DOSE_CORRECTED=exp_data$DV/exp_data$DOSE #Calculate dose-corrected DV
if(is.na(LIMIT)==FALSE){
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=LIMIT/exp_data[exp_data[,LLOQFLAG_col]==1,"DOSE"] # DV_DOSE_CORRECTED: use the BLQ LIMIT value if it has been entered. So the LIMIT gets dose-corrected
}else{
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=exp_data[exp_data[,LLOQFLAG_col]==1,LLOQ_col]/exp_data[exp_data[,LLOQFLAG_col]==1,"DOSE"] # DV_DOSE_CORRECTED
}
}
if(LLOQFLAG_col!=FALSE & DV_col=="DV"){
if(is.na(LIMIT)==FALSE){
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=LIMIT # DV
}else{
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=exp_data[exp_data[,LLOQFLAG_col]==1,LLOQ_col] # DV
}
}
if(LLOQFLAG_col==FALSE & DV_col=="DV_DOSE_CORRECTED"){
exp_data$DV_DOSE_CORRECTED=exp_data$DV/exp_data$DOSE # DV_DOSE_CORRECTED
}
if(DV_col=="DV_DOSE_CORRECTED"){
# correct any calculations of DV_DOSE_CORRECTED that have problems due to zero dose
exp_data[,DV_col][exp_data$DOSE==0]=exp_data$DV[exp_data$DOSE==0]
}
zero_obs=exp_data[,DV_col][exp_data$EVID==0 & exp_data[,DV_col]<=0]
if(length(zero_obs)>0){
print('WARNING...WARNING...WARNING...')
print(paste0(length(zero_obs),' observations have values less than or equal to zero. This will cause problems if you have requested a logarithmic y axis'))
print('If you are not expecting any observations with value of zero, then check that LLOQFLAG_col and LLOQ_col have been correctly defined')
}
# before TAD and TIME get altered, save the exp_data as all_exp_data, which will be returned by the function at the end
all_exp_data=exp_data
#------------------------- Replace TAD and TIME = 0 with a small value (to avoid problems on a log scale)----------------
### AARON 1 START ### here we make the adjustments to TIME and TAD to enable plotting with log x-axis. The adjustments are actually used with linear x-axis too, which is a bit lazy really
#TIME
#minimum=min(exp_data[exp_data[,"TIME"]!=0,"TIME"]) # we are assuming here that TIME shouldn't normally be negative
#if(minimum>=0.1){
# exp_data[exp_data[,"TIME"]==0,"TIME"]=0.1
#}else{
# exp_data[exp_data[,"TIME"]==0,"TIME"]=minimum
#}
#
#TAD
#if(ind_var_col=="TAD"){
# minimum=min(exp_data[exp_data[,"TAD"]>0,"TAD"]) # remember that TAD could be negative (eg for observations before first dose). For datasets which have records with TAD<0 then log x-axis will be a problem unless the records with negative TAD are filtered before plotting
# if(minimum>=0.1){
# exp_data[exp_data[,"TAD"]==0,"TAD"]=0.1
# }else{
# exp_data[exp_data[,"TAD"]==0,"TAD"]=minimum
# }
#}
### AARON 1 END ###
#-------------------------------------------------------
dose_time=exp_data[exp_data$EVID%in%c(1,4),] #this set of data will be used later to plot vertical lines when dose was administrated
### now check for any subjects who have no observations and remove them
subject=unique(exp_data[,ID_col])
IDS_NULL=NULL
for(g in 1:length(subject)){
obs_recs=exp_data[exp_data[,ID_col]==subject[g] & exp_data$EVID==0,DV_col]
if(length(obs_recs)==0){
IDS_NULL=c(subject[g],IDS_NULL)
}
}
if(length(IDS_NULL)>0){
num_null=length(IDS_NULL)
print(paste0('The following ',num_null,' IDs have no observations and will not be plotted:'))
print(paste0(IDS_NULL,collapse=' '))
exp_data=exp_data[!exp_data[,ID_col]%in%IDS_NULL,]
}
#---------------------------------------------------- Plot data -------------------------------------------------------------
#---------------------------------------------------- Single plot or panel of 8 ---------------------------------------------
subject=unique(exp_data[,ID_col])
print('Generating plots...')
#Plots
if (panel==TRUE){ #panel of 8 plots
loops=ceiling(length(subject)/8) #number of panels each consisiting of 8 plots
N=c(2,4) #2 rows and 4 columns
M=c(4,4,0,0)
plot=8
WIDTH=1500 #width of a figure generated
xl='' #for panels labels will be placed using mtext function
yl=''
}else{ #single plot
loops=length(subject)
N=c(1,1)
M=c(0,0,0,0)
plot=1
WIDTH=800
xl=x_label
yl=y_label
}
for(s in 1:loops){ #this loop is necessary when panels of 8 plots are generated; "loops" is equal to the number of subjects when a single plot has to be generated for each subject
if(s!=loops){ #to ensure that panels include in total as many plots as number of subjects and not the extra empty ones
PLOT=plot #8 plots
} else{
PLOT=length(subject) #in the last panel as many plots as the number of remaining subjects
}
if(occ!=FALSE){ #to always display occasion number in the title of the plot
occ_label=occ
}else{
occ_label='all'
}
if (panel==TRUE){
filename=paste0('panel_', 'IDs ',subject[1], '-',subject[PLOT],' ', filter,' ', 'OCC_',occ_label) #the name has to be changed accordingly depending if panels or single plots are generated
}else{
filename=paste0('ID_', subject[1],' ', filter,' ','OCC_',occ_label)
}
fnme=paste0(save_path,filename,'.png')
png(fnme,width=WIDTH,height=800)
par(mfrow=N)
par(mar=c(4.5,4.5,3,1))
par(oma=M)
par(cex.lab=1.2)
par(cex.main=1.2)
par(cex.axis=1.2)
for(i in 1:PLOT){ #The previous loop is not continued but the new one is started because the changes to TAD column (TAD=0) have to be saved for the whole database to be later plotted
subjectunique=exp_data[exp_data[,ID_col]==subject[i], ]
# add additional plot title label if required
additional_title_label=''
if(sum(title_label_col!='')){
for(L in 1:length(title_label_col)){
additional_title_label=paste0(' ',title_label_col[L],'=',subjectunique[,title_label_col[L]][1]) #to later make a vector of titles of disfferent columns and its values for each subject
if (occ!=FALSE){ #multiple dose: dose might change with each occasion
if(title_label_col[L]=='AMT'){
amt_unique=dose_time[dose_time[,ID_col]==subject[i] & dose_time$OCC==occ,'AMT']
additional_title_label=paste0(' ',title_label_col[L],'=',amt_unique[1])
}
if(title_label_col[L]=='DOSE'){
dose_unique=dose_time[dose_time[,ID_col]==subject[i] & dose_time$OCC==occ,'DOSE']
additional_title_label=paste0(' ',title_label_col[L],'=',dose_unique[1])
}
if(title_label_col[L]=='OCC'){
occ_unique=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,'OCC']
additional_title_label=paste0(' ',title_label_col[L],'=',occ_unique[1])
}
}else{
if(title_label_col[L]=='AMT'){
amt_unique=dose_time[dose_time[,ID_col]==subject[i],'AMT']
additional_title_label=paste0(' ',title_label_col[L],'=',amt_unique[1])
}
if(title_label_col[L]=='DOSE'){
dose_unique=dose_time[dose_time[,ID_col]==subject[i],'DOSE']
additional_title_label=paste0(' ',title_label_col[L],'=',dose_unique[1])
}
if(title_label_col[L]=='OCC'){
additional_title_label=paste0(' ',title_label_col[L],'=','all')
}
}
if (L==1){
additional_title_list=c(paste0('ID=',subject[i]), additional_title_label) #to make a list of column titles to be displayed at the top of each plot
}else{
additional_title_list=c(additional_title_list,additional_title_label)
}
}
}else{
additional_title_list=paste0('ID=',subject[i])
}
title_main=paste(as.character(additional_title_list), collapse=", ") #so that column titles are displayed in one line at the top of each plot
#set colours for observation records: red for BLQ, black for all others
if(LLOQFLAG_col %in% names(subjectunique)){
subjectunique[(subjectunique$EVID==0 & subjectunique[,LLOQFLAG_col]==0),'COLOUR']='blue'
subjectunique[(subjectunique$EVID==0 & subjectunique[,LLOQFLAG_col]!=0),'COLOUR']='green'
} else {
subjectunique[subjectunique$EVID==0,'COLOUR']='blue'
}
#Plot data
if(colours_column==TRUE){
COLOUR_occ=exp_data[,'ALL_COLOURS'][exp_data$OCC==occ & exp_data$EVID==0]
COLOUR_all=exp_data[,'ALL_COLOURS'][exp_data$EVID==0]
}
id_dat=subjectunique[subjectunique$EVID==0,] # get observation records for current patient
t_last=max(exp_data$TIME)
if(keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
t_last=max(subjectunique[subjectunique$EVID==0,ind_var_col]) # get TIME of last observation for current patient
}
# set y-axis limits
min_y=min(exp_data[exp_data$EVID==0,DV_col])
if(NONMEM_table==TRUE){
min_y=min(min_y,exp_data[exp_data$EVID==0,IPRED_col])
}
if(all_same_scale_y==FALSE & NONMEM_table==FALSE){
ylims=c(min(id_dat[,DV_col]),max(id_dat[,DV_col]))
}
if(all_same_scale_y==FALSE & NONMEM_table==TRUE){
ylims=c(min(min(id_dat[,DV_col]),min(id_dat[,IPRED_col])),max(max(id_dat[,DV_col]),max(id_dat[,IPRED_col])))
}
if(all_same_scale_y==TRUE & NONMEM_table==FALSE){
ylims=c(min_y,max(exp_data[exp_data$EVID==0,DV_col]))
}
if(all_same_scale_y==TRUE & NONMEM_table==TRUE){
ylims=c(min_y,max(max(exp_data[,DV_col]),max(exp_data[,IPRED_col])))
}
### AARON 2 START ### here we set x-axis limits.
# set x-axis limits
if(occ==FALSE){
if(all_same_scale_x==FALSE){
xlims=c(min(exp_data[exp_data[,ID_col]==subject[i],ind_var_col]),max(exp_data[exp_data[,ID_col]==subject[i],ind_var_col]))
if(LOG)
{
xlims=range(exp_data[exp_data[,ID_col]==subject[i],ind_var_col][exp_data[exp_data[,ID_col]==subject[i],ind_var_col]>0])
}
}
if(all_same_scale_x==TRUE){
xlims=c(min(exp_data[,ind_var_col]),max(exp_data[,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[,ind_var_col][exp_data[,ind_var_col]>0])
}
}
if(all_same_scale_x==FALSE & keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
# adjust the x-axis maximum if we dont want to include doses after last observation
xlims[2]=max(exp_data[exp_data[,ID_col]==subject[i] & exp_data[,ind_var_col]<=t_last,ind_var_col]) # on reading this again I think xlims[2]=t_last is all that is needed
}
if(all_same_scale_x==TRUE & keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
max_obs_time=max(exp_data[exp_data$EVID==0,ind_var_col])
xlims[2]=max_obs_time
}
}
if(occ!=FALSE & all_same_scale_x==FALSE){
## take care here. This is all to do with plotting data for a particular value of OCC. If a patient has no data for that value of OCC then they will get an empty plot so x-limits still have to be set
fix=0
xl1=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]
xl2=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ & exp_data$EVID==0,ind_var_col]
if(length(xl1)==0){
# no doses or observations for this OCC for current subject
xl1=0.1
fix=1
}
if(length(xl2)==0){
# no observations for this OCC for current subject
xl2=xl1*2
fix=1
}
if(fix==0){
xlims=c(min(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]),max(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ & exp_data$EVID==0,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col][exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]>0])
}
} else {
xlims=c(xl1,xl2)
}
}
if(occ!=FALSE & all_same_scale_x==TRUE){
xlims=c(min(exp_data[,ind_var_col]),max(exp_data[exp_data$EVID==0,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[,ind_var_col][exp_data[,ind_var_col]>0])
}
}
# now make adjustments if axis min = axis max
if(xlims[2]==xlims[1]){
xlims[2]=xlims[1] + 0.1*abs(xlims[1]) + 0.001 # add 0.001 just in case xlims[1]=0
}
if(ylims[2]==ylims[1]){
ylims[2]=ylims[1] + 0.1*abs(ylims[1]) + 0.001 # add 0.001 just in case ylims[1]=0
}
# print(paste0(subject[i],' : ',xlims[1],' ',xlims[2],' : ',ylims[1],' ',ylims[2]))
### AARON 2 END ###
### AARON 3 START ### here we build the plots
if(occ!=FALSE){
if (background_data==TRUE){
X=exp_data[,ind_var_col][exp_data$OCC==occ & exp_data$EVID==0]
x0=plot0(x=X,y=exp_data[,DV_col][exp_data$OCC==occ & exp_data$EVID==0],type='p',col=COLOUR_occ,pch=19, cex=1.0, xlimi=xlims, ylim=ylims, ylab=yl, xlab=xl, main=NULL,...)
title(title_main, line=2)
points0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col=subjectunique[subjectunique$OCC==occ & subjectunique$EVID==0,'COLOUR'],pch=19,cex=1.7)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
if(NONMEM_table==TRUE){
points0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
if(sum(subjectunique$OCC==occ)){ #there might be situation that there are no data records for a certain subject and certain occasion; in this case an empty plot is generated if background data is not requested
X=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0]
x0=plot0(x=X,y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],type='p',col=subjectunique[subjectunique$OCC==occ & subjectunique$EVID==0,'COLOUR'],pch=19,cex=1.7, xlimi=xlims, ylim=ylims, ylab=yl, xlab=xl, main=NULL,...)
title(title_main, line=2)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
if(NONMEM_table==TRUE){
points0(x=X,y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.7)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
plot0(0, type='n', xlab="", ylab="", xaxt='n', yaxt='n') #if there is no data present for a specific subject and occasion then an empty window is shown
title(title_main, line=2)
}
}
}else{
if (background_data==TRUE){
X=exp_data[,ind_var_col][exp_data$EVID==0]
x0=plot0(x=X,y=exp_data[,DV_col][exp_data$EVID==0],type='p',col=COLOUR_all, pch=19, cex=1.0, xlimi=xlims, ylim=ylims, ylab=yl,xlab=xl,main=NULL,...)
title(title_main, line=2)
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='p',col=subjectunique[,'COLOUR'][subjectunique$EVID==0],pch=19,cex=1.7)
if(ind_var_col=='TIME'){
# add lines to DV if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
}
if(NONMEM_table==TRUE){
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
}
if(NONMEM_table==TRUE & ind_var_col=='TIME'){
# add lines to IPREDS if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
X=subjectunique[,ind_var_col][subjectunique$EVID==0]
x0=plot0(x=X,y=subjectunique[,DV_col][subjectunique$EVID==0],type='p',col=subjectunique[,'COLOUR'][subjectunique$EVID==0],pch=19,cex=1.7, xlimi=xlims, ylim=ylims, ylab=yl,xlab=xl,main=NULL,...)
title(title_main, line=2)
if(ind_var_col=='TIME'){
# add lines to DV if TIME is x axis
points0(x=X,y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
}
if(NONMEM_table==TRUE){
points0(x=X,y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
}
if(NONMEM_table==TRUE & ind_var_col=='TIME'){
# add lines to IPREDS if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],x0=x0,y=subjectunique[,IPRED_col][subjectunique$EVID==0],type='l',col='orange',lty=2)
}
}
}
if (dose_marker==TRUE & length(dose_time[dose_time[,ID_col]==subject[i],"ALL_COLOURS"])>0){ #to plot vertical line when the dose was administrated
colours=dose_time[dose_time[,ID_col]==subject[i],"ALL_COLOURS"] #dose colours
id_dose_times=dose_time[dose_time[,ID_col]==subject[i],ind_var_col] #times when the dose was administarted
if(LOG) id_dose_times[id_dose_times<=0]=x0
mark=id_dose_times
for(c in 1:length(mark)){ #to make sure that tickmarks have different colours depending on the dose administrated
axis(side=3, at=mark[c], col.ticks=colours[c], labels=FALSE, lwd.ticks=3, lwd=0, tck=-0.03) #tck - the length of the tickmark (if negative tickamark is plotted outside the plot)
}
if(extend_dose_marker==TRUE){
abline(v=id_dose_times,lty=2,lwd=1,col='gray60')
}
}
if (panel==TRUE){
mtext(y_label, side=2, outer=TRUE, cex=1.5) #to ensure that the plot labels in the panel of 6 plots are only on the outside
mtext(x_label, side=1, outer=TRUE, cex=1.5)
}
}
plot_counter=plot_counter+1
dev.off()
if (panel==TRUE){ #to remove subjects from the vector of all subjects each time a panel of 8 is generated so that new panels are generated for remaining subjects only
subject=subject[!subject%in%subject[1:8]]
}else{
subject=subject[!subject%in%subject[1]] #and the same for single plots
}
if(plot_dose_marker==TRUE & legend_counter==0){
#Plot legend as a separate plot
legend_counter=1
fnme=paste0(save_path,'Legend.png')
png(fnme,width=800,height=400)
par(mfrow=c(2,1))
par(cex.lab=1.2)
par(cex.axis=1.2)
plot(1,type='n',axes=F,ann=F)
if(is.na(LIMIT)==TRUE){
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set according to data in column: ',LLOQ_col)
}else{
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set to: ',LIMIT)
}
print(leg1)
legend('topleft',c('Quantified observation',leg1,'IPRED from NONMEM model'),pt.cex=c(2,2,2),pch=19,col=c('blue','green','orange'),bty='n')
plot(x=dose_unique_sorted,y=rep(1,length(all_dose_colours)),col=all_dose_colours,cex=3,pch=19, ylab='',yaxt='n',xlab="DOSE",main='Legend for colouring of doses',xaxt='n',xlim=range(pretty(dose_unique_sorted,n=8)))
axis(side=1,at=pretty(dose_unique_sorted,n=8),labels=as.character(pretty(dose_unique_sorted,n=8)))
dev.off()
}
if(plot_dose_marker==FALSE & legend_counter==0){
#Plot legend as a separate plot
legend_counter=1
fnme=paste0(save_path,'Legend_no_dose_key.png')
png(fnme,width=800,height=250)
par(mfrow=c(1,1))
par(cex.lab=1.2)
par(cex.axis=1.2)
plot(1,type='n',axes=F,ann=F)
if(is.na(LIMIT)==TRUE){
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set according to data in column: ',LLOQ_col)
}else{
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set to: ',LIMIT)
}
print(leg1)
legend('topleft',c('Quantified observation',leg1,'IPRED from NONMEM model'),pt.cex=c(2,2,2),pch=19,col=c('blue','green','orange'),bty='n')
dev.off()
}
}
print(paste0(plot_counter,' files saved to ',save_path))
return(invisible(all_exp_data))
#-----------------------------------------------------THE END of FUNCTION ----------------------------------------------------
}
######################################################################################################################################
######################################################################################################################################
jgrevel/BAST1-R-Library documentation built on May 21, 2019, 10:11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions NONMEM_plot
Documented in NONMEM_plot
#' Plots data from NONMEM files
#'
#' The tool provides a good overview of the available data for each patient and helps to identify severe outliers or programming errors. This can be useful when
#' doing QC work upon building a new NONMEM data file. The tool was created to visualise any NONMEM data file that contains continuous type DV. The observations
#' can be PK or PD type data, and dosing records don't necessarily have to be present. To identify potential outliers, observations for each individual can be
#' superimposed on observations for all patients.
#'
#' @author Katarzyna Nurzynska
#'
#' @param data data frame containing data in NONMEM format.
#' @param save_path path to the folder where plots will be saved.
#' @param ID_col name of the column in NONMEM file which contains patient identifiers.
#' @param ind_var_col name of the column in NONMEM file which contains the variable for use as the x-axis in plots, usually 'TIME' or 'TAD'.
#' @param DV_col name of the column in NONMEM file which contains observations for plotting on y-axis. If dose corrected DV is to be used then set DV_col='DV_DOSE_CORRECTED'.
#' @param x_unit label for units on x-axis.
#' @param y_unit label for units on y-axis.
#' @param dose_marker logical. If TRUE, add dose marker at the top of the plot. If the AMT or DOSE column is present in the data file this marker will be coloured
#' according to the dose size. A separate plot Legend.png will be generated with a key to colouring of different dose sizes.
#' @param panel logical. If FALSE, generate single plots for each ID, or if TRUE, panels of 8 plots.
#' @param keep_doses_after_last_obs logical. If TRUE, the x-axis of each plot will be scaled such that doses after the last observation are visualised. If FALSE, the
#' x-axis will only reach as far as the last observation.
#' @param calc_dose logical. If TRUE, add DOSE column (showing size of most recent non-zero AMT) to the data. Required if DOSE is not already in the data file
#' and you set either dose_marker=TRUE or DV_col=DV_DOSE_CORRECTED.
#' @param calc_tad logical. If TRUE, add TAD column to the data. Required if TAD is not already in the data file and ind_var_col='TAD'.
#' @param calc_occ logical. If TRUE, add OCC column to the data. Required if OCC is not already in the data file and you want to plot data from different occasions using
#' occ=1 or occ=2 etc.
#' @param occ numeric. Plot data for a particular occasion. Default is FALSE which leads to plots including data from all occasions. If the argument is equal to e.g. 2
#' only data from second occasion is plotted.
#' @param LLOQFLAG_col name of the LLOQFLAG column in NONMEM data file, which indicates if an observation is BLQ (1) or has been quantified (0). BLQ observations are
#' plotted as green circles. Set to FALSE if a LLOQFLAG column is not in the data file.
#' @param LLOQ_col name of the LLOQ column in NONMEM data file, which contains the values of the limit of quantification. If the argument is a number (rather than text
#' such as 'LLOQ') then the lower limit of quantification for every sample will be set to that number, and the BLQ samples will be plotted as green circles at that
#' concentration.
#' @param title_label_col vector of character strings, containing names of the columns in NONMEM data file which will be used for supplying some additional information
#' to plot titles. If title_label_col='' then no additional labels are added. For multiple dose studies, the dose shown in the title will correspond to the selected
#' occasion. If all occasions are plotted, the value of the dose in the title will correspond to the first occasion.
#' @param colour_shading numeric. Sets the relationship between the dose amounts and the colour shades. Important when the range of doses is wide. Try experimenting
#' with different values to see effect on Legend.png and the dose markers. It should be a positive number. If set to 0 then it will instead be calculated as
#' \eqn{sqrt(max(AMT)/min(AMT))/2}.
#' @param colour_range vector of colour names. Sets the colour transitions for dose markers and background data. A minimum of 2 colours is required and plot colours will
#' be interpolated between the members of colour_range, depending on the size of DOSE.
#' @param backround_data logical. If TRUE, background data (data from all other individuals) are added to the plot.
#' @param plot_IPRED logical. If TRUE, then IPRED data are added to the plots. Obviously the data file will have to be a NONMEM $TABLE file such that IPRED data are
#' available.
#' @param all_same_scale_x logical. If TRUE, x-axis range should be equal for all individuals. If FALSE, x-axis range should be different for each individual.
#' @param all_same_scale_y logical. If TRUE, y-axis range should be equal for all individuals. If FALSE, y-axis range should be different for each individual.
#' @param file_text character. Free text that will be added to each filename. Useful when creating lots of different visualisations from the same data.
#' @param ind_var_conversion numeric. Factor to be used for converting the data in the x-axis. For example set to 1/24 to convert TIME from hours to days.
#' @param extend_dose_marker logical. If TRUE, the dose markers will be extended as grey dashed lines from top to bottom of each plot. This can help determine if
#' observations close in time to a dose actually lie before or after the dosing event.
#'
#' @return The tool will plot observations and dosing records vs time, saving one plot for each individual or a panel of eight plots for subject-rich data files.
#' @details User should consult the latest Best Practice guide at
#'
#' \code{\\\\BAST_SERVER\\BAST_Shared_Files\\QA Documents\\BAST Computer Systems Validation\\final_documents\\Procedure_for_visualizing_NONMEM_files\\}
#'
#' for more in-depth documentation on the use of this function.
#' @note Only records with EVID=0 are plotted.
#' @note Function will accept standard arguments for the default R plot function, such as log='y'.
#' @note Columns 'ID', 'EVID', 'DV' and 'TIME' are expected to be present. Columns 'TAD', 'AMT', DOSE' and 'OCC' (occasion) are optional.
#' @note It is assumed that DV and IPRED are on a linear scale. If not, then the user should correct this before using the NONMEM_plot function.
#' @note Before running the function, the user is expected to filter data file to records of interest e.g. placebo or biomarker type, through filtering based
#' on columns such as CMT or DVID.
#' @export
############################################################################################################################################################################################################################################################################################################################################
#FUNCTION: NONMEM_visualising_tool
#PURPOSE: The tool provides a good overview of the available data for each patient and helps to identify severe outliers or programming errors.
# This can be useful when doing QC work upon building a new NONMEM data file.
# The tool was created to visualise any NONMEM data file that contains continuous type DV.
# The observations could be PK or PD type data, and dosing records don't necessarily have to be present.
# To identify potential outliers, observations for each individual can be superimposed on observations for all patients.
# The tool will plot observations and dosing records vs time, saving one plot for each individual or a panel of eight plots for subject-rich data files.
############################################################################################################################################################################################################################################################################################################################################
#
# Only records with EVID=0 are plotted
# Function will accept standard arguments for the default R plot command, such as log='y'
# Columns 'ID', 'EVID', 'DV' and 'TIME' are expected to be present. Columns 'TAD', 'AMT', DOSE' and 'OCC' (occasion) are optional.
# 'TAD', 'DOSE' and 'OCC' can be added to the data file before plotting by using the following arguments: calc_tad=TRUE, calc_dose=TRUE, calc_occ=TRUE.
# It is assumed that DV and IPRED are on a linear scale. If not, then the user should correct this before using the NONMEM_plot function
# Before running the function, the user is expected to filter data file to records of interest e.g. placebo or biomarker type, through filtering based on columns such as CMT or DVID
#
# Minimum function inputs:
# NAME DEFAULT VALUE DESCRIPTION
# data NONE name of the data frame containing the NONMEM-formatted data
# save_path NONE path to the folder where plots will be saved
#
# Optional function inputs:
# NAME DEFAULT VALUE DESCRIPTION
# ID_col 'ID' Name of the column in NONMEM file which contains patient identifiers.
# ind_var_col 'TIME' Name of the column in NONMEM file which contains the variable for use as the x-axis in plots. 'TAD' is the other main alternative for ind_var_col.
# DV_col 'DV' Name of the column in NONMEM file which contains observations for plotting on y-axis. If dose corrected DV is to be used then set DV_col='DV_DOSE_CORRECTED'.
# x_unit 'hours' Label for units on x-axis.
# y_unit 'ug/L' Label for units on y-axis.
# dose_marker TRUE Add dose marker at the top of the plot. If the AMT or DOSE column is present in the data file this marker will be coloured according to the dose size. A separate plot Legend.png will be generated with a key to colourin g of different dose sizes.
# panel TRUE Generate single plots for each ID when the argument is set to FALSE, or panels of 8 plots when the argument is set to TRUE.
# keep_doses_after_last_obs TRUE When set to TRUE the x-axis of each plot will be scaled such that doses after the last observation are visualised. When set to FALSE the x-xis will only reach as far as the last observation
# calc_dose TRUE Add DOSE column (showing size of most recent non-zero AMT) to the data when the argument is TRUE. Required if DOSE is not already in the data file and you set either dose_marker=TRUE or DV_col=DV_DOSE_CORRECTED.
# calc_tad TRUE Add TAD column if the argument is equal to TRUE. Required if TAD is not already in the data file and ind_var_col='TAD'
# calc_occ TRUE Add OCC column if the argument is equal to TRUE. Required if OCC is not already in the data file and you want to plot data from different occasions using occ=1 or occ=2 etc.
# occ FALSE Plot data for a particular occasion. Default is FALSE which leads to plots including data from all occasions. If the argument is equal to e.g. 2 only data from second occasion is plotted.
# LLOQFLAG_col FALSE Name of the column in NONMEM file which indicates if an observation is BLQ (1) or has been quantified (0). BLQ observations are plotted as green circles. Set to FALSE if a LLOQFLAG column is not in the data file.
# LLOQ_col FALSE Name of the column in NONMEM file which contains the values of the limit of quantification. If the argument is a number (rather than text such as 'LLOQ') then the lower limit of quantification for every sample will be set to that number, and the BLQ samples will be plotted as green circles at that concentration.
# title_label_col c('DOSE','OCC') Names of the columns in NONMEM file which will be used for supplying some additional information to plot titles. If title_label_col='' then no additional labels are added. For multiple dose studies the dose shown in the title will correspond to the selected occasion. If all occasions are plotted the value of the dose in the title will correspond to the first occasion.
# colour_shading 0 Sets the relationship between the dose amounts and the colour shades. Important when the range of doses is wide. Try experimenting with different values to see effect on Legend.png and the dose markers. It should be a positive number. If set to 0 then it will instead be calculated as sqrt(max(AMT)/min(AMT))/2.
# colour_range c('mistyrose','pink','red','red4','black') Sets the colour transitions for dose markers and background data. A minimum of 2 colours is required and plot colours will be interpolated between the members of colour_range, depending on the size of DOSE.
# backround_data FALSE If TRUE background data (data from all other individuals) are added to the plot.
# plot_IPRED FALSE If set to TRUE then IPRED data are added to the plots. Obviously the data file will have to be a NONMEM $TABLE file such that IPRED data are available.
# all_same_scale_x FALSE Determines if x-axis range should be same for all individuals (TRUE) or different for each individual (FALSE).
# all_same_scale_y TRUE Determines if y-axis range should be the same for all individuals (TRUE) or different for each individual (FALSE).
# file_text '' Free text that will be added to each filename. Useful when creating lots of different visualisations from the same data.
# ind_var_conversion 1 Factor to be used for converting the data in the x-axis. For example set to 1/24 to convert TIME from hours to days.
# extend_dose_marker FALSE When set to TRUE the dose markers will be extended as grey dashed lines from top to bottom of each plot. This can help determine if observations close in to time to a dose actually lie before or after the dosing event.
##############################################################
##############################################################
NONMEM_plot=function(data,save_path,ID_col='ID',ind_var_col='TIME',DV_col='DV',x_unit='hours',y_unit='ug/L',dose_marker=TRUE,panel=TRUE,
keep_doses_after_last_obs=TRUE,calc_dose=TRUE,calc_occ=TRUE,calc_tad=TRUE,occ=FALSE,LLOQFLAG_col=FALSE,LLOQ_col=FALSE,
title_label_col=c('DOSE','OCC'),colour_shading=0,
colour_range=c('mistyrose','pink','red','red4','black'),background_data=FALSE,plot_IPRED=FALSE,all_same_scale_x=FALSE,
all_same_scale_y=TRUE,file_text='',ind_var_conversion=1,extend_dose_marker=FALSE,...){
##############
LOG=FALSE # check if x axis is logged
if(!is.null(list(...)$log))
{
LOG=grepl('x',list(...)$log)
}
plot0=function(x,y,xlimi=NULL,...) ## A function added by Aaron for doing nice looking things with log scale on x axis
{
LOG=FALSE
x0=0
if(!is.null(list(...)$log))
{
LOG=grepl('x',list(...)$log)
}
snapx=function(x)
{
usr=par()$usr
opar <- par('xlog','ylog')
par(xlog=FALSE)
par(ylog=FALSE)
w=0.5*strwidth('E')
h=0.5*strheight('W')
polygon(x+c(w,w,-w,-w),usr[3]+c(h,-h,-h,h),col='white',lty=0,xpd=NA)
lines(x+c(w,w,NA,-w,-w),usr[3]+c(h,-h,NA,-h,h),xpd=NA)
par(opar)
}
if(LOG) # & any(x<=0) # this part has been excluded from use here so results have consistency
{
if(!is.null(xlimi))
{
r=log(xlimi,10)
}else{
r=range(log(x[x>0],10),na.rm=TRUE)
}
x0=r[1]-0.2*diff(r)
xr=x
xr[xr<=0]=10**x0
plot(xr,y,xaxt='n',xlim=10**c(x0,r[2]),...)
ax=axisTicks(r+diff(r)*0.05*c(-1,1),log=TRUE)
axis(1,at=c(10**x0,ax),labels=c(0,ax))
snapx(0.5*r[1]+0.5*x0)
}else{
if(is.null(xlimi)) xlimi=range(x,na.rm=TRUE)
plot(x,y,xlim=xlimi,...)
}
return(invisible(10**x0))
}
points0=function(x,y,x0,...)
{
if(par()$xlog)
{
x[x<=0]=x0
}
points(x,y,...)
}
lines0=function(x,y,x0,...)
{
if(par()$xlog)
{
x[x<=0]=x0
}
lines(x,y,...)
}
##############
exp_data=data
# set axes labels
if(ind_var_col=='TIME'){
x_label=paste0('Time (',x_unit,')')
}
if(ind_var_col=='TAD'){
x_label=paste0('Time after dose (',x_unit,')')
}
if(ind_var_col!='TIME' & ind_var_col!='TAD'){
x_label=paste0(ind_var_col,' (',x_unit,')')
}
y_label=paste0(DV_col,' (',y_unit,')')
# rename some variables
all_subjects_colours=colour_range
filter=file_text
NONMEM_table=plot_IPRED
if(LLOQFLAG_col!=FALSE & LLOQ_col==FALSE){
stop("LLOQFLAG_col has been set to a column name but LLOQ_col has been set to FALSE. When LLOQFLAG_col has a name then LLOQ_col must also have a name")
}
if(LLOQFLAG_col==FALSE){
LLOQ_col=FALSE
}
# if LLOQ_col is a number then it is not a column name but it is an LLOQ value to be applied to all samples where LLOQFLAG==1
LIMIT=NA
if(is.numeric(LLOQ_col)==TRUE){
LIMIT=LLOQ_col
}
# if user has specified a particular values of occ to plot, then it makes no sence to have set keep_doses_after_last_obs to TRUE, so set it to FALSE
if(occ!=FALSE){
keep_doses_after_last_obs=FALSE
}
if(!DV_col%in%c('DV','DV_DOSE_CORRECTED')){
stop('The only permitted values for DV_col are: DV or DV_DOSE_CORRECTED')
}
plot_counter=0
legend_counter=0
#-----------------------------------------------------------------------------------------------------------------
if(calc_dose==TRUE & !"AMT"%in%names(exp_data)){
stop("You have requested calculation of DOSE column but the is no AMT column in the data file. DOSE cannot be calculated")
}
# check if dose column is required
if((dose_marker==TRUE | "DOSE"%in%title_label_col | DV_col=="DV_DOSE_CORRECTED") & !"DOSE"%in%names(exp_data)){
# dose column is required and is not already in the data
if(calc_dose==FALSE){
stop("The function arguments require that a DOSE column available. DOSE is not in the data, so you much set calc_dose to TRUE to continue")
}
if(calc_dose==TRUE & !"AMT"%in%names(exp_data)){
stop("You have requested calculation of DOSE column but the is no AMT column in the data file. DOSE cannot be calculated")
}
}
#---------------------------------------------------- Option to add DOSE column to the database ------------------------------
if (calc_dose==TRUE){
exp_data=calc_dose(exp_data,ID_col=ID_col)
}
#------------------------------------------------------ DOSE column ----------------------------------
if("DOSE"%in%names(exp_data) & DV_col=="DV"){
colours_column=TRUE
plot_dose_marker=TRUE
}
if("DOSE"%in%names(exp_data) & DV_col=="DV_DOSE_CORRECTED"){
colours_column=TRUE
plot_dose_marker=TRUE
}
if(!"DOSE"%in%names(exp_data) & DV_col=="DV"){ #in this case function cannot add colours to the plot based on the dose size
colours_column=FALSE
COLOUR_occ='gray'
COLOUR_all='gray'
plot_dose_marker=FALSE
}
#---------------------------------------------------- Option to add TAD column --------------------------
if(calc_tad==TRUE & "EVID"%in%names(exp_data) & ID_col%in%names(exp_data) & "TIME"%in%names(exp_data)){
exp_data=calc_tad(exp_data,ID_col=ID_col)
}
if(ind_var_col=="TAD" & !"TAD"%in%names(exp_data)){
stop("Function arguments require that a TAD column is in the data file. TAD is not in the data file. Try setting calc_tad to TRUE")
}
exp_data[,ind_var_col]=exp_data[,ind_var_col]*ind_var_conversion # Convert independent variable. Often this will be conversion of TIME from hours to days
#---------------------------------------------------- Option to add OCC column to the databse -------------------------------
#Each occasion consisits of data records related to administrated dose and all concentrations which were measured afterwards
#Occasion is equal to zero before administration of the first dose
if(occ!=FALSE & calc_occ==FALSE & !"OCC"%in%names(exp_data)){
stop("Function arguments require OCC column to be present in the data file. Try setting calc_occ to TRUE")
}
if(calc_occ==TRUE & !"AMT"%in%names(exp_data)){
stop("Calculation of OCC column has been requested but AMT column is not in the data file. Cannot continue")
}
if(calc_occ==TRUE & "AMT"%in%names(exp_data)){
exp_data=calc_occ(exp_data,ID_col=ID_col,calc_method='by_obs')
}
#---------------------------------------------------- IPRED --------------------------
if (NONMEM_table==TRUE){
IPRED_col="IPRED"
if (DV_col=='DV_DOSE_CORRECTED'){
exp_data$IPRED_DOSE_CORRECTED=exp_data$IPRED/exp_data$DOSE #Calculate dose-corrected DV
IPRED_col="IPRED_DOSE_CORRECTED"
}
}
#------------------------------------------------ Add column with different colours for different doses ---------------------
if(colours_column==TRUE){
dose_unique=unique(exp_data[,"DOSE"])
dose_unique_sorted=sort(dose_unique)
if(colour_shading<0){
colour_shading=0
}
if(colour_shading==0){
colour_shading=sqrt(max(dose_unique)/min(dose_unique))/2
}
#Add column with with different colours for different doses to be later shown on the plot as e.g. background data
colfunc=colorRampPalette(all_subjects_colours, bias=colour_shading) #Data in the backround to be coloured as shades of pink to identify increase in dose
all_dose_colours=colfunc(length(dose_unique_sorted))
exp_data[,"ALL_COLOURS"]='9999'
for(col in 1:length(all_dose_colours)){
exp_data[which(exp_data[,"DOSE"]==dose_unique_sorted[col]),'ALL_COLOURS']=all_dose_colours[col]
}
}
if("OCC"%in%names(exp_data)){
exp_data[,"OCC_ORIG"]=exp_data[,"OCC"] #the original occasion column with OCC=0 remaining unchanged
exp_data[exp_data$OCC==0,"OCC"]=1 #occasion column used for plotting: records before the first dose are plotted as occasion 1
}
#---------------------------------------------------- Calculate dose-corrected DV -------------------------------------------
#----------------------------------------------- Handling data below the limit of quantification ---------------------------
#For data records below the limit of quantification DV_col is equal to limit of quantification; if column with the value of LOQ does not exist then the limit is equal to LIMIT
if(LLOQFLAG_col!=FALSE & DV_col=="DV_DOSE_CORRECTED"){
exp_data$DV_DOSE_CORRECTED=exp_data$DV/exp_data$DOSE #Calculate dose-corrected DV
if(is.na(LIMIT)==FALSE){
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=LIMIT/exp_data[exp_data[,LLOQFLAG_col]==1,"DOSE"] # DV_DOSE_CORRECTED: use the BLQ LIMIT value if it has been entered. So the LIMIT gets dose-corrected
}else{
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=exp_data[exp_data[,LLOQFLAG_col]==1,LLOQ_col]/exp_data[exp_data[,LLOQFLAG_col]==1,"DOSE"] # DV_DOSE_CORRECTED
}
}
if(LLOQFLAG_col!=FALSE & DV_col=="DV"){
if(is.na(LIMIT)==FALSE){
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=LIMIT # DV
}else{
exp_data[exp_data[,LLOQFLAG_col]==1,DV_col]=exp_data[exp_data[,LLOQFLAG_col]==1,LLOQ_col] # DV
}
}
if(LLOQFLAG_col==FALSE & DV_col=="DV_DOSE_CORRECTED"){
exp_data$DV_DOSE_CORRECTED=exp_data$DV/exp_data$DOSE # DV_DOSE_CORRECTED
}
if(DV_col=="DV_DOSE_CORRECTED"){
# correct any calculations of DV_DOSE_CORRECTED that have problems due to zero dose
exp_data[,DV_col][exp_data$DOSE==0]=exp_data$DV[exp_data$DOSE==0]
}
zero_obs=exp_data[,DV_col][exp_data$EVID==0 & exp_data[,DV_col]<=0]
if(length(zero_obs)>0){
print('WARNING...WARNING...WARNING...')
print(paste0(length(zero_obs),' observations have values less than or equal to zero. This will cause problems if you have requested a logarithmic y axis'))
print('If you are not expecting any observations with value of zero, then check that LLOQFLAG_col and LLOQ_col have been correctly defined')
}
# before TAD and TIME get altered, save the exp_data as all_exp_data, which will be returned by the function at the end
all_exp_data=exp_data
#------------------------- Replace TAD and TIME = 0 with a small value (to avoid problems on a log scale)----------------
### AARON 1 START ### here we make the adjustments to TIME and TAD to enable plotting with log x-axis. The adjustments are actually used with linear x-axis too, which is a bit lazy really
#TIME
#minimum=min(exp_data[exp_data[,"TIME"]!=0,"TIME"]) # we are assuming here that TIME shouldn't normally be negative
#if(minimum>=0.1){
# exp_data[exp_data[,"TIME"]==0,"TIME"]=0.1
#}else{
# exp_data[exp_data[,"TIME"]==0,"TIME"]=minimum
#}
#
#TAD
#if(ind_var_col=="TAD"){
# minimum=min(exp_data[exp_data[,"TAD"]>0,"TAD"]) # remember that TAD could be negative (eg for observations before first dose). For datasets which have records with TAD<0 then log x-axis will be a problem unless the records with negative TAD are filtered before plotting
# if(minimum>=0.1){
# exp_data[exp_data[,"TAD"]==0,"TAD"]=0.1
# }else{
# exp_data[exp_data[,"TAD"]==0,"TAD"]=minimum
# }
#}
### AARON 1 END ###
#-------------------------------------------------------
dose_time=exp_data[exp_data$EVID%in%c(1,4),] #this set of data will be used later to plot vertical lines when dose was administrated
### now check for any subjects who have no observations and remove them
subject=unique(exp_data[,ID_col])
IDS_NULL=NULL
for(g in 1:length(subject)){
obs_recs=exp_data[exp_data[,ID_col]==subject[g] & exp_data$EVID==0,DV_col]
if(length(obs_recs)==0){
IDS_NULL=c(subject[g],IDS_NULL)
}
}
if(length(IDS_NULL)>0){
num_null=length(IDS_NULL)
print(paste0('The following ',num_null,' IDs have no observations and will not be plotted:'))
print(paste0(IDS_NULL,collapse=' '))
exp_data=exp_data[!exp_data[,ID_col]%in%IDS_NULL,]
}
#---------------------------------------------------- Plot data -------------------------------------------------------------
#---------------------------------------------------- Single plot or panel of 8 ---------------------------------------------
subject=unique(exp_data[,ID_col])
print('Generating plots...')
#Plots
if (panel==TRUE){ #panel of 8 plots
loops=ceiling(length(subject)/8) #number of panels each consisiting of 8 plots
N=c(2,4) #2 rows and 4 columns
M=c(4,4,0,0)
plot=8
WIDTH=1500 #width of a figure generated
xl='' #for panels labels will be placed using mtext function
yl=''
}else{ #single plot
loops=length(subject)
N=c(1,1)
M=c(0,0,0,0)
plot=1
WIDTH=800
xl=x_label
yl=y_label
}
for(s in 1:loops){ #this loop is necessary when panels of 8 plots are generated; "loops" is equal to the number of subjects when a single plot has to be generated for each subject
if(s!=loops){ #to ensure that panels include in total as many plots as number of subjects and not the extra empty ones
PLOT=plot #8 plots
} else{
PLOT=length(subject) #in the last panel as many plots as the number of remaining subjects
}
if(occ!=FALSE){ #to always display occasion number in the title of the plot
occ_label=occ
}else{
occ_label='all'
}
if (panel==TRUE){
filename=paste0('panel_', 'IDs ',subject[1], '-',subject[PLOT],' ', filter,' ', 'OCC_',occ_label) #the name has to be changed accordingly depending if panels or single plots are generated
}else{
filename=paste0('ID_', subject[1],' ', filter,' ','OCC_',occ_label)
}
fnme=paste0(save_path,filename,'.png')
png(fnme,width=WIDTH,height=800)
par(mfrow=N)
par(mar=c(4.5,4.5,3,1))
par(oma=M)
par(cex.lab=1.2)
par(cex.main=1.2)
par(cex.axis=1.2)
for(i in 1:PLOT){ #The previous loop is not continued but the new one is started because the changes to TAD column (TAD=0) have to be saved for the whole database to be later plotted
subjectunique=exp_data[exp_data[,ID_col]==subject[i], ]
# add additional plot title label if required
additional_title_label=''
if(sum(title_label_col!='')){
for(L in 1:length(title_label_col)){
additional_title_label=paste0(' ',title_label_col[L],'=',subjectunique[,title_label_col[L]][1]) #to later make a vector of titles of disfferent columns and its values for each subject
if (occ!=FALSE){ #multiple dose: dose might change with each occasion
if(title_label_col[L]=='AMT'){
amt_unique=dose_time[dose_time[,ID_col]==subject[i] & dose_time$OCC==occ,'AMT']
additional_title_label=paste0(' ',title_label_col[L],'=',amt_unique[1])
}
if(title_label_col[L]=='DOSE'){
dose_unique=dose_time[dose_time[,ID_col]==subject[i] & dose_time$OCC==occ,'DOSE']
additional_title_label=paste0(' ',title_label_col[L],'=',dose_unique[1])
}
if(title_label_col[L]=='OCC'){
occ_unique=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,'OCC']
additional_title_label=paste0(' ',title_label_col[L],'=',occ_unique[1])
}
}else{
if(title_label_col[L]=='AMT'){
amt_unique=dose_time[dose_time[,ID_col]==subject[i],'AMT']
additional_title_label=paste0(' ',title_label_col[L],'=',amt_unique[1])
}
if(title_label_col[L]=='DOSE'){
dose_unique=dose_time[dose_time[,ID_col]==subject[i],'DOSE']
additional_title_label=paste0(' ',title_label_col[L],'=',dose_unique[1])
}
if(title_label_col[L]=='OCC'){
additional_title_label=paste0(' ',title_label_col[L],'=','all')
}
}
if (L==1){
additional_title_list=c(paste0('ID=',subject[i]), additional_title_label) #to make a list of column titles to be displayed at the top of each plot
}else{
additional_title_list=c(additional_title_list,additional_title_label)
}
}
}else{
additional_title_list=paste0('ID=',subject[i])
}
title_main=paste(as.character(additional_title_list), collapse=", ") #so that column titles are displayed in one line at the top of each plot
#set colours for observation records: red for BLQ, black for all others
if(LLOQFLAG_col %in% names(subjectunique)){
subjectunique[(subjectunique$EVID==0 & subjectunique[,LLOQFLAG_col]==0),'COLOUR']='blue'
subjectunique[(subjectunique$EVID==0 & subjectunique[,LLOQFLAG_col]!=0),'COLOUR']='green'
} else {
subjectunique[subjectunique$EVID==0,'COLOUR']='blue'
}
#Plot data
if(colours_column==TRUE){
COLOUR_occ=exp_data[,'ALL_COLOURS'][exp_data$OCC==occ & exp_data$EVID==0]
COLOUR_all=exp_data[,'ALL_COLOURS'][exp_data$EVID==0]
}
id_dat=subjectunique[subjectunique$EVID==0,] # get observation records for current patient
t_last=max(exp_data$TIME)
if(keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
t_last=max(subjectunique[subjectunique$EVID==0,ind_var_col]) # get TIME of last observation for current patient
}
# set y-axis limits
min_y=min(exp_data[exp_data$EVID==0,DV_col])
if(NONMEM_table==TRUE){
min_y=min(min_y,exp_data[exp_data$EVID==0,IPRED_col])
}
if(all_same_scale_y==FALSE & NONMEM_table==FALSE){
ylims=c(min(id_dat[,DV_col]),max(id_dat[,DV_col]))
}
if(all_same_scale_y==FALSE & NONMEM_table==TRUE){
ylims=c(min(min(id_dat[,DV_col]),min(id_dat[,IPRED_col])),max(max(id_dat[,DV_col]),max(id_dat[,IPRED_col])))
}
if(all_same_scale_y==TRUE & NONMEM_table==FALSE){
ylims=c(min_y,max(exp_data[exp_data$EVID==0,DV_col]))
}
if(all_same_scale_y==TRUE & NONMEM_table==TRUE){
ylims=c(min_y,max(max(exp_data[,DV_col]),max(exp_data[,IPRED_col])))
}
### AARON 2 START ### here we set x-axis limits.
# set x-axis limits
if(occ==FALSE){
if(all_same_scale_x==FALSE){
xlims=c(min(exp_data[exp_data[,ID_col]==subject[i],ind_var_col]),max(exp_data[exp_data[,ID_col]==subject[i],ind_var_col]))
if(LOG)
{
xlims=range(exp_data[exp_data[,ID_col]==subject[i],ind_var_col][exp_data[exp_data[,ID_col]==subject[i],ind_var_col]>0])
}
}
if(all_same_scale_x==TRUE){
xlims=c(min(exp_data[,ind_var_col]),max(exp_data[,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[,ind_var_col][exp_data[,ind_var_col]>0])
}
}
if(all_same_scale_x==FALSE & keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
# adjust the x-axis maximum if we dont want to include doses after last observation
xlims[2]=max(exp_data[exp_data[,ID_col]==subject[i] & exp_data[,ind_var_col]<=t_last,ind_var_col]) # on reading this again I think xlims[2]=t_last is all that is needed
}
if(all_same_scale_x==TRUE & keep_doses_after_last_obs==FALSE & ind_var_col=='TIME'){
max_obs_time=max(exp_data[exp_data$EVID==0,ind_var_col])
xlims[2]=max_obs_time
}
}
if(occ!=FALSE & all_same_scale_x==FALSE){
## take care here. This is all to do with plotting data for a particular value of OCC. If a patient has no data for that value of OCC then they will get an empty plot so x-limits still have to be set
fix=0
xl1=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]
xl2=exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ & exp_data$EVID==0,ind_var_col]
if(length(xl1)==0){
# no doses or observations for this OCC for current subject
xl1=0.1
fix=1
}
if(length(xl2)==0){
# no observations for this OCC for current subject
xl2=xl1*2
fix=1
}
if(fix==0){
xlims=c(min(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]),max(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ & exp_data$EVID==0,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col][exp_data[exp_data[,ID_col]==subject[i] & exp_data$OCC==occ,ind_var_col]>0])
}
} else {
xlims=c(xl1,xl2)
}
}
if(occ!=FALSE & all_same_scale_x==TRUE){
xlims=c(min(exp_data[,ind_var_col]),max(exp_data[exp_data$EVID==0,ind_var_col]))
if(LOG)
{
xlims=range(exp_data[,ind_var_col][exp_data[,ind_var_col]>0])
}
}
# now make adjustments if axis min = axis max
if(xlims[2]==xlims[1]){
xlims[2]=xlims[1] + 0.1*abs(xlims[1]) + 0.001 # add 0.001 just in case xlims[1]=0
}
if(ylims[2]==ylims[1]){
ylims[2]=ylims[1] + 0.1*abs(ylims[1]) + 0.001 # add 0.001 just in case ylims[1]=0
}
# print(paste0(subject[i],' : ',xlims[1],' ',xlims[2],' : ',ylims[1],' ',ylims[2]))
### AARON 2 END ###
### AARON 3 START ### here we build the plots
if(occ!=FALSE){
if (background_data==TRUE){
X=exp_data[,ind_var_col][exp_data$OCC==occ & exp_data$EVID==0]
x0=plot0(x=X,y=exp_data[,DV_col][exp_data$OCC==occ & exp_data$EVID==0],type='p',col=COLOUR_occ,pch=19, cex=1.0, xlimi=xlims, ylim=ylims, ylab=yl, xlab=xl, main=NULL,...)
title(title_main, line=2)
points0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col=subjectunique[subjectunique$OCC==occ & subjectunique$EVID==0,'COLOUR'],pch=19,cex=1.7)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
if(NONMEM_table==TRUE){
points0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
if(sum(subjectunique$OCC==occ)){ #there might be situation that there are no data records for a certain subject and certain occasion; in this case an empty plot is generated if background data is not requested
X=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0]
x0=plot0(x=X,y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],type='p',col=subjectunique[subjectunique$OCC==occ & subjectunique$EVID==0,'COLOUR'],pch=19,cex=1.7, xlimi=xlims, ylim=ylims, ylab=yl, xlab=xl, main=NULL,...)
title(title_main, line=2)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
if(NONMEM_table==TRUE){
points0(x=X,y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.7)
lines0(x=subjectunique[,ind_var_col][subjectunique$OCC==occ & subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$OCC==occ & subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
plot0(0, type='n', xlab="", ylab="", xaxt='n', yaxt='n') #if there is no data present for a specific subject and occasion then an empty window is shown
title(title_main, line=2)
}
}
}else{
if (background_data==TRUE){
X=exp_data[,ind_var_col][exp_data$EVID==0]
x0=plot0(x=X,y=exp_data[,DV_col][exp_data$EVID==0],type='p',col=COLOUR_all, pch=19, cex=1.0, xlimi=xlims, ylim=ylims, ylab=yl,xlab=xl,main=NULL,...)
title(title_main, line=2)
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='p',col=subjectunique[,'COLOUR'][subjectunique$EVID==0],pch=19,cex=1.7)
if(ind_var_col=='TIME'){
# add lines to DV if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
}
if(NONMEM_table==TRUE){
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
}
if(NONMEM_table==TRUE & ind_var_col=='TIME'){
# add lines to IPREDS if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='l',col='orange',lty=2)
}
}else{
X=subjectunique[,ind_var_col][subjectunique$EVID==0]
x0=plot0(x=X,y=subjectunique[,DV_col][subjectunique$EVID==0],type='p',col=subjectunique[,'COLOUR'][subjectunique$EVID==0],pch=19,cex=1.7, xlimi=xlims, ylim=ylims, ylab=yl,xlab=xl,main=NULL,...)
title(title_main, line=2)
if(ind_var_col=='TIME'){
# add lines to DV if TIME is x axis
points0(x=X,y=subjectunique[,DV_col][subjectunique$EVID==0],x0=x0,type='l',col='black',lty=2)
}
if(NONMEM_table==TRUE){
points0(x=X,y=subjectunique[,IPRED_col][subjectunique$EVID==0],x0=x0,type='p',col='orange',pch=19,cex=1.2)
}
if(NONMEM_table==TRUE & ind_var_col=='TIME'){
# add lines to IPREDS if TIME is x axis
points0(x=subjectunique[,ind_var_col][subjectunique$EVID==0],x0=x0,y=subjectunique[,IPRED_col][subjectunique$EVID==0],type='l',col='orange',lty=2)
}
}
}
if (dose_marker==TRUE & length(dose_time[dose_time[,ID_col]==subject[i],"ALL_COLOURS"])>0){ #to plot vertical line when the dose was administrated
colours=dose_time[dose_time[,ID_col]==subject[i],"ALL_COLOURS"] #dose colours
id_dose_times=dose_time[dose_time[,ID_col]==subject[i],ind_var_col] #times when the dose was administarted
if(LOG) id_dose_times[id_dose_times<=0]=x0
mark=id_dose_times
for(c in 1:length(mark)){ #to make sure that tickmarks have different colours depending on the dose administrated
axis(side=3, at=mark[c], col.ticks=colours[c], labels=FALSE, lwd.ticks=3, lwd=0, tck=-0.03) #tck - the length of the tickmark (if negative tickamark is plotted outside the plot)
}
if(extend_dose_marker==TRUE){
abline(v=id_dose_times,lty=2,lwd=1,col='gray60')
}
}
if (panel==TRUE){
mtext(y_label, side=2, outer=TRUE, cex=1.5) #to ensure that the plot labels in the panel of 6 plots are only on the outside
mtext(x_label, side=1, outer=TRUE, cex=1.5)
}
}
plot_counter=plot_counter+1
dev.off()
if (panel==TRUE){ #to remove subjects from the vector of all subjects each time a panel of 8 is generated so that new panels are generated for remaining subjects only
subject=subject[!subject%in%subject[1:8]]
}else{
subject=subject[!subject%in%subject[1]] #and the same for single plots
}
if(plot_dose_marker==TRUE & legend_counter==0){
#Plot legend as a separate plot
legend_counter=1
fnme=paste0(save_path,'Legend.png')
png(fnme,width=800,height=400)
par(mfrow=c(2,1))
par(cex.lab=1.2)
par(cex.axis=1.2)
plot(1,type='n',axes=F,ann=F)
if(is.na(LIMIT)==TRUE){
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set according to data in column: ',LLOQ_col)
}else{
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set to: ',LIMIT)
}
print(leg1)
legend('topleft',c('Quantified observation',leg1,'IPRED from NONMEM model'),pt.cex=c(2,2,2),pch=19,col=c('blue','green','orange'),bty='n')
plot(x=dose_unique_sorted,y=rep(1,length(all_dose_colours)),col=all_dose_colours,cex=3,pch=19, ylab='',yaxt='n',xlab="DOSE",main='Legend for colouring of doses',xaxt='n',xlim=range(pretty(dose_unique_sorted,n=8)))
axis(side=1,at=pretty(dose_unique_sorted,n=8),labels=as.character(pretty(dose_unique_sorted,n=8)))
dev.off()
}
if(plot_dose_marker==FALSE & legend_counter==0){
#Plot legend as a separate plot
legend_counter=1
fnme=paste0(save_path,'Legend_no_dose_key.png')
png(fnme,width=800,height=250)
par(mfrow=c(1,1))
par(cex.lab=1.2)
par(cex.axis=1.2)
plot(1,type='n',axes=F,ann=F)
if(is.na(LIMIT)==TRUE){
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set according to data in column: ',LLOQ_col)
}else{
leg1=paste0('For visualisation of observations below limit of quantification the DV values are set to: ',LIMIT)
}
print(leg1)
legend('topleft',c('Quantified observation',leg1,'IPRED from NONMEM model'),pt.cex=c(2,2,2),pch=19,col=c('blue','green','orange'),bty='n')
dev.off()
}
}
print(paste0(plot_counter,' files saved to ',save_path))
return(invisible(all_exp_data))
#-----------------------------------------------------THE END of FUNCTION ----------------------------------------------------
}
######################################################################################################################################
######################################################################################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.