---
title: "nonmem2R VPC plot functions"
author: "Magnus Åstrand"
date: "February 2020"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{nonmem2R VPC plot functions}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r,echo=FALSE,message=FALSE,warning=FALSE}
library(nonmem2R)
rm(list=ls())
```

## Introduction
Herein is a demo of functions for VPC plots in the nonmem2R R-package. 
The functions creates VPC plots using output files generated by the vpc command in Pearl Speaks NONMEM (PsN). There are 3 functions; `vpcfig`, `vpcfig2`, and `vpcfig3`. The `vpcfig` function use lattice graphics for plotting and have essentially the same arguments and functionality as the newer function `vpcfig2`. The `vpcfig3` function is an experimental test function loading the raw simulated data and not the by psn aggregated summary of the simulations. The primary focus herein is on vpcfig2 using ggplot.

## vpcfig2
Similar to the GOF functions in nonmem2R `vpcfig2` returns an ggplot object that can be further modified by adding ggplot formatting, e.g. changing axis labels by adding `+labs(x="Time after dose")`.

Formatting of plot symbols, lines, and model regions can be controlled by the arguments 
`col.data`, `cex.data`, `pch.data`, `alpha.data`, `col.line`, `lwd.line`, `lty.line`, `col.segm`, and `alpha.segm`. 

`vpcfig2` will by default add a caption in the bottom of the graph to indicate date when figure was generated together with the full path to the script. _However_, this require that the script name has been set using `set.script.name`, and then the script name will be pasted together with the date and path as derived from `getwd`. If script name has not been set no caption will be added. 


```{r,echo=FALSE,fig.width=7,fig.height=5}
# Get path to the example files included in nonmem2R package
file1 <- system.file("extdata", "vpctab004.dat", package = "nonmem2R")
path1 <-gsub("vpctab004.dat","",file1)
path2 <- system.file("extdata", "vpc001", package = "nonmem2R")
```


Below `vpcfig2` is used with default options without any extra formating.
Here the psn generated files (vpctab-file and vpc-results file) are specified as separate arguments. Alternatively and more convenient it's possible to just use the path to the folder containing these two files, see the following examples. 

```{r,echo=TRUE,fig.width=7,fig.height=5}
# Get path to the example files included in nonmem2R package
file1 <- system.file("extdata", "vpctab004.dat", package = "nonmem2R")
file2 <- system.file("extdata", "vpc_results.csv", package = "nonmem2R")
# produce VPC with default setting, here specifying both vpctab and vpcresult
vpcfig2(vpctab=file1,vpcresult=file2)
```

### Adding ggplot formating

In the below example ggplot formatting have been added with modified axis lables, logarithmic y-scale, and using the ggplot black and white theme.
Also, the input source is specified using the folder path containing the vpctab and vpc results file, path1 is set to `r path1`. 

```{r,echo=TRUE,fig.width=7,fig.height=5}
vpcfig2(path1)+
  labs(x="Time after dose (hours)",y="Plasma concentration (nmol/L)")+
  scale_y_log10()+
  theme_bw()
```

### Using the fy argument for logarithmic scale 
An alternative way to get logarithmic scale is to use the `fy` argument to `vpcfig2`. Below example use a semi-logarithmic scale with an offset of 1 so that zero can be included. To get proper y-axis labels `scale_y_continuous` is added using the `fy` function. The example also show how the placement of the legend can be modified using the built in ggplot functionality.
```{r,echo=TRUE,fig.width=7,fig.height=5}
fy<-function(y){log(y+1)}
y.ticks<-c(0,1,10,100,1000)
vpcfig2(path1,fy=fy)+
  labs(x="Time after dose (hours)",y="Plasma concentration (nmol/L)")+
  scale_y_continuous(breaks=fy(y.ticks),labels=y.ticks,minor_breaks=NULL)+
  theme(legend.position="top")
```


### Re-naming, subsetting and or re-order of strata's
The names for strata's are by defult those produced by psn, `TRTPN == 1` and `TRTPN == 3` in the above examples. The `strata.names` argument can be used to set more informative labels by providing a character-vector with new-labels for each of the stratas. 

The `strata.subset` argument can be used to re-order the stratas and or select a subset of stratas that should be included by providing a integer-vector with indexes which of stratas should be included and in what order. For example `strata.subset=1` will show only the `TRTPN == 1` strata wheares `strata.subset=2:1` will show both strata's byt in the oposite order.

NOTE: length of `strata.names` must allways be equal to the total number of stratas in the VPC _before_ any subsetting using `strata.subset` and in the original order.

See below example using both the `strata.names` and the `strata.subset` argument.
  
```{r,echo=TRUE,fig.width=7,fig.height=6}
vpcfig2(path1,strata.subset=2:1,strata.names=c("Cohort 1","Cohort 2"))
```


### VPC for BLQ data
Nonmem2R version 0.2.2 now also can handle VPC for data below quatification (BLQ) when the PSN VPC command was executed with the `-lloq` option. `vpcfig2` will then display the VPC for observed data and proportion of BLQ in separate panels. Below example show a 2 strata VPC with BLQ data in both strata.

The console will display important information on the faceting y-axis breaks/labels to guide how to add/modify formatting. In particular the y-axis breaks/labels may be challenging. In order to separate the range of data from the range of proportion of BLQ, and to get the panel for BLQ data smaller than the panel for data (>LOQ), the proportions y-data are transformed.

The example below use the fy argument for log-scale of data, and the vpcfig2 default y-axis is modified. Note how the breaks and labels for the BLQ panel is just copied from the output in the console. Also note that it's highly recommended __NOT__ to use `scale_y_log10()` when plotting BLQ VPC's. Simarly it is recommeneded to  use the `ylab` and `censoring.labels` for setting y-axis lables.

```{r,echo=TRUE,fig.width=7,fig.height=6}
f1<-vpcfig2(path2,fy=log,strata.names=c("710mg","1000mg"),ylab="Plasma conc (umol/L)",
            censoring.labels=c("a","BLQ(%)"),xlab="Time after dose (hrs)")
##breaks for DV
bry<-c(0.1,0.2,0.4,0.6,1,2,4,6,10,20,40,60)
mbry<-c(1:10,(1:10)/10,(1:10)*10)
##breaks for BLQ, copied from the vpcfif2 console output
brb<-c(5.55058263829153,6.2985767486289,7.04657085896627)
lbb<-c(0,50,100)
## Get minor breaks for BLQ in between the main breaks
mbrb<-(brb[-1]+brb[length(brb)])/2
f1+scale_y_continuous(breaks=c(brb,log(bry)),
                      labels=c(lbb,bry),minor_breaks=c(mbrb,log(mbry)))
```