/
intro-poped.Rmd
321 lines (243 loc) · 12.5 KB
/
intro-poped.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
---
title: "Introduction to PopED"
author: "Andrew Hooker"
#date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
rmarkdown::html_vignette:
#html_document:
toc: true
#toc_depth: 3
#number_sections: true
fig_width: 6
#fig_height: 5
#toc_float:
# collapsed: false
# smooth_scroll: true
vignette: >
%\VignetteIndexEntry{Introduction to PopED}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
if(Sys.getenv("LOGNAME")=="ancho179") devtools::load_all("~/Documents/_PROJECTS/PopED/repo/PopED/")
set.seed(1234)
knitr::opts_chunk$set(
collapse = TRUE
, comment = "#>"
#, fig.width=6
, cache = TRUE
)
```
PopED computes optimal experimental designs for both population and individual studies based on nonlinear mixed-effect models. Often this is based on a computation of the Fisher Information Matrix (FIM).
To get started you need to define
* A model.
* An initial design (and design space if you want to optimize)
* The tasks to perform.
There are a number of functions to help you with these tasks. See `?poped` for more information.
Below is an example to introduce the package. Several other examples are available as r-scripts in the "examples" folder in the
PopED installation directory located at:
```{r,eval=FALSE}
system.file("examples", package="PopED")
```
You can view a list of the example files using the commands:
```{r,eval=FALSE}
ex_dir <- system.file("examples", package="PopED")
list.files(ex_dir)
```
You can then open one of the examples (for example, `ex.1.a.PK.1.comp.oral.md.intro.R`, is very similar to the code found in this vignette) using the following code
```{r,eval=FALSE}
file_name <- "ex.1.a.PK.1.comp.oral.md.intro.R"
ex_file <- system.file("examples",file_name,package="PopED")
file.copy(ex_file,tempdir(),overwrite = T)
file.edit(file.path(tempdir(),file_name))
```
In addition, there is another vignette called "Examples" that explores the new features in each example.
# Define a model
Here we define a one-compartment pharmacokinetic model with linear absorption using an analytical solution. In this case the solution is applicable for both multiple and single dosing. Note that this function is also predefined in PopED as `ff.PK.1.comp.oral.md.CL` (see `?ff.PK.1.comp.oral.md.CL` for more information).
```{r,eval=TRUE}
library(PopED)
packageVersion("PopED")
```
```{r struct_model}
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
N = floor(xt/TAU)+1
f=(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)))
return(list( f=f,poped.db=poped.db))
})
}
```
Next we define the parameters of this function, in this case the between-subject variability (BSV) for each parameter is log-normally distributed (parameter `Favail` is assumed not to have BSV). `DOSE` and `TAU` are defined as covariates (in vector `a`) so that we can optimize their values later.
```{r}
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])
}
```
Now we define the residual unexplained variability (RUV) function, in this case the RUV has both an additive and proportional component.
```{r}
feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- ff(model_switch,xt,parameters,poped.db)
f <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = f*(1+epsi[,1])+epsi[,2]
return(list(y=y,poped.db=poped.db))
}
```
# Create a PopED database
We create a poped database to link the model defined above with a set of model parameters, the initial design and design space for optimization.
In this example, the parameter values are defined for the fixed effects (`bpop`), the between-subject variability variances (`d`) and the residual variability variances (`sigma`). We also fix the parameter `Favail` using `notfixed_bpop`, since we have only oral dosing and the parameter is not identifiable. Fixing a parameter means that we assume the parameter will not be estimated (and is know without uncertainty). In addition, we fix the small additive RUV term, as this term is reflecting the higher error expected at low concentration measurements (limit of quantification measurements) and would typically be calculated from analytical assay methods (for example, the standard deviation of the parameter might be 20% of the limit of quantification).
For the initial design, we define two groups (`m=2`) of 20 individuals (`groupsize=20`), with doses of 20 mg or 40 mg every 24 hours (`a`). The initial design has 5 sample times per individual (`xt`).
For the design space, which can be searched during optimization, we define a potential dose range of between 0 and 200 mg (`mina` and `maxa`), and a range of potential sample times between 0 and 10 hours for the first three samples and between 240 and 248 hours for the last two samples (`minxt` and `maxxt`). Finally, we fix the two groups of subjects to have the same sample times (`bUseGrouped_xt=TRUE`).
```{r}
poped.db <- create.poped.database(
# Model
ff_fun=ff,
fg_fun=sfg,
fError_fun=feps,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(prop=0.04,add=5e-6),
notfixed_sigma=c(1,0),
# Design
m=2,
groupsize=20,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24),
xt=c( 1,2,8,240,245),
# Design space
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
bUseGrouped_xt=TRUE)
```
# Design simulation
First it may make sense to check your model and design to make sure you get what you expect when simulating data. Here we plot the model typical values:
```{r simulate_without_BSV}
plot_model_prediction(poped.db, model_num_points = 300)
```
Next, we plot the expected prediction interval (by default a 95% PI) of the data taking into account the BSV and RUV using the option `PI=TRUE`. This option makes predictions based on first-order approximations to the model variance and a normality assumption of that variance. Better (and slower) computations are possible with the `DV=T`, `IPRED=T` and `groupsize_sim = some large number` options.
```{r simulate_with_BSV}
plot_model_prediction(poped.db,
PI=TRUE,
separate.groups=T,
model_num_points = 300,
sample.times = FALSE)
```
We can get these predictions numerically as well:
```{r}
dat <- model_prediction(poped.db,DV=TRUE)
head(dat,n=5);tail(dat,n=5)
```
# Design evaluation
Next, we evaluate the initial design
```{r}
(ds1 <- evaluate_design(poped.db))
```
We see that the fixed-effect and residual variability parameters are relatively well estimated with this design, but the between-subject variability parameters are less well estimated.
## Evaluate alternative design
What about an alternative design with sparse sampling? For example, what if each individual only has 3 time points at 1, 2 and 245 hours:
```{r}
poped.db.new <- create.poped.database(
# Model
ff_fun=ff,
fg_fun=sfg,
fError_fun=feps,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(prop=0.04,add=5e-6),
notfixed_sigma=c(1,0),
# Design
m=2,
groupsize=20,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24),
xt=c( 1,2,245),
# Design space
minxt=c(0,0,240),
maxxt=c(10,10,248),
bUseGrouped_xt=TRUE)
```
```{r}
(ds2 <- evaluate_design(poped.db.new))
```
## Comparison of designs
The precision on CL is similar with the alternative design but the other parameters are less well estimated.
```{r,results='hide'}
(design_eval <- round(data.frame("Design 1"=ds1$rse,"Design 2"=ds2$rse)))
```
```{r design_summary,echo=FALSE}
knitr::kable(design_eval) #%>%
#kableExtra::kable_styling("striped",full_width = FALSE)
```
Comparing the objective function value (OFV), we see that the alternative design (less samples per subject) has a smaller OFV (=worse). We can compare the two OFVs using efficiency, which tells us the proportion extra individuals that are needed in the alternative design to have the same information content as the original design (around 4 times more individuals than are currently in the design).
```{r}
efficiency(ds2$ofv,ds1$ofv,poped.db)
```
# Design optimization
Now we can optimize the sample times of the original design by maximizing the OFV^[Tip: to make the optimization run faster use the option `parallel = TRUE` in the `poped_optim` command.].
```{r optimize,message = FALSE,results='hide'}
output <- poped_optim(poped.db, opt_xt=TRUE)
```
```{r simulate_optimal_design}
summary(output)
plot_model_prediction(output$poped.db)
```
We see that there are four distinct sample times for this design. This means that for this model, with these exact parameter values, that the most information from the study to inform the parameter estimation is with these sample times.
## Examine efficiency of sampling windows
Of course, this means that there are multiple samples at some of these time points. We can explore a more practical design by looking at the loss of efficiency if we spread out sample times in a uniform distribution around these optimal points ($\pm 30$ minutes).
```{r simulate_efficiency_windows,fig.width=6,fig.height=6}
plot_efficiency_of_windows(output$poped.db,xt_windows=0.5)
```
Here we see the efficiency ($(|FIM_{optimized}|/|FIM_{initial}|)^{1/npar}$) drops below 80% in some cases, which is mostly caused by an increase in the parameter uncertainty of the BSV parameter on absorption (om_KA). Smaller windows or different windowing on different samples might be needed. To investigate see `?plot_efficiency_of_windows`.
## Optimize over a discrete design space
In the previous example we optimized over a continuous design space (sample times could be optimized to be any value between a lower and an upper limit). We could also limit the search to only "allowed" values, for example, only samples taken on the hour are allowed.
```{r, message = FALSE,results='hide'}
poped.db.discrete <- create.poped.database(poped.db,discrete_xt = list(c(0:10,240:248)))
output_discrete <- poped_optim(poped.db.discrete, opt_xt=TRUE)
```
```{r simulate_discrete_optimization}
summary(output_discrete)
plot_model_prediction(output_discrete$poped.db, model_num_points = 300)
```
Here we see that the optimization ran somewhat quicker, but gave a less efficient design.
## Optimize 'Other' design variables
One could also optimize over dose, to see if a different dose could help in parameter estimation .
```{r optimize_dose,message = FALSE,results='hide', eval=FALSE}
output_dose_opt <- poped_optim(output$poped.db, opt_xt=TRUE, opt_a=TRUE)
```
In this case the results are predictable ... higher doses give observations with somewhat lower absolute residual variability leading to both groups at the highest allowed dose levels (200 mg in this case).
## Cost function to optimize dose
Optimizing the dose of a study just to have better model parameter estimates may be somewhat implausible. Instead, let's use a cost function to optimize dose based on some sort of target concentration ... perhaps typical population trough concentrations of 0.2 and 0.35 for the two groups of patients at 240 hours.
First we define the criteria we use to optimize the doses, here a least squares minimization.
```{r}
crit_fcn <- function(poped.db,...){
pred_df <- model_prediction(poped.db)
sum((pred_df[pred_df["Time"]==240,"PRED"] - c(0.2,0.35))^2)
}
crit_fcn(output$poped.db)
```
Now we minimize the cost function
```{r cost_optimization, message = FALSE,results='hide',cache=TRUE}
output_cost <- poped_optim(poped.db, opt_a = TRUE, opt_xt = FALSE,
ofv_fun=crit_fcn,
maximize = FALSE)
```
```{r simulate_cost_optmization}
summary(output_cost)
```
We see that the optimal doses are `r round(output_cost$poped.db$design$a[1,"DOSE"],1)` and `r round(output_cost$poped.db$design$a[2,"DOSE"],1)` for the two groups. This leads to population trough concentrations of 0.2 and 0.35 for the two groups of patients at 240 hours:
```{r}
library(ggplot2)
plot_model_prediction(output_cost$poped.db, model_num_points = 300)+
coord_cartesian(xlim=c(230,250))
```