In [130]:
# Settings
options(repr.plot.width=5, repr.plot.height=5)

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.

Fit model in Monolix

Fit the data in Monolix using the model file given at 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

Plot confidence ellipse

Get covariance matrix from Monolix

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 ellipse confidence intervals

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]:
also installing the dependency ‘segmented’

The downloaded binary packages are in

Then load the package

In [63]:
mixtools package, version 1.1.0, Released 2017-03-10
This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.

Use the help or Google to find how to use the ellipse function and plot the said ellipses.

Fisher information matrix from the model jacobian

Install the package numDeriv that will be used to compute the jacobian

Load the package

Load data

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

In [35]:
url_data = ""
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

Define model functions

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){
In [192]:

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

In [193]:
model_function = function(time, param){
In [194]:
model_function(c(0, 1, 2), 2)
  1. 1
  2. 7.38905609893065
  3. 54.5981500331442

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]

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]

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 data and fit

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.

Compute covariance matrix from jacobian

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

Optimal design

Test different sampling times

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.

Optimal design

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

In [ ]: