Case Study I: Building a Covariate Model

Author

Patrick Kofod Mogensen 

using Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
using CSV

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

The dataset is available from the PharmaDatasets package. 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.

pkdata = CSV.read(dataset("pumas/event_data/CS1_IV1EST_PAR", String), DataFrame)
first(pkdata, 10)
10×11 DataFrame
Row ID 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:

  • ID: 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, :WT => :WEIGHT, :ISM => :SEX)

Notice, that we rename the pkdata DataFrame column names in pairs of Symbols. The Symbols 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 Populations using the read_pumas function. In our case, the population is read using the following code:

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

Notice, that some inputs are listed as Symbols (a symbol with name “NAME” is written with a colon in front:NAME) and some are vectors of Symbols. 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 Symbols 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 RealDomains, 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.

iv1cmt = @model begin
    @metadata begin
        desc = "INTRAVENOUS BOLUS STUDY"
        timeu = u"hr" # hour
    end

    @param begin
        θcl ∈ RealDomain(lower = 0.1, init = 1.0)
        θvc ∈ RealDomain(lower = 1.0, upper = 20.0, init = 10.0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(lower = 0.0, init = 1.0^2)
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.09^2)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @pre begin
        CL = θcl * exp(η[1])
        Vc = θvc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
        conc_model := @. Central / Vc
        CONC ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model * σ_prop)^2))
    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 = (
    θcl = 0.9,
    θvc = 10.0,
    Ω = Diagonal([0.1, 0.1]),
    σ_add = sqrt(0.09),
    σ_prop = sqrt(0.01),
)

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.

iv1cmt_results = fit(iv1cmt, population, initial_est_iv1cmt, FOCE())
[ 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.023588180541992188
     1     4.945610e+03     5.427910e+02
 * time: 0.8952591419219971
     2     4.870609e+03     4.624623e+02
 * time: 1.0240471363067627
     3     4.746284e+03     2.387450e+02
 * time: 1.0929920673370361
     4     4.722871e+03     2.277545e+02
 * time: 1.1591191291809082
     5     4.697752e+03     2.695471e+02
 * time: 1.2197740077972412
     6     4.670279e+03     2.251210e+02
 * time: 1.3132569789886475
     7     4.653649e+03     7.845640e+01
 * time: 1.3691980838775635
     8     4.651567e+03     3.842645e+01
 * time: 1.4162371158599854
     9     4.650461e+03     3.919094e+01
 * time: 1.468696117401123
    10     4.649212e+03     3.992715e+01
 * time: 1.521078109741211
    11     4.647847e+03     4.563022e+01
 * time: 1.5793890953063965
    12     4.644574e+03     4.320284e+01
 * time: 1.6603360176086426
    13     4.641136e+03     3.574904e+01
 * time: 1.7143971920013428
    14     4.639193e+03     3.808485e+01
 * time: 1.765308141708374
    15     4.638535e+03     3.905231e+01
 * time: 1.8157060146331787
    16     4.637660e+03     4.078589e+01
 * time: 1.8714220523834229
    17     4.635157e+03     4.722472e+01
 * time: 1.9285809993743896
    18     4.620341e+03     1.466178e+02
 * time: 2.0025269985198975
    19     4.597789e+03     3.472506e+02
 * time: 2.0511929988861084
    20     4.567135e+03     3.126339e+02
 * time: 2.093830108642578
    21     4.491345e+03     1.786065e+02
 * time: 2.171344041824341
    22     4.468647e+03     1.802908e+02
 * time: 2.220393180847168
    23     4.440016e+03     8.571120e+01
 * time: 2.2681610584259033
    24     4.429365e+03     9.173855e+01
 * time: 2.3347151279449463
    25     4.413616e+03     5.147760e+01
 * time: 2.377368211746216
    26     4.409551e+03     3.372222e+01
 * time: 2.4180359840393066
    27     4.406285e+03     1.628709e+01
 * time: 2.4599971771240234
    28     4.406165e+03     1.665758e+01
 * time: 2.5035250186920166
    29     4.406084e+03     1.661883e+01
 * time: 2.5444161891937256
    30     4.405484e+03     1.562746e+01
 * time: 2.5877861976623535
    31     4.404601e+03     1.315858e+01
 * time: 2.6480090618133545
    32     4.403168e+03     1.780770e+01
 * time: 2.6877071857452393
    33     4.402151e+03     1.436821e+01
 * time: 2.7281761169433594
    34     4.401710e+03     5.107809e+00
 * time: 2.7653450965881348
    35     4.401646e+03     8.160373e-01
 * time: 2.8052830696105957
    36     4.401643e+03     4.611581e-01
 * time: 2.846492052078247
    37     4.401643e+03     4.575021e-01
 * time: 2.886507987976074
    38     4.401643e+03     4.526048e-01
 * time: 2.9442222118377686
    39     4.401643e+03     4.350247e-01
 * time: 2.9776041507720947
    40     4.401642e+03     4.353443e-01
 * time: 3.010867118835449
    41     4.401642e+03     7.252686e-01
 * time: 3.043915033340454
    42     4.401641e+03     8.811155e-01
 * time: 3.0764381885528564
    43     4.401640e+03     6.601170e-01
 * time: 3.113746166229248
    44     4.401639e+03     2.317088e-01
 * time: 3.1524240970611572
    45     4.401639e+03     2.594974e-02
 * time: 3.1900391578674316
    46     4.401639e+03     1.773922e-03
 * time: 3.2255501747131348
    47     4.401639e+03     1.773922e-03
 * time: 3.3259501457214355
    48     4.401639e+03     1.773922e-03
 * time: 3.403231143951416
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
Comparison to NONMEM

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.

Viewing results of the model fit

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

  1. coef(iv1cmt_results) which gives a named tuple of the coefficient values
  2. 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

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

For small variances, an approximation is often used

(CV_CL = sqrt(0.1714) * 100, CV_Vc = sqrt(0.19817) * 100)
(CV_CL = 41.400483088968905,
 CV_Vc = 44.51628915352222,)
Variances and standard deviations in the specification of random effects

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.

iv1cmt_inspect = inspect(iv1cmt_results)
[ 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_warnings = true, 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.

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

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.

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

iv1cmt_covar_results = fit(
    iv1cmt_covar,
    population,
    initial_est_iv1cmt_covar,
    Pumas.FOCE();
    optim_options = (show_trace = false,),
)
[ 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)
7×3 DataFrame
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.

iv1cmt_covar_inspect = inspect(iv1cmt_covar_results)
[ 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_warnings = true, 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_covar = infer(iv1cmt_covar_results)
infer_table = coeftable(infer_covar)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
7×6 DataFrame
Row parameter estimate se relative_se ci_lower ci_upper
String Float64 Float64 Float64 Float64 Float64
1 θcl 0.89494 0.0591972 0.0661466 0.778916 1.01096
2 θCLCR 0.889519 0.0743711 0.0836083 0.743754 1.03528
3 θvc 9.27417 0.293658 0.0316641 8.69861 9.84973
4 Ω₁,₁ 0.0562111 0.00981012 0.174523 0.0369836 0.0754385
5 Ω₂,₂ 0.0924213 0.012244 0.13248 0.0684235 0.116419
6 σ_add 2.91678 0.100189 0.0343491 2.72041 3.11314
7 σ_prop 0.121523 0.00657897 0.0541376 0.108629 0.134418

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

infer_table.RSE = infer_table.se ./ infer_table.estimate .* 100
infer_table
7×7 DataFrame
Row parameter estimate se relative_se ci_lower ci_upper RSE
String Float64 Float64 Float64 Float64 Float64 Float64
1 θcl 0.89494 0.0591972 0.0661466 0.778916 1.01096 6.61466
2 θCLCR 0.889519 0.0743711 0.0836083 0.743754 1.03528 8.36083
3 θvc 9.27417 0.293658 0.0316641 8.69861 9.84973 3.16641
4 Ω₁,₁ 0.0562111 0.00981012 0.174523 0.0369836 0.0754385 17.4523
5 Ω₂,₂ 0.0924213 0.012244 0.13248 0.0684235 0.116419 13.248
6 σ_add 2.91678 0.100189 0.0343491 2.72041 3.11314 3.43491
7 σ_prop 0.121523 0.00657897 0.0541376 0.108629 0.134418 5.41376

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.