Last time, we focused on the percentile VPC, which was already quite an improvement to the scatter VPC but one may wonder when looking at a percentile VPC, how certain are we about the simulated median and the 10-90% distributions?
In reality, there is a confidence interval around these percentile distributions (e.g. how sure are we about the median) and when is a deviation from the observations acceptable? Since every time we run the Monte Carlo simulation, the samples taken will be different. In order to visualize this, a confidence interval VPC can be generated. Luckily, the majority of the steps will be equal to the ones that we have already performed for the percentile VPC.
Our confidence VPC will be developed on the same data as we did for the percentile VPC. For this example, data was generated in which opportunistic sampling around a time point (e.g. 2, 4, 6, 8, 12, 16h) was performed. This results in a scatter of observations around each time point and we will again use binning during this analysis.
Analysis outline
When we use the PsN function vpc, our dataset used for modelling will be used as the basis for the simulation. This will result in every observation of the dataset being simulated 500/1000 times, depending on the chosen number of samples. The no_of_bins and bin_by_count will not be used in this tutorial since we will create our own binning function. In the end, we want to generate the 80% confidence/prediction interval of the simulated data with the 95% confidence interval around these simulated percentiles.
This is the PsN function to create the simulated dataset from Model1.mod:
vpc -samples=500 -no_of_bins=8 -bin_by_count=1 -dir=vpc_Model1 Model1.mod
Getting started
In R, we first want to define the working directory of our ‘final’ model and the prediction and confidence interval we want to use at a later stage:
##################################################################################
############################ CODING BLOCK ########################################
##################################################################################
# Run-Time Environment: R version 3.4.0
# Short title: Confidence VPC using ggplot2
#
#
##################################################################################
##################################################################################
################
# Clean workspace
rm(list=ls(all=TRUE))
## Load libraries
library(reshape)
library(ggplot2)
library(xpose)
library(dplyr)
#########################################
# Set working directory of model file
setwd("...")
Number_of_simulations_performed <- 500 ## Used as 'samples' argument in PsN vpc function. Needs to be a round number.
#########################################
# Set the percentile intervals
PI <- 80 #specify prediction interval in percentage. Common is 80% interval (10-90%)
CI <- 95 #specify confidence interval in percentage. Common is 95%
# Make a vector for the upper and lower percentile
perc_PI <- c(0+(1-PI/100)/2, 1-(1-PI/100)/2)
perc_CI <- c(0+(1-CI/100)/2, 1-(1-CI/100)/2)
Binning in R
Then, we need to define what kind of bins we want to use (how we specify which observations should be grouped together). We can, for example, divide the bins on different time intervals or by an equal number of observations. In this example we take the easy route, by manually defining the boundaries of the bins around each sampling time point.
######################################### # Specify the bin times manually bin_times <- c(0,3,5,7,10,14,18) Number_of_bins <- length(bin_times)-1
This results in bins between 0-3h, 3-5h, 5-7h, etc. and creates an object which includes the number of bins.
Distribution of observations
Since nothing has changed compared to the percentile VPC, the following code will be used to calculate the distribution of the observations. We again use the piping function from dplyr to calculate the statistics per bin:
############################################################
######### Read dataset with original observations
Obs <- read.csv("Original_dataset.csv",na.strings = ".")
Obs<- Obs[Obs$CMT == 2,] # Remove dosing lines from the dataset
#########################################
######### Create bins in original dataset
### Set bins by backward for loop
for(i in Number_of_bins:1){
Obs$BIN[Obs$TIME <= bin_times[i+1]] <-i
}
############################################
## Calculate confidence intervals of the observations
obs_vpc <- Obs %>%
group_by(BIN) %>% ## Set the stratification identifiers (e.g. TIME, DOSE, CMT, BIN)
summarize(C_median = median(DV, na.rm = T), C_lower = quantile(DV, perc_PI[1], na.rm = T), C_upper = quantile(DV, perc_PI[2], na.rm = T))
## Get median time of each bin of observations for plotting purposes
bin_middle <- Obs %>%
group_by(BIN) %>%
summarize(bin_middle = median(TIME, na.rm = T))
# As a vector
bin_middle <- bin_middle$bin_middle
obs_vpc$TIME <- bin_middle
Distribution of Simulated replicates
In order to create the confidence interval around our percentiles, we first need to calculate the percentiles for each iteration separately. We therefore load in the data, and add a new column with the value of the replicate.
#########################################
######### Load in simulated data
## Search for the generated file by PsN
files <- list.files(pattern = "1.npctab.dta", recursive = TRUE, include.dirs = TRUE)
# This reads the VPC npctab.dta file and save it in dataframe_simulations
dataframe_simulations <- read_nm_tables(paste('.\\',files,sep=""))
dataframe_simulations <- dataframe_simulations[dataframe_simulations$MDV == 0,]
## Set the replicate number to the simulated dataset as a new column
dataframe_simulations$replicate <- rep(1:Number_of_simulations_performed,each=nrow(dataframe_simulations)/Number_of_simulations_performed)
# Set vector with unique replicates
Replicate_vector <- unique(dataframe_simulations$replicate)
We can then again specify which observations are in which bin, using the user specified bin limits.
#########################################
######### Create bins in simulated dataset
Set bins by backward for-loop
for(i in Number_of_bins:1){
dataframe_simulations$BIN[dataframe_simulations$TIME <= bin_times[i+1]] <-i
}
Then, we loop trough all the replicates and calculate the distributions of the predictions which will be saved in the sim_PI dataframe. This dataframe now consists of the summary statistics per replicate (e.g. 500 different median values).
sim_PI <- NULL
## Calculate predictions intervals per bin
for(i in Replicate_vector){
# Run this for each replicate
sim_vpc_ci <- dataframe_simulations %>%
filter(replicate %in% i) %>% # Select an individual replicate
group_by(BIN) %>% # Calculate everything per bin
summarize(C_median = median(DV),
C_lower = quantile(DV, perc_PI[1]),
C_upper = quantile(DV, perc_PI[2])) %>% # Calculate prediction intervals
mutate(replicate = i) # Include replicate number
sim_PI <- rbind(sim_PI, sim_vpc_ci)
}
This dataframe can then be used to calculate the 95% confidence interval of each percentile (10%, 50% and 90%) between all the replicates. We also add the time of the boundaries of each bin to be able to create the rectangles in ggplot later on.
###########################
# Calculate confidence intervals around these prediction intervals calculated with each replicate
sim_CI <- sim_PI %>%
group_by(BIN) %>%
summarize(C_median_CI_lwr = quantile(C_median, perc_CI[1]), C_median_CI_upr = quantile(C_median, perc_CI[2]), # Median
C_low_lwr = quantile(C_lower, perc_CI[1]), C_low_upr = quantile(C_lower, perc_CI[2]), # Lower percentages
C_up_lwr = quantile(C_upper, perc_CI[1]), C_up_upr = quantile(C_upper, perc_CI[2]) # High percentages
)
### Set bin boundaries in dataset
sim_CI$x1 <- NA
sim_CI$x2 <- NA
for(i in 1:Number_of_bins){
sim_CI$x1[sim_CI$BIN == i] <-bin_times[i]
sim_CI$x2[sim_CI$BIN == i] <-bin_times[i+1]
}
Graphics
Instead of depicting the simulated intervals with lines, we use the geom_rect function of ggplot to plot rectangles in our figure. This rectangle is the 95% confidence interval around the corresponding value.
For the observations, we add the distributions with geom_line (black), and add all observations as separate points. And we of course want to save it as a high resolution .png image.
#################################################
############### Plotting #########################
### VPC on a linear scale
CI_VPC_lin <- ggplot() +
## Set the rectangles from the simulations
geom_rect(data=sim_CI, mapping=aes(xmin=x1, xmax=x2, ymin=C_median_CI_lwr, ymax=C_median_CI_upr), fill='red', alpha=0.25) + # Median
geom_rect(data=sim_CI, mapping=aes(xmin=x1, xmax=x2, ymin=C_low_lwr, ymax=C_low_upr), fill='blue', alpha=0.25) + # 10% percentile
geom_rect(data=sim_CI, mapping=aes(xmin=x1, xmax=x2, ymin=C_up_lwr, ymax=C_up_upr), fill='blue', alpha=0.25) + # 90% percentile
# Lines of the observations
geom_line(data = obs_vpc, aes(TIME, C_median), col = 'black', linetype = 'dashed', size = 1.25) +
geom_line(data = obs_vpc, aes(TIME, C_lower), col = 'black', linetype = 'dashed', size = 1.25) +
geom_line(data = obs_vpc, aes(TIME, C_upper), col = 'black', linetype = 'dashed', size = 1.25) +
###### Add observations
geom_point(data = Obs,aes(x=TIME,y=DV),color='black')+
####### Set ggplot2 theme and settings
theme_bw()+
theme(axis.title=element_text(size=9.0))+
theme(strip.background = element_blank(),
strip.text.x = element_blank(),legend.position="none") +
# Set axis labels
labs(x="Time (h)",y="Concentration (mg/L)")+
# Add vertical lines to indicate dosing
geom_vline(xintercept = 0, linetype="dashed", size=1) +
# Set axis
scale_y_continuous(expand=c(0.01,0))+
scale_x_continuous(expand=c(0.01,0))+
## Add title and subtitle
ggtitle("Percentile VPC","")
#####################################
## Save as png
png("Figure_Simulation_Confidence_VPC_with_Observations.png",width = 10, height = 5, units='in',res=600)
print(CI_VPC_lin)
dev.off()
This analysis will result in the following confidence interval VPC for our developed model. We can clearly see that the lines of the observations are well within the simulated distributions, indicating adequate model performance.
COMMENT
Any suggestions or typo’s? Leave a comment or contact me at info@pmxsolutions.com!
This is incredibly helpful. Thank you very much for the step-by-step for plotting VPC without reliance on xpose or psn.
Thank you so much for sharing! It is very helpful!