# 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
install.packages("mixtools")
Then load the package
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
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
f = function(x){
x^2
}
f(2)
We can do more complex functions that depend on several inputs, say time and parameters
model_function = function(time, param){
exp(-time*param)
}
model_function(c(0, 1, 2), 2)
We can also define a function that depends on a vector of parameters
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$
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