R: Plot model predictions

plot_model_prediction {PopED}

R Documentation

Plot model predictions

Description

Function plots model predictions for the typical value in the population, individual predictions and data predictions.

Usage

plot_model_prediction(
  poped.db,
  model_num_points = 100,
  groupsize_sim = 100,
  separate.groups = F,
  sample.times = T,
  sample.times.IPRED = F,
  sample.times.DV = F,
  PRED = T,
  IPRED = F,
  IPRED.lines = F,
  IPRED.lines.pctls = F,
  alpha.IPRED.lines = 0.1,
  alpha.IPRED = 0.3,
  sample.times.size = 4,
  DV = F,
  alpha.DV = 0.3,
  DV.lines = F,
  DV.points = F,
  alpha.DV.lines = 0.3,
  alpha.DV.points = 0.3,
  sample.times.DV.points = F,
  sample.times.DV.lines = F,
  alpha.sample.times.DV.points = 0.3,
  alpha.sample.times.DV.lines = 0.3,
  y_lab = "Model Predictions",
  facet_scales = "fixed",
  facet_label_names = T,
  model.names = NULL,
  DV.mean.sd = FALSE,
  PI = FALSE,
  PI_alpha = 0.3,
  ...
)

Arguments

`poped.db`	A PopED database.
`model_num_points`	How many extra observation rows should be created in the data frame for each group or individual per model. If used then the points are placed evenly between `model_minxt` and `model_maxxt`. This option is used by `plot_model_prediction` to simulate the response of the model on a finer grid then the defined design. If `NULL` then only the input design is used. Can be a single value or a vector the same length as the number of models.
`groupsize_sim`	How many individuals per group should be simulated when DV=TRUE or IPRED=TRUE to create prediction intervals?
`separate.groups`	Should there be separate plots for each group.
`sample.times`	Should sample times be shown on the plots.
`sample.times.IPRED`	Should sample times be shown based on the IPRED y-values.
`sample.times.DV`	Should sample times be shown based on the DV y-values.
`PRED`	Should a PRED line be drawn.
`IPRED`	Should we simulate individual predictions?
`IPRED.lines`	Should IPRED lines be drawn?
`IPRED.lines.pctls`	Should lines be drawn at the chosen percentiles of the IPRED values?
`alpha.IPRED.lines`	What should the transparency for the IPRED.lines be?
`alpha.IPRED`	What should the transparency of the IPRED CI?
`sample.times.size`	What should the size of the sample.times be?
`DV`	should we simulate observations?
`alpha.DV`	What should the transparency of the DV CI?
`DV.lines`	Should DV lines be drawn?
`DV.points`	Should DV points be drawn?
`alpha.DV.lines`	What should the transparency for the DV.lines be?
`alpha.DV.points`	What should the transparency for the DV.points be?
`sample.times.DV.points`	TRUE or FALSE.
`sample.times.DV.lines`	TRUE or FALSE.
`alpha.sample.times.DV.points`	What should the transparency for the sample.times.DV.points be?
`alpha.sample.times.DV.lines`	What should the transparency for the sample.times.DV.lines be?
`y_lab`	The label of the y-axis.
`facet_scales`	Can be "free", "fixed", "free_x" or "free_y"
`facet_label_names`	TRUE or FALSE
`model.names`	A vector of names of the response model/s (the length of the vector should be equal to the number of response models). It is Null by default.
`DV.mean.sd`	Plot the mean and standard deviation of simulated observations.
`PI`	Plot prediction intervals for the expected data given the model. Predictions are based on first-order approximations to the model variance and a normality assumption of that variance. As such these computations are more approximate than using `DV=T` and `groupsize_sim = some large number`.
`PI_alpha`	The transparency of the PI.
`...`	Additional arguments passed to the `model_prediction` function.

Value

A ggplot object. If you would like to further edit this plot don't forget to load the ggplot2 library using library(ggplot2).

Examples

## Warfarin example from software comparison in:
## Nyberg et al., "Methods and software tools for design evaluation 
##   for population pharmacokinetics-pharmacodynamics studies", 
##   Br. J. Clin. Pharm., 2014. 

library(PopED)

## find the parameters that are needed to define from the structural model
ff.PK.1.comp.oral.md.CL

## -- parameter definition function 
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c(CL=bpop[1]*exp(b[1]),
               V=bpop[2]*exp(b[2]),
               KA=bpop[3]*exp(b[3]),
               Favail=bpop[4],
               DOSE=a[1])
    return(parameters) 
}

## -- Define initial design  and design space
poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL",
                                  fg_file="sfg",
                                  fError_file="feps.prop",
                                  bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), 
                                  notfixed_bpop=c(1,1,1,0),
                                  d=c(CL=0.07, V=0.02, KA=0.6), 
                                  sigma=0.01,
                                  groupsize=32,
                                  xt=c( 0.5,1,2,6,24,36,72,120),
                                  minxt=0,
                                  maxxt=120,
                                  a=70)

##  create plot of model without variability 
plot_model_prediction(poped.db)

##  create plot of model with variability by simulating from OMEGA and SIGMA
plot_model_prediction(poped.db,IPRED=TRUE,DV=TRUE)

##  create plot of model with variability by 
##  computing the expected variance (using an FO approximation) 
##  and then computing a prediction interval 
##  based on an assumption of normality
##  computation is faster but less accurate 
##  compared to using DV=TRUE (and groupsize_sim = 500)
plot_model_prediction(poped.db,PI=TRUE)

##-- Model: One comp first order absorption + inhibitory imax
## -- works for both mutiple and single dosing  
ff <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    
    y=xt
    MS <- model_switch
    
    # PK model
    N = floor(xt/TAU)+1
    CONC=(DOSE*Favail/V)*(KA/(KA - CL/V)) * 
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - 
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))  
    
    # PD model
    EFF = E0*(1 - CONC*IMAX/(IC50 + CONC))
    
    y[MS==1] = CONC[MS==1]
    y[MS==2] = EFF[MS==2]
    
    return(list( y= y,poped.db=poped.db))
  })
}

## -- parameter definition function 
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( V=bpop[1]*exp(b[1]),
                KA=bpop[2]*exp(b[2]),
                CL=bpop[3]*exp(b[3]),
                Favail=bpop[4],
                DOSE=a[1],
                TAU = a[2],
                E0=bpop[5]*exp(b[4]),
                IMAX=bpop[6],
                IC50=bpop[7])
  return( parameters ) 
}


## -- Residual Error function
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- ff(model_switch,xt,parameters,poped.db) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  
  MS <- model_switch
  
  pk.dv <- y*(1+epsi[,1])+epsi[,2]
  pd.dv <-  y*(1+epsi[,3])+epsi[,4]
  
  y[MS==1] = pk.dv[MS==1]
  y[MS==2] = pd.dv[MS==2]
  
  return(list( y= y,poped.db =poped.db )) 
}

poped.db <- 
  create.poped.database(
    ff_fun="ff",
    fError_fun="feps",
    fg_fun="sfg",
    groupsize=20,
    m=3,
    bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,
           E0=1120,IMAX=0.807,IC50=0.0993),  
    notfixed_bpop=c(1,1,1,0,1,1,1),
    d=c(V=0.09,KA=0.09,CL=0.25^2,E0=0.09), 
    sigma=c(0.04,5e-6,0.09,100),
    notfixed_sigma=c(0,0,0,0),
    xt=c( 1,2,8,240,240,1,2,8,240,240),
    minxt=c(0,0,0,240,240,0,0,0,240,240),
    maxxt=c(10,10,10,248,248,10,10,10,248,248),
    discrete_xt = list(0:248),
    G_xt=c(1,2,3,4,5,1,2,3,4,5),
    bUseGrouped_xt=1,
    model_switch=c(1,1,1,1,1,2,2,2,2,2),
    a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24),c(DOSE=0, TAU=24)),
    maxa=c(DOSE=200,TAU=40),
    mina=c(DOSE=0,TAU=2),
    ourzero=0)

##  create plot of model and design 
plot_model_prediction(poped.db,facet_scales="free",
                      model.names = c("PK","PD"))

##  create plot of model with variability by  
##  computing the expected variance (using an FO approximation) 
##  and then computing a prediction interval 
##  based on an assumption of normality
##  computation is faster but less accurate 
##  compared to using DV=TRUE (and groupsize_sim = 500)
plot_model_prediction(poped.db,facet_scales="free",
                      model.names = c("PK","PD"),
                      PI=TRUE,
                      separate.groups = TRUE)