```
using Pumas
using DataFramesMeta
using PumasUtilities
```

# Case Study I - Building a Covariate Model

## 1 First of three case studies

This is the first of three case studies based on the material from UMB Center for Translational Medicine found here. The case studies goes through some basic PK and PKPD workflows from reading data to mapping to population, fitting, diagnostics generation, and much more.

This first case study covers basic functions needed to perform a compartmental PK project. The data is from an IV injection study to keep the model simple and focus on workflows and concepts. The next two case studies feature more advanced concepts.

## 2 Getting the data

Please download the data file for this tutorial and place it in the same location as this tutorial file. The study data consists of pharmacokinetic profiles from 100 individuals, randomized to receive either a 100mg or 250mg IV dose. In the Pumas framework, this data can be represented as a file, a DataFrame, or a JuliaHub DataSet. Here, we load a `DataFrame`

by reading the file using `CSV.read`

. We use the `first(input, number_of_elements)`

function to show the first 10 rows of the `DataFrame`

.

```
using CSV
= CSV.read("CS1_IV1EST_PAR.csv", DataFrame; header = 4)
pkdata first(pkdata, 10)
```

Row | CID | TIME | CONC | AMT | DOSE | MDV | AGE | WT | SCR | ISM | CLCR |
---|---|---|---|---|---|---|---|---|---|---|---|

Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | Float64 | Float64 | Float64 | Int64 | Float64 | |

1 | 1 | 0.0 | 0.0 | 100 | 100 | 1 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

2 | 1 | 0.25 | 13.026 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

3 | 1 | 0.5 | 14.984 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

4 | 1 | 0.75 | 14.16 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

5 | 1 | 1.0 | 19.316 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

6 | 1 | 1.5 | 13.146 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

7 | 1 | 2.0 | 12.921 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

8 | 1 | 2.5 | 8.485 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

9 | 1 | 3.0 | 16.437 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

10 | 1 | 4.0 | 10.724 | 0 | 100 | 0 | 34.823 | 38.212 | 1.1129 | 0 | 42.635 |

The file contains 11 columns:

- CID: subject identification number.
- TIME: Blood sampling time. Needs to be monotonically increasing.
- CONC: The dependent variable. In this case, a plasma concentration.
- AMT: Amount of drug administered at dosing time in the TIME column or zero for observation rows
- DOSE: Dose group. Included for post-processing and asset generation (tables and plots).
- MDV: Missing dependent variable. Can be used to flag observations that should be ignored.
- AGE: Age in years
- WT: Weight in Kg
- SCR: Serum Creatinine (sCr) in mg / dL
- ISM: IS Male (1) or female (0)
- CLCR: Creatinine clearance in ml / min

The dataset variable names reflect its origin from a NONMEM case study. Several conventions used are historical and relate specifically to NONMEM’s functionality. We don’t need to follow these conventions, so we will rename some columns to make them easier to read. First, let’s rename the columns as follows:

`rename!(pkdata, :CID => :ID, :WT => :WEIGHT, :ISM => :SEX)`

Notice, that we rename the `pkdata`

DataFrame column names in pairs of `Symbol`

s. The `Symbol`

s are the column names with `:`

in front. If we omit the `:`

we will get an error. This is because Pumas will then think a variable name is being referred to. We renamed `ISM`

to `SEX`

and we would also like to use self-documenting values instead of 0 and 1. In the following, we overwrite 1 with `"male"`

and 0 with `"female"`

:

`@rtransform! pkdata :SEX = :SEX == 1 ? "male" : "female"`

The final step before our analysis data is ready is to add event columns. These indicate whether a row is an event or an observation. If the row specifies an event it needs to indicate which compartment the dose applies to. The `EVID`

column indicates the type. If it is 0 the row is an observation and if it is 1 the row in a dosing event.

```
@rtransform! pkdata :EVID = :AMT == 0 ? 0 : 1
@rtransform! pkdata :CMT = :AMT == 0 ? missing : Symbol("Central")
```

Instead of using the MDV functionality that would typically be used in NONMEM we set the `CONC`

column to `missing`

for event records.

```
allowmissing!(pkdata, :CONC)
@rtransform! pkdata :CONC = :EVID == 1 ? missing : :CONC
```

## 3 Mapping to a population

In Pumas, we take tabular datasets and convert them into `Population`

s using the `read_pumas`

function. In our case, the population is read using the following code:

```
= read_pumas(
population
pkdata;= :ID,
id = :TIME,
time = [:CONC],
observations = :EVID,
evid = :AMT,
amt = :CMT,
cmt = [:WEIGHT, :CLCR, :AGE, :SEX],
covariates )
```

```
Population
Subjects: 100
Covariates: WEIGHT, CLCR, AGE, SEX
Observations: CONC
```

Notice, that some inputs are listed as `Symbol`

s (a symbol with name “NAME” is written with a colon in front`:NAME`

) and some are vectors of `Symbol`

s. This is because the inputs `observations`

and `covariates`

can potentially have one or multiple names in the `observations`

or zero, one or many in the `covariates`

case. For consistency, we just always ask for a `Vector`

of `Symbol`

s even if one input is given. Once created we see some summary information about how many subjects there are and what kind of information they have.

## 4 Model code for Base Model

To start the pharmacometric analysis we will have to build a model. We will use a simple, one-compartment, IV bolus model with first-order elimination from the central compartment. Such a model has a closed-form solution which we implement using the clinically relevant parameters of clearance (`CL`

) and volume of distribution (`Vc`

). For random effects, we’ll include between-subject variability (BSV) on both CL and Vc, and use a combined (additive + proportional) residual error model.

To define a `PumasModel`

we need to specify the following model parts:

`@param`

is used to define the model’s fixed effects and their corresponding domains. Pumas accepts unicode characters (e.g.,`θcl`

,`σ_add`

) which helps bridge pharmacometric theory and application. Domain functions allow you to specify, lower and upper bounds along with and initial (`init`

) estimate.`@random`

is used to define random effects and their corresponding distributions.`@pre`

is used to pre-process variables for the dynamic system and statistical specification.`@dynamics`

is used to define the dynamics system for the model. The differential equuation for this model would be`Central' = -CL/Vc*Central`

, but this system has a closed-form solution which Pumas implements using`Central1`

. In general, you should use closed-form solutions like`Central1`

when available since they come with a performance boost.`@derived`

defines the error model for the dependent variables.`~`

is used when specifying a distribution.`:=`

is used for intermediate calculations, the results of which are not included in the final output.

In the `derived`

block `@.`

is the “dot call” macro, a full explanation of which, is beyond the scope of this document. In short, the calculations within the `@derived`

block are performed on vectors of numbers and the `@.`

ensures that the operators within a function or expression are “vectorized” appropriately.

In this case, the dynamic model is called `Central1`

and has to be specified in the `@dynamics`

part of the model code. Individual parameters are defined with `RealDomain`

s, the variance matrix `Ω`

is specified as `PDiagDomain`

(so no off-diagonal elements). There are two random effects and the assumed error model is a combined additive and proportional normal error model.

```
= @model begin
iv1cmt @metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
@param begin
∈ RealDomain(lower = 0.1, init = 1.0)
θcl ∈ RealDomain(lower = 1.0, upper = 20.0, init = 10.0)
θvc ∈ PDiagDomain(2)
Ω ∈ RealDomain(lower = 0.0, init = 1.0^2)
σ_add ∈ RealDomain(lower = 0.0, init = 0.09^2)
σ_prop end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= θcl * exp(η[1])
CL = θvc * exp(η[2])
Vc end
@dynamics Central1
@derived begin
:= @. Central / Vc
conc_model ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model * σ_prop)^2))
CONC end
end
```

```
PumasModel
Parameters: θcl, θvc, Ω, σ_add, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: CONC
Observed: CONC
```

Besides the data and model we have to set initial parameters. This is done by specifying a `NamedTuple`

using the syntax below.

```
= (
initial_est_iv1cmt = 0.9,
θcl = 10.0,
θvc = Diagonal([0.1, 0.1]),
Ω = sqrt(0.09),
σ_add = sqrt(0.01),
σ_prop )
```

If all parameters were given `init`

values in `@param`

we could also have used those initial values using the `init_params(iv1cmt)`

function call.

## 5 Fit base model

To fit the model to the data we use the `fit`

function.

`= fit(iv1cmt, population, initial_est_iv1cmt, FOCE()) iv1cmt_results `

```
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.220817e+03 7.627370e+03
* time: 0.018917083740234375
1 4.945610e+03 5.427910e+02
* time: 0.37656712532043457
2 4.870609e+03 4.624623e+02
* time: 0.4175240993499756
3 4.746284e+03 2.387450e+02
* time: 0.4337170124053955
4 4.722871e+03 2.277545e+02
* time: 0.46802401542663574
5 4.697752e+03 2.695471e+02
* time: 0.48253417015075684
6 4.670279e+03 2.251210e+02
* time: 0.5049490928649902
7 4.653649e+03 7.845640e+01
* time: 0.5179400444030762
8 4.651567e+03 3.842645e+01
* time: 0.5302901268005371
9 4.650461e+03 3.919094e+01
* time: 0.5508852005004883
10 4.649212e+03 3.992715e+01
* time: 0.5637540817260742
11 4.647847e+03 4.563022e+01
* time: 0.5844972133636475
12 4.644574e+03 4.320284e+01
* time: 0.5975570678710938
13 4.641136e+03 3.574904e+01
* time: 0.6191802024841309
14 4.639193e+03 3.808485e+01
* time: 0.6309762001037598
15 4.638535e+03 3.905231e+01
* time: 0.6439120769500732
16 4.637660e+03 4.078589e+01
* time: 0.663905143737793
17 4.635157e+03 4.722472e+01
* time: 0.6762912273406982
18 4.620341e+03 1.466178e+02
* time: 0.6967120170593262
19 4.597789e+03 3.472506e+02
* time: 0.7084691524505615
20 4.567135e+03 3.126339e+02
* time: 0.7282130718231201
21 4.491345e+03 1.786065e+02
* time: 0.7449431419372559
22 4.468647e+03 1.802908e+02
* time: 0.7556710243225098
23 4.440016e+03 8.571120e+01
* time: 0.7738850116729736
24 4.429365e+03 9.173855e+01
* time: 0.7840280532836914
25 4.413616e+03 5.147760e+01
* time: 0.8036291599273682
26 4.409551e+03 3.372222e+01
* time: 0.8134231567382812
27 4.406285e+03 1.628709e+01
* time: 0.8245680332183838
28 4.406165e+03 1.665758e+01
* time: 0.8427350521087646
29 4.406084e+03 1.661883e+01
* time: 0.8513450622558594
30 4.405484e+03 1.562746e+01
* time: 0.8622941970825195
31 4.404601e+03 1.315858e+01
* time: 0.8799600601196289
32 4.403168e+03 1.780770e+01
* time: 0.8894641399383545
33 4.402151e+03 1.436821e+01
* time: 0.9008150100708008
34 4.401710e+03 5.107809e+00
* time: 0.9178640842437744
35 4.401646e+03 8.160373e-01
* time: 0.9269471168518066
36 4.401643e+03 4.611581e-01
* time: 0.9372940063476562
37 4.401643e+03 4.575021e-01
* time: 0.9538631439208984
38 4.401643e+03 4.526048e-01
* time: 0.9620401859283447
39 4.401643e+03 4.350247e-01
* time: 0.9719321727752686
40 4.401642e+03 4.353443e-01
* time: 0.9887781143188477
41 4.401642e+03 7.252686e-01
* time: 0.9970190525054932
42 4.401641e+03 8.811155e-01
* time: 1.0069630146026611
43 4.401640e+03 6.601170e-01
* time: 1.0243091583251953
44 4.401639e+03 2.317088e-01
* time: 1.032717227935791
45 4.401639e+03 2.594974e-02
* time: 1.0421290397644043
46 4.401639e+03 1.773922e-03
* time: 1.0593252182006836
47 4.401639e+03 1.773922e-03
* time: 1.0777430534362793
48 4.401639e+03 1.773922e-03
* time: 1.108795166015625
```

```
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -4401.639
Number of subjects: 100
Number of parameters: Fixed Optimized
0 6
Observation records: Active Missing
CONC: 1500 0
Total: 1500 0
--------------------
Estimate
--------------------
θcl 0.41707
θvc 7.1609
Ω₁,₁ 0.1714
Ω₂,₂ 0.19817
σ_add 2.9282
σ_prop 0.12113
--------------------
```

### 5.1 Understanding the output

From the final output we see several lines of information. We see the information about the likelihood approximation, dynamical model solution method and so on. An important line to know about is the `Log-likelihood`

line.

` Log-likelihood value: -4401.639`

During fitting we maximize this value, and it can also be queried programmatically with the `loglikelihood`

function.

`loglikelihood(iv1cmt_results)`

`-4401.638994836392`

If you are comparing NONMEM and Pumas results it is useful to know what to compare the output of `loglikelihood`

to. In the `.lst`

or output file of a NONMEM run you can find a section like the one below.

```
TOTAL DATA POINTS NORMALLY DISTRIBUTED (N): 1500
N*LOG(2PI) CONSTANT TO OBJECTIVE FUNCTION: 2756.264670795735
OBJECTIVE FUNCTION VALUE WITHOUT CONSTANT: 6047.013314836777
OBJECTIVE FUNCTION VALUE WITH CONSTANT: 8803.277985632512
```

We can compare this to the Pumas results by running

`-2loglikelihood(iv1cmt_results)`

`8803.277989672784`

And compare it to the `OBJECTIVE FUNCTION VALUE WITH CONSTANT`

line to verify that the fits match between the two pieces of software.

Besides high level information about the model fit we also see individual parameter estimates. For example, we have the first two lines of parameters that show the population level or fixed effects of the clearance and volume of distribution estimates.

```
θcl 0.41707
θvc 7.1609
```

This implies a clearance of approximately 0.417 L/hr and a volume of distribution of 7.161 L.

There are two ways to view the results of the model fit.

`coef(iv1cmt_results)`

which gives a named tuple of the coefficient values`coeftable(iv1cmt_results)`

which gives a dataframe of the parameters and their estimates.

We can also look at the estimated variances (often referred to as omegas or Ωs) of the random effects (often referred to as etas or ηs).

```
Ω₁,₁ 0.1714
Ω₂,₂ 0.19817
```

The subscripts refer to the position in the variance-covariance matrix Ω. When the row (first number) and column (second number) is identical we are on the diagonal of the matrix. This means we are looking at variances. To calculate the coefficient of variation (%) often used to describe the variability we can use the formula for the Log-Normal distribution

`= sqrt(exp(0.1714) - 1) * 100, CV_Vc = sqrt(exp(0.19817) - 1) * 100) (CV_CL `

```
(CV_CL = 43.23950048893227,
CV_Vc = 46.815556713938285,)
```

For small variances, an approximation is often used

`= sqrt(0.1714) * 100, CV_Vc = sqrt(0.19817) * 100) (CV_CL `

```
(CV_CL = 41.400483088968905,
CV_Vc = 44.51628915352222,)
```

It is possible to specify random effects either using a vector-matrix notation

`η ~ MvNormal(Ω)`

where Ω is specified in the `@param`

block as a `PSDDomain`

(covariances are all present) or `PDiagDomain`

(covariances are zero). It is also possible to specify individual random effects as scalar (single number) distributions

`ηCL ~ Normal(0, ωCL)`

Notice, that the scalar version is specified using lower-case “omega” or ω using unicode characters because scalar Normal distributions are parameterized by the mean and *standard deviation* not mean and variance. You also need to specify the mean as 0. Omitting the 0 incorrectly specifies a variable with a *unit variance* normal distribution with *mean* ωCL. We strongly suggest using the uppercase / lowercase distinction between variance-covariance matrices and scalar standard deviations. If you prefer to use the scalar version above but still estimate a variance in the `@param`

block such that all output is in variances, you can use the “squared” unicode character ² and write

`ωCL² ∈ RealDomain(lower=0.0)`

and

`ηCL ~ Normal(0, sqrt(ωCL²))`

The last two parameters refer to the two components of the error model.

```
σ_add 2.9282
σ_prop 0.12113
```

Again, these are standard deviations which is why we had to square them in the error model specification.

`CONC ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model*σ_prop)^2))`

This concludes looking at the results from the `fit`

function call.

## 6 Create additional output and interpret it

To get access to diagnostic such as predictions and residuals we need to use the `inspect`

function.

`= inspect(iv1cmt_results) iv1cmt_inspect `

```
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.
```

```
FittedPumasModelInspection
Fitting was successful: true
Likehood approximation used for weighted residuals : Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_every = 1, time_limit = NaN, )
)
```

### 6.1 Diagnostic plots for structural model and interpretation

We can look at some basic plots that will help us determine the quality of the fit. We can look at the observations with individual predictions and population predictions overlaid. We will generate them with the `subject_fits`

function and turn on `separate = true`

to show individual subjects on individual plots. To collect them in grids of 9 (default value, can be changed with the `limit`

keyword) plots we use the `paginate = true`

option. To look at the first page of 9 plots we index into the generated variable that holds the plots.

```
= subject_fits(iv1cmt_inspect; paginate = true, separate = true)
all_subject_fits 1] all_subject_fits[
```

Another good set of plots to look at would be `goodness_of_fit`

.

`goodness_of_fit(iv1cmt_inspect)`

Looking at the first plot specifically we see that something must be missing in our model. We are not able to predict the data well if we ignore individual variability. This means that our model will have a hard time extrapolating into unseen data, because the top right plot only looks as good as it does because we could use individual data to calculate individual empirical bayes estimates for the random effects.

## 7 Exploring covariate relationships

To improve our model, we will now look for covariate relationships with empirical bayes estimates (EBEs). In essence, the EBEs are able to adjust to any missing information in the specification of `CL`

and `Vc`

up to some degree. In order to obtain the ability to extrapolate better outside the data, we need to see if we have any measured information about the subjects that could stand in place of the unexplained and unmeasured random effect.

`empirical_bayes_vs_covariates(iv1cmt_inspect)`

It is quite clear that our EBEs seem to pick up some information that is also present in covariates. At this point the “art of modeling” enters the stage, and you should use whatever knowledge you may have of the physiology, pharmacology, or any standard modeling strategies that may apply to your situation.

### 7.1 Covariate model

To improve the model we will add a weight effect on volume of distribution and creatinine clearance effect on clearance.

```
= @model begin
iv1cmt_covar @param begin
∈ RealDomain(lower = 0.1)
θcl ∈ RealDomain(lower = 0.0)
θCLCR ∈ RealDomain(lower = 1.0)
θvc ∈ PDiagDomain(2)
Ω ∈ RealDomain(lower = 0.0)
σ_add ∈ RealDomain(lower = 0.0)
σ_prop end
@covariates CLCR WEIGHT
@random begin
~ MvNormal(Ω)
η end
@pre begin
= θcl * (CLCR / 120)^θCLCR * exp(η[1])
CL = θvc * (WEIGHT / 70)^0.75 * exp(η[2])
Vc end
@dynamics Central1
@derived begin
:= @. Central / Vc
conc_model ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model * σ_prop)^2))
CONC end
end
```

```
PumasModel
Parameters: θcl, θCLCR, θvc, Ω, σ_add, σ_prop
Random effects: η
Covariates: CLCR, WEIGHT
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: CONC
Observed: CONC
```

and we set initial parameters

```
= (
initial_est_iv1cmt_covar = 0.9,
θcl = 0.1,
θCLCR = 10.0,
θvc = Diagonal([0.1, 0.1]),
Ω = sqrt(9.0),
σ_add = sqrt(0.01),
σ_prop )
```

```
(θcl = 0.9,
θCLCR = 0.1,
θvc = 10.0,
Ω = [0.1 0.0; 0.0 0.1],
σ_add = 3.0,
σ_prop = 0.1,)
```

### 7.2 Covariate model results and interpretation

Now, let us fit the covariate model. This time, we set the `optim_options`

keyword to `(show_trace = false,)`

to suppress the progress information while `fit`

is running.

```
= fit(
iv1cmt_covar_results
iv1cmt_covar,
population,
initial_est_iv1cmt_covar,FOCE();
Pumas.= (show_trace = false,),
optim_options )
```

```
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
```

```
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -4325.4154
Number of subjects: 100
Number of parameters: Fixed Optimized
0 7
Observation records: Active Missing
CONC: 1500 0
Total: 1500 0
---------------------
Estimate
---------------------
θcl 0.89494
θCLCR 0.88952
θvc 9.2742
Ω₁,₁ 0.056211
Ω₂,₂ 0.092421
σ_add 2.9168
σ_prop 0.12152
---------------------
```

We immediately see that the variances have decrease over the base model, but if we don’t remember the values we can summarize them using the `compare_estimates`

function.

`compare_estimates(base_model = iv1cmt_results, covariate_model = iv1cmt_covar_results)`

Row | parameter | base_model | covariate_model |
---|---|---|---|

String | Float64? | Float64? | |

1 | θcl | 0.41707 | 0.89494 |

2 | θvc | 7.16092 | 9.27417 |

3 | Ω₁,₁ | 0.171395 | 0.0562111 |

4 | Ω₂,₂ | 0.198171 | 0.0924213 |

5 | σ_add | 2.92824 | 2.91678 |

6 | σ_prop | 0.121131 | 0.121523 |

7 | θCLCR | missing | 0.889519 |

We should note several things here. First, the new parameter is at the bottom such that we can compare parameters with the same name on the same line. Second, the estimates of CL and Vc changed dramatically, but they cannot be compared directly, because they represent different quantities. For example, in the base model 0.417 L/hr is the population estimate for CL, while in the covariate model, 0.894 L/hr is the typical value of CL in a patient weighing 70 kg with a CLCR of 120 mL/min. Third, BSV decreased significantly! Fourth, the error model standard deviations basically didn’t change. This means that we have not really changed how much residual error is left, but we have managed to move some unexplained differences between subjects into explained differences through measured covariates. This is likely to improve the external validity or ability of the model to extrapolate under the assumption that we’ve identified the correct covariate model.

Let us look at the goodness of fit plots to verify that population predictions are now more in line with the observations.

`= inspect(iv1cmt_covar_results) iv1cmt_covar_inspect `

```
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.
```

```
FittedPumasModelInspection
Fitting was successful: true
Likehood approximation used for weighted residuals : Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_every = 1, time_limit = NaN, )
)
```

`goodness_of_fit(iv1cmt_covar_inspect)`

Additionally, we look at the correlations between covariates and EBEs.

`empirical_bayes_vs_covariates(iv1cmt_covar_inspect)`

Everything looks fine and we’ll accept our covariate model as our final model. To complete the analysis, we move on to calculating standard errors and relative standard errors (RSEs).

## 8 Inference

To calculate standard errors we use the `infer`

function.

```
= infer(iv1cmt_covar_results)
infer_covar = coeftable(infer_covar) infer_table
```

```
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
```

Row | parameter | estimate | se | ci_lower | ci_upper |
---|---|---|---|---|---|

String | Float64 | Float64 | Float64 | Float64 | |

1 | θcl | 0.89494 | 0.0591972 | 0.778916 | 1.01096 |

2 | θCLCR | 0.889519 | 0.0743712 | 0.743754 | 1.03528 |

3 | θvc | 9.27417 | 0.293658 | 8.69861 | 9.84973 |

4 | Ω₁,₁ | 0.0562111 | 0.00980985 | 0.0369841 | 0.075438 |

5 | Ω₂,₂ | 0.0924213 | 0.012244 | 0.0684234 | 0.116419 |

6 | σ_add | 2.91678 | 0.100189 | 2.72041 | 3.11314 |

7 | σ_prop | 0.121523 | 0.00657896 | 0.108629 | 0.134418 |

To calculate RSEs we can simply create a new column in the DataFrame.

```
= infer_table.se ./ infer_table.estimate .* 100
infer_table.RSE infer_table
```

Row | parameter | estimate | se | ci_lower | ci_upper | RSE |
---|---|---|---|---|---|---|

String | Float64 | Float64 | Float64 | Float64 | Float64 | |

1 | θcl | 0.89494 | 0.0591972 | 0.778916 | 1.01096 | 6.61466 |

2 | θCLCR | 0.889519 | 0.0743712 | 0.743754 | 1.03528 | 8.36084 |

3 | θvc | 9.27417 | 0.293658 | 8.69861 | 9.84973 | 3.16641 |

4 | Ω₁,₁ | 0.0562111 | 0.00980985 | 0.0369841 | 0.075438 | 17.4518 |

5 | Ω₂,₂ | 0.0924213 | 0.012244 | 0.0684234 | 0.116419 | 13.2481 |

6 | σ_add | 2.91678 | 0.100189 | 2.72041 | 3.11314 | 3.43491 |

7 | σ_prop | 0.121523 | 0.00657896 | 0.108629 | 0.134418 | 5.41375 |

## 9 Summarizing results

## 10 Conclusion

In this tutorial, we saw how to map tabular data to a `Population`

, how to formulate a base model and fit it, how to use built-in functionality to assess the quality of the formulated model, how to use the insight to build at covariate model, and how to finalize the analysis with diagnostic verification and inference on estimated parameters.