Home PK intro R NONMEM Calculators Data checker

Code snippets

‹ Back

R code snippets


Click on lines to open the code

Database creation and modifications

Quick database evaluation 

library(readxl)
library(tidyverse)
library(wrapr)
library(lubridate)


"""
Check columns and date time formats. Otherwise it will not work. 
Also check import (na = '.', skiprows etc)

Required columns and naming:
- ID
- DV
- TAD
- AMT
- MDV
- EVID
- Date (format='%m/%d/%y')
- TIME (format='%H:%M:%S')

colist = list of covariates to test matching the covariate column names within the dataset

"""


# SETTINGS:

colist = c('GNDR','HT','WT','AGE','ALB','CRP')
date_present = 'y'  # 'y' or 'n' calendar dates
clock_times = 'y' # 'y' or 'n'  clock times as in H:m:S
stats = 'y' # 'y' or 'n'  print stats of columns

# Read data and format
db <- read_excel("test.xlsx",na = ".")

if(clock_times == 'y'){
  db <- db %>% mutate(TIME = format(TIME,'%H:%M:%S'))
}

# Create doc en start writing 
zz <- file("LOG.txt","w")

writeLines('-----------------------------------------------------',con=zz,sep="\n")
writeLines("LOG FILE: CHECKING DATABASE FOR MISSING VALUES",con=zz,sep="\n")
writeLines("SDTSassen v23.06.20",con=zz,sep="\n")
writeLines('-----------------------------------------------------',con=zz,sep="\n")
writeLines("",con=zz,sep="\n")

if(date_present == 'y'){
  CMT.count = 0
  writeLines("CMT missing:",con=zz,sep="\n")
  for(i in 1:nrow(db)){
    if(is.na(db$CMT[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      CMT.count =+ 1
    }
  }
  if(CMT.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("DV missing (MDV=0):",con=zz,sep="\n")
  DV.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$DV[i]) & db$MDV[i]==0 & !is.na(db$MDV[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      DV.count =+ 1
    }
  }
  if(DV.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("AMT missing (EVID=1 or 4):",con=zz,sep="\n")
  AMT.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$AMT[i]) & (db$EVID[i]==1 | db$EVID[i]==4) & !is.na(db$EVID[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      AMT.count =+ 1
    }
  }
  if(AMT.count == 0){writeLines("-None-",con=zz,sep="\n")}

  
  writeLines('',con=zz,sep="\n")
  writeLines("RATE missing (EVID=1 or 4):",con=zz,sep="\n")
  RATE.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$RATE[i]) & (db$EVID[i]==1 | db$EVID[i]==4) & !is.na(db$EVID[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      RATE.count =+ 1
    }
  }
  if(RATE.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  
  EVID.count = 0
  writeLines('',con=zz,sep="\n")
  writeLines("EVID missing:",con=zz,sep="\n")
  for(i in 1:nrow(db)){
    if(is.na(db$EVID[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      EVID.count =+ 1
    }
  }
  if(EVID.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("ID missing:",con=zz,sep="\n")
  ID.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$ID[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      ID.count =+ 1
    }
  }
  if(ID.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("CMT > 3:",con=zz,sep="\n")
  CMT2.count = 0
  for(i in 1:nrow(db)){
    if(!is.na(db$CMT[i]) & db$CMT[i] >3){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i], db$CMT[i])
      writeLines(err,con=zz,sep="\n")
      CMT2.count =+ 1
    }
  }
  if(CMT2.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("DATE missing:",con=zz,sep="\n")
  DATE.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$DATE[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      DATE.count =+ 1
    }
  }
  if(DATE.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("TIME missing:",con=zz,sep="\n")
  TIME.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$TIME[i])){
      err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      TIME.count =+ 1
    }
  }
  if(TIME.count == 0){writeLines("-None-",con=zz,sep="\n")}
}


if(date_present == 'n'){
  CMT.count = 0
  writeLines("CMT missing:",con=zz,sep="\n")
  for(i in 1:nrow(db)){
    if(is.na(db$CMT[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      CMT.count =+ 1
    }
  }
  if(CMT.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("DV missing (MDV=0):",con=zz,sep="\n")
  DV.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$DV[i]) & db$MDV[i]==0 & !is.na(db$MDV[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      DV.count =+ 1
    }
  }
  if(DV.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("AMT missing (EVID=1 or 4):",con=zz,sep="\n")
  AMT.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$AMT[i]) & (db$EVID[i]==1 | db$EVID[i]==4) & !is.na(db$EVID[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      AMT.count =+ 1
    }
  }
  if(AMT.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  
  writeLines('',con=zz,sep="\n")
  writeLines("RATE missing (EVID=1 or 4):",con=zz,sep="\n")
  RATE.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$RATE[i]) & (db$EVID[i]==1 | db$EVID[i]==4) & !is.na(db$EVID[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      RATE.count =+ 1
    }
  }
  if(RATE.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  EVID.count = 0
  writeLines('',con=zz,sep="\n")
  writeLines("EVID missing:",con=zz,sep="\n")
  for(i in 1:nrow(db)){
    if(is.na(db$EVID[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      EVID.count =+ 1
    }
  }
  if(EVID.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("ID missing:",con=zz,sep="\n")
  ID.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$ID[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      ID.count =+ 1
    }
  }
  if(ID.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("CMT > 3:",con=zz,sep="\n")
  CMT2.count = 0
  for(i in 1:nrow(db)){
    if(!is.na(db$CMT[i]) & db$CMT[i] >3){
      err = paste('-',db$ID[i], db$TIME[i], db$CMT[i])
      writeLines(err,con=zz,sep="\n")
      CMT2.count =+ 1
    }
  }
  if(CMT2.count == 0){writeLines("-None-",con=zz,sep="\n")}
  
  writeLines('',con=zz,sep="\n")
  writeLines("TIME missing:",con=zz,sep="\n")
  TIME.count = 0
  for(i in 1:nrow(db)){
    if(is.na(db$TIME[i])){
      err = paste('-',db$ID[i], db$TIME[i])
      writeLines(err,con=zz,sep="\n")
      TIME.count =+ 1
    }
  }
  if(TIME.count == 0){writeLines("-None-",con=zz,sep="\n")}
}


# Check chronological order

if(date_present == 'y'){
  db <- db %>% mutate(DATE = format(DATE,format='%m/%d/%y'))
  db$DATE2 <- mdy(db$DATE)
  dbdt <- db %>% mutate(TIME = format(TIME, format='%H:%M:%S')) %>% filter(!is.na(TIME))
  dbdt$TIME2 <- hms(dbdt$TIME)
  dbdt$DTTM <- dbdt$DATE2 + dbdt$TIME2
  dbdt <- dbdt %>% filter(!is.na(DTTM))
  
  writeLines('',con=zz,sep="\n")
  writeLines("Check chronological order:",con=zz,sep="\n")
  for(i in 2:nrow(dbdt)){
    if(dbdt$ID[i] == dbdt$ID[i-1]){
      if(dbdt$DTTM[i] < dbdt$DTTM[i-1]){
        err = paste('-',db$ID[i-1],db$DATE[i-1], db$TIME[i-1],'gevolgd door', db$ID[i],db$DATE[i], db$TIME[i])
        writeLines(err,con=zz,sep="\n")
      }
    }
  }
}



# Other values

writeLines('',con=zz,sep="\n")
writeLines('',con=zz,sep="\n")
writeLines('',con=zz,sep="\n")
writeLines('-----------------------------------------------------',con=zz,sep="\n")
writeLines('Checking other selected columns for missing values',con=zz,sep="\n")
writeLines('-----------------------------------------------------',con=zz,sep="\n")

if(date_present == 'y'){
  for(a in colist){
    writeLines('',con=zz,sep="\n")
    writeLines(a,con=zz,sep="\n")
    loop.count = 0
    for(i in 1:nrow(db)){
      if(is.na(db[i,a][])){
        err = paste('-',db$ID[i],db$DATE[i], db$TIME[i])
        writeLines(err,con=zz,sep="\n")
        loop.count =+ 1
      }
    }
    if(loop.count == 0){writeLines("-None-",con=zz,sep="\n")}
  }
} else {
  for(a in colist){
    writeLines('',con=zz,sep="\n")
    writeLines(a,con=zz,sep="\n")
    loop.count = 0
    for(i in 1:nrow(db)){
      if(is.na(db[i,a][])){
        err = paste('-',db$ID[i],db$TIME[i])
        writeLines(err,con=zz,sep="\n")
        loop.count =+ 1
      }
    }
    if(loop.count == 0){writeLines("-None-",con=zz,sep="\n")}
  }
}


# Iterate over columns and summarize quick statistics

writeLines('',con=zz,sep="\n")
writeLines('-----------------------------------------------------',con=zz,sep="\n")
writeLines('Brief raw summary of columns',con=zz,sep="\n")
writeLines('-----------------------------------------------------',con=zz,sep="\n")

if(stats == 'y'){
  for(i in 1:ncol(db)){
    if(!is.character(db[[i]])){
      nm <- colnames(db[i])
      summ <- summary(db[,i])
      writeLines(nm,con=zz,sep="\n")
      writeLines(summ,con=zz,sep="\n")
      writeLines('',con=zz,sep="\n")
    }
  }
} else {
    writeLines('',con=zz,sep="\n")
    writeLines('Print stats was set to no',con=zz,sep="\n")
}

close(zz)

# rm(AMT.count, CMT.count, CMT2.count, EVID.count, DV.count, ID.count, 
#    loop.count, RATE.count, summ, i,  DATE.count, a, colist, err, TIME.count, zz, dbdt)
Time after dose 

library(tidyverse)
library(readxl)
library(writexl)

"""
# Import dataset with ADDL and II DV MDV and date time columns
# Check ADDL and II. Also check new individuals
# 'dttm' is the date time column
"""

db <- read_excel("............", na = ".")

TAD <- c()

for(i in 1:nrow(db)){
  if(db$MDV[i] == 1){
    if(!is.na(db$ADDL[i])){
      ldos <- db$dttm[i]
      TAD <- c(TAD,0)
      II <- db$II[i]
    }else{
      ldos <- db$dttm[i]
      TAD <- c(TAD,0)
      II <- 0
    }
  } else if(db$MDV[i]==0 & II==0) {
    TAD <- c(TAD,(as.numeric(difftime(db$dttm[i],ldos, units = c('hours')))))
  } else {
    TAD <- c(TAD,(as.numeric(difftime(db$dttm[i],ldos, units = c('hours')))%%II))
  }
}

newdb <- cbind(db,TAD)
newdb <- newdb %>% select(ID,DOSE,DV,dttm,TAD)


## Write to file ##
write_xlsx(newdb,  "TAD.xlsx")

## Open WD ## 
shell.exec(getwd())
Time after dose (simple version) 

''' 
This code assumes TIME as time column (no clock times) and no ADDL and II. 
Dose is in the AMT column and NA values in the dataset are '.'
'''

# load dataset e.g. a csv file or excel as in the example above
db <- read.csv(file = 'dataset.csv')

# create empty vector to store TAD values for each row
TAD <- c()

# Loop over all rows in database. If it is a dose TAD = 0 and 
# lastdose (ldose) = TIME, which is your time column.
# If there is no dose on that row it takes the TIME of that row minus the ldose (time of last dose)

for(i in 1:nrow(db)){
  if(db$AMT[i] != '.'){
    ldose = db$TIME[i]
    TAD <- c(TAD,0)
  }else{
    TAD <- c(TAD, db$TIME[i] - ldose)
  }
}

# combine time after dose vector with dataset
newdb <- cbind(db,TAD)

# write to csv 
write.csv(newdb, 'ceftriaxon_tad.csv')
Time after start 

'''
Script to calculate cumulative time and time differences per ID
'''

require(tidyverse)
require(lubridate)

## Import dataset with DATE and TIME fields
## This step might require some adjustments for correct representation
## It depends on the format of the DATE and TIME columns

db <- "INPUT DATASET"

db$date <- as_date(db$DATE)
db$time <- hms::as_hms(db$TIME)
db$dttm <- ymd_hms(paste(db$date, db$time, sep = '-'))


## Calculate the time difference between consecutive rows
## It resets for every new ID

timedif <- c(0)
start <- db$dttm[1]

for(i in 2:dim(db)[1]){
  if(db$ID[i] == db$ID[i-1]){
    timedif <- c(timedif, time_length(db$dttm[i] - start, unit='hours'))}
  else {
    timedif <- c(timedif, 0)
    start <- db$dttm[i]}
}

new_db <- cbind(db, timedif)
Check TAD with ADDL 

'''
Script to check whether ADDL is correctly done by calculating the TADs
'''

IF (NEWIND.LT.2.OR.EVID.EQ.3) THEN
                TDOSA=-999 ; TIME OF MOST RECENT DOSE. -999 IF NO PREVIOUS DOSE.
        TADA=0.0   ; TIME AFTER DOSE ADMINISTERED
ENDIF

IF (DOSTIM==0) THEN ; IGNORE NON-EVENT DOSE TIMES (DOSTIM>0)
        IF(EVID.EQ.1.OR.EVID.EQ.4) THEN ; NEW DOSE EVENT RECORD
           TADA=0.0
           DIV=II
           TDOSA=TIME
           TLAST=TDOSA+ADDL*II ; TLAST IS THE TIME OF THE FINAL IMPLIED DOSE
        ENDIF
                IF (TDOSA.GE.0.AND.EVID.NE.1.AND.EVID.NE.4) THEN ; THERE WAS AN EARLIER DOSE
           IF (TIME>TLAST) THEN
               TADA=TIME-TLAST ; CURRENT TIME IS PAST THE TIME OF THE FINAL IMPLIED DOSE
           ELSE
              DIFF=TIME-TDOSA
              TADA=MOD(DIFF,DIV) ; COMPUTES TIME OF THE MOST RECENT IMPLIED DOSE
           ENDIF
        ENDIF ; NOT A DOSE
   ENDIF ; END OF DOSTIM==0
Add ADDL for NONMEM (additional dosages). 
library(lubridate)
library(tidyverse)

'''

## Note: for datasets with fixed II and databases with only sample and dose rows.
## Requires column 'TAS' which is time after start for each patients (cumulative time) and MDV

''' 

# First, specify your dose interval
II <- 24

db1 <- db1 %>% mutate('lastdos' = case_when(lag(MDV==1)~lag(TAS))) %>%
  fill(lastdos) %>% mutate('dosingdiff' = case_when(MDV==1~(TAS-lastdos)/II)) %>% 
  mutate('iADDL' = round(dosingdiff-1)) %>%
  mutate('iADDL' = ifelse(MDV==1 & ID != lag(ID),'.', iADDL),'ADDL' = 0)
  
# We have now created a column with all ADDL values.
# Now we have to get them to the right spot (Notice that every ADDL values belongs to the prior dosing event)

# Create a list with all iADDL values without the 'NA' values
dADDL <- na.omit(db1$iADDL)

# Check whether the line is a dosing event and give the right ADDL value to the column 'ADDL' 
l <- 1
for(k in 1:nrow(db1)){
  if(db1$MDV[k]==1){
    db1$ADDL[k]<-dADDL[l]
    l <-l+1
  }else{
    db1$ADDL[k]<-NA}
}

# Specify the dosing interval in your dataset.
# Create ADDL so dosages will be given up to last observation

db2<- db1 %>% 
  mutate('II' = case_when(!is.na(ADDL)&ADDL!=0~II)) %>% group_by(ID) %>% mutate('lastsamp' = max(TAS, na.rm=T)) %>%
  mutate('newADDL' = case_when(ADDL=='.'~round(((lastsamp-TAS)/II)-1,digits = 0))) %>%
  mutate('ADDL' = ifelse(!is.na(newADDL)&newADDL>0.5,newADDL,ADDL))

# Remove temporary columns
db3 <- db2 %>% select(-c(lastdos,dosingdiff,newADDL,iADDL,lastsamp))
Export to excel 

'''
######## Export to excel #############
'''

require(writexl)

# Get date and add to filename
today <- Sys.Date()
today <- format(today, format="%d%m%y")

# give it a name in this case 'Export_' add date and '.xlsx'
file <- paste("Export_", today, ".xlsx", sep="")

# Write to file. 'db' is the name of the database you want to export to excel.
write_xlsx(db, file)


Diagnostics

Goodness of fit plots (GoF)* 

require(gridExtra)
require(lattice)
require(xpose4)

"""
# Make sure you are in the right working directory containing the tables and check run nummer

# To change 'idv' use: change.xvardef(newdb,var='idv') <- 'TAD'
# To stratify use: dv.vs.pred(newdb, by='GENDER',...)
# To take a subset use: dv.vs.pred(newdb, subset='BLOCK==1',...)
"""
	
new.runno <- '085' 
newdb <- xpose.data(new.runno)

basic.gof(newdb)


#-------- Seperate graphs -------------#

# dv_vs_pred:
dv.vs.pred(newdb, col="black", smooth = T, pch=1, type ="p",
	cex=1.2, xlb=list("Population predictions"), ylb=list("Observations"), main = NULL)

# dv_vs_ipred:
dv.vs.ipred(newdb, col="black", smooth = T, pch=1, type ="p",
	cex=1.2, xlb=list("Individual predictions"), ylb=list("Observations"), main = NULL)

# |IWRES|:
absval.iwres.vs.idv(newdb, col="black", smooth = T, pch=1, type ="p",
cex=1.2, xlb=list("Time after dose (h)"), ylb=list("|IWRES|"), main = NULL)

# CWRES:
cwres.vs.idv(newdb, col="black", smooth = T, pch=1, type ="p",
 cex=1.2, xlb=list("Time after dose (h)"), ylb=list("Conditional weighted residuals"), main = NULL)


#-------- Combine graphs into one figure -------------#

plot1 <- dv.vs.pred(newdb, col="black", smooth = T, pch=1, type ="p",
 	cex=1.2, xlb=list("Population predictions"), ylb=list("Observations"), main = NULL)
plot2 <- dv.vs.ipred(newdb, col="black", smooth = T, pch=1, type ="p", 
	cex=1.2, xlb=list("Individual predictions"), ylb=list("Observations"), main = NULL)
plot3 <- absval.iwres.vs.idv(newdb, col="black", smooth = T, pch=1, type ="p", 
	cex=1.2, xlb=list("Time after dose (h)"), ylb=list("|IWRES|"), main = NULL)
plot4 <- cwres.vs.idv(newdb, col="black", smooth = T, pch=1, type ="p", 
	cex=1.2, xlb=list("Time after dose (h)"), ylb=list("Conditional weighted residuals"), main = NULL)

grid.arrange(plot1, plot2, plot3, plot4, ncol=2, nrow=2)


# For log scale add: ,scales = list(y = list(log = 10), x = list(log=10))
# For limits add: ylim=c(-3,3)
Normalised prediction distribution errors (info NPDE [npde_2.0.tar.gz]) 

'''
# Make sure you are in the right working directory containing the tables containing the observed data and simulated data
'''

require(npde)

setwd("")

obs <- read.table('obs.csv', header = FALSE, sep=",", dec=".")
sim <- read.table('sim.csv', header = FALSE, sep=",", dec=".")
autonpde(obs,sim, iid=1, ix=2, iy=3, imdv=4, namsav="npderun", type.graph="pdf")

# iid = column with id in obs
# ix = column with x in obs (time)
# iy = column with y in obs (dv)
# imdv = column with mdv values in obs
Visual predictive check (VPC)*
Using xpose4 (lattice): 
'''
# First run the VPC in e.g. pirana using PsN: 
# VPC -samples=1000 -bin_array=2.5,5,10,15,20 -bin_by_count=0 -idv=TAD -predcorr
# Setwd sets the directory that contains the results and tab file of the VPC.
'''

setwd("C:/Users/Documents/NONMEM/VPCs/vpc001")
library(xpose4)
xpose.VPC(vpc.info = "vpc_results.csv", vpctab = "vpctab", PI.ci = "area", PI.real= T, type ="p")
# Add labels for multiple plots. For example for four VPC plots

xpose.VPC(vpc.info = "vpc_results.csv", vpctab = "vpctab", PI.ci = "area", PI.real= T, type ="p", 
	main.sub = c('A','B','C','D'), main = 'VPC main title', ylb = 'label y-axis', xlb = 'label x axis')
	
# Tweak layout
xpose.VPC(vpc.info = "vpc_results.csv", vpctab = "vpctab", PI.ci = "area", pch = 16, type ="p", PI.real=T,
	PI.ci.area.smooth =TRUE, PI.ci.down.arcol = "azure4",PI.ci.up.arcol = "azure4",
	PI.ci.med.arcol = "darkcyan", main=NULL, main.sub = c('left plot','right plot'))

Using xpose (ggplot): 
'''
# First run the VPC in e.g. pirana using PsN: 
# VPC -samples=1000 -bin_array=2.5,5,10,15,20 -bin_by_count=0 -idv=TAD -predcorr
# Setwd sets the directory that contains the VPC folder (not inside folder)
'''

setwd("C:/Users/Documents/NONMEM/VPCs/")

library(tidyverse)
library(xpose)

xpdb <- xpose_data(runno = '3000', prefix = NULL)

xpdb %>% vpc_data(opt = vpc_opt(pi = c(0.05, 0.95), bins = "data",n_bins = 6), 
                  psn_folder = 'VPC_3000') %>% vpc(smooth = T)


*More information, functions and options for xpose4 can be found here: xpose4 documentation