In [130]:

```
# Settings
options(repr.plot.width=5, repr.plot.height=5)
options(digits=3)
```

In this hands-on session, we will investigate precision of parameter estimates in a one-compartmental absorption model applied to plasma pharmacokinetics data of the orally administered anti-asthmatic drug theophylline.

Download the data file at http://benzekry.perso.math.cnrs.fr/DONNEES/TP_R/calculus/theophylline_data.txt

Fit the data in `Monolix`

using the model file given at http://benzekry.perso.math.cnrs.fr/DONNEES/TP_R/calculus/one_comp_abs.txt. It is a one-compartmental model with first order linear absorption kinetics. Note that there is only one individual. Thus, remove the random effects and assume normal distribution of the parameters. Use a constant error model

Report the values of the parameters in a vector named `params`

In `R`

, define the correlation matrix from the values given in `Monolix`

.

Using the correlation and standard errors, compute the covariance matrix of the estimates. Hint: $Corr=D\cdot Cov\cdot D$ with $D$ the diagonal matrix with standard errors on the diagonal

Plot the 95\% confidence ellipse between parameters $V$ and $k$. And between parameters $k$ and $k_a$. To do so, you will need the `ellipse`

function from the `mixtools`

package. To install this package, do

In [62]:

```
install.packages("mixtools")
```

Then load the package

In [63]:

```
library(mixtools)
```

Use the help or Google to find how to use the `ellipse`

function and plot the said ellipses.

Install the package `numDeriv`

that will be used to compute the jacobian

Load the package

Load the data. Use the code below to skip the first line but still read the header

In [35]:

```
url_data = "http://benzekry.perso.math.cnrs.fr/DONNEES/TP_R/calculus/theophylline_data.txt"
header = read.table(url_data, nrows = 1, stringsAsFactors = FALSE)
df = read.table(url_data, skip=2, header = FALSE)
colnames(df) = header
```

Put the time and concentration data into two vectors `t`

and `conc`

In the following we will need to use the concept of `function`

. We have previously used many built-in functions from `R`

(e.g. `log()`

, `exp()`

, `seq()`

or the `ellipse()`

function above). Let's see how to define our own functions. For instance, the function $f: x\mapsto x^2$ is defined by

In [191]:

```
f = function(x){
x^2
}
```

In [192]:

```
f(2)
```

We can do more complex functions that depend on several inputs, say time and parameters

In [193]:

```
model_function = function(time, param){
exp(-time*param)
}
```

In [194]:

```
model_function(c(0, 1, 2), 2)
```

We can also define a function that depends on a vector of parameters

In [ ]:

```
model_function = function(time, params){
V = params[1]
k = params[2]
1/V*exp(-k*time)
}
```

Finally, the function can depend on an external parameter, or constant, defined in the global environment, such as the dose $D$

In [197]:

```
D = 4.4
model_function = function(time, params){
V = params[1]
k = params[2]
D/V*exp(-k*time)
}
```

Define a function `one_comp_abs_t`

that takes as input a time vector `time`

and a parameter vector `params`

(containing values for $V, k$ and $k_a$) and computes the concentration in the central compartment of a one-compartmental model with oral absorption (see slides of the class for the formula).

Plot the data and the model solution as a continuous curve (i.e. at some discretized values with a large number of points), on the same graph.

To compute the jacobian of the model with respect to the parameters, define a function `one_comp_abs_vect`

of the parameters only (with time fixed to the time vector `t`

).

Using the `jacobian`

function, compute the jacobian matrix of the model evaluated in the parameter estimates $S$. What should be its dimension?

What is the rank of $S$? Use the `rankMatrix()`

function from the `Matrix`

package to compute it. What do you conclude in terms of identifiability?

Look at the slides of the class to get the formula for the theoretical Fisher information matrix from the model jacobian in the case of nonlinear least-squares estimation. Compute the covariance matrix of the estimates from this formula. You will need to compute the inverse of a matrix, which you can do using the `inv()`

function of the `pracma`

package. You will also need to use the value of $\sigma$, which corresponds to $a$ in `Monolix`

.

Compare the covariance matrix that you have computed and the one from `Monolix`

.

Using the `cov2cor()`

function compute the correlation matrix

Using the `eigen()`

function, look at the eigenvalues of your correlation matrix. Compare with `Monolix`

. Compute the `eigenvalues`

from the covariance matrix. Comment

Define a function `cov_matrix_f`

that takes as input a time vector and returns the theoretical covariance matrix computed from the model jacobian.

The goal of optimal design of the sampling times is to minimize the area of the confidence region for the parameter estimates. A surrogate metrix of this area is given by the determinant of the covariance matrix.

Compute the determinant and relative standard errors for the sampling times in the data

Test the effect of having equally spaced sampling times, with the same number of measurements.

Try increasing the number of sampling measurements

What happens if you put all the sampling times after 5 hours? Or all before 5 hours?

Try to play to beat the sampling of the data. With the same number of measurements, try to find a repartition that gives a smaller value of the determinant.

For given number of sampling measurements (say 3, 5 or 10), find the sampling times that minimize the value of the determinant

In [ ]:

```
```