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
pkdata = CSV.read("CS1_IV1EST_PAR.csv", DataFrame; header = 4)
first(pkdata, 10)10×11 DataFrame
| 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 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 : :CONC3 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: CONCNotice, 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:
@paramis 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.@randomis used to define random effects and their corresponding distributions.@preis used to pre-process variables for the dynamic system and statistical specification.@dynamicsis used to define the dynamics system for the model. The differential equuation for this model would beCentral' = -CL/Vc*Central, but this system has a closed-form solution which Pumas implements usingCentral1. In general, you should use closed-form solutions likeCentral1when available since they come with a performance boost.@deriveddefines 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
endPumasModel
Parameters: θcl, θvc, Ω, σ_add, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: CONC
Observed: CONCBesides 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.07010698318481445
1 4.945610e+03 5.427910e+02
* time: 0.6771280765533447
2 4.870609e+03 4.624623e+02
* time: 0.7377200126647949
3 4.746284e+03 2.387450e+02
* time: 0.842108964920044
4 4.722871e+03 2.277545e+02
* time: 0.9394049644470215
5 4.697752e+03 2.695471e+02
* time: 1.0829501152038574
6 4.670279e+03 2.251210e+02
* time: 1.1757969856262207
7 4.653649e+03 7.845640e+01
* time: 1.2931931018829346
8 4.651567e+03 3.842645e+01
* time: 1.3638160228729248
9 4.650461e+03 3.919094e+01
* time: 1.403784990310669
10 4.649212e+03 3.992715e+01
* time: 1.4281740188598633
11 4.647847e+03 4.563022e+01
* time: 1.4533920288085938
12 4.644574e+03 4.320284e+01
* time: 1.4899060726165771
13 4.641136e+03 3.574904e+01
* time: 1.5154099464416504
14 4.639193e+03 3.808485e+01
* time: 1.5501139163970947
15 4.638535e+03 3.905231e+01
* time: 1.5737171173095703
16 4.637660e+03 4.078589e+01
* time: 1.6058070659637451
17 4.635157e+03 4.722472e+01
* time: 1.6298561096191406
18 4.620341e+03 1.466178e+02
* time: 1.654383897781372
19 4.597789e+03 3.472506e+02
* time: 1.6896231174468994
20 4.567135e+03 3.126339e+02
* time: 1.7113149166107178
21 4.491345e+03 1.786065e+02
* time: 1.7580299377441406
22 4.468647e+03 1.802908e+02
* time: 1.779642105102539
23 4.440016e+03 8.571120e+01
* time: 1.8011469841003418
24 4.429365e+03 9.173855e+01
* time: 1.831618070602417
25 4.413616e+03 5.147760e+01
* time: 1.8538250923156738
26 4.409551e+03 3.372222e+01
* time: 1.885369062423706
27 4.406285e+03 1.628709e+01
* time: 1.9071829319000244
28 4.406165e+03 1.665758e+01
* time: 1.9263310432434082
29 4.406084e+03 1.661883e+01
* time: 1.954474925994873
30 4.405484e+03 1.562746e+01
* time: 1.9739649295806885
31 4.404601e+03 1.315858e+01
* time: 1.9931390285491943
32 4.403168e+03 1.780770e+01
* time: 2.0232810974121094
33 4.402151e+03 1.436821e+01
* time: 2.0433619022369385
34 4.401710e+03 5.107809e+00
* time: 2.0618340969085693
35 4.401646e+03 8.160373e-01
* time: 2.0898399353027344
36 4.401643e+03 4.611581e-01
* time: 2.10862398147583
37 4.401643e+03 4.575021e-01
* time: 2.125473976135254
38 4.401643e+03 4.526048e-01
* time: 2.1510701179504395
39 4.401643e+03 4.350247e-01
* time: 2.1683549880981445
40 4.401642e+03 4.353443e-01
* time: 2.1851890087127686
41 4.401642e+03 7.252686e-01
* time: 2.2105789184570312
42 4.401641e+03 8.811155e-01
* time: 2.2278900146484375
43 4.401640e+03 6.601170e-01
* time: 2.2449910640716553
44 4.401639e+03 2.317088e-01
* time: 2.2707130908966064
45 4.401639e+03 2.594974e-02
* time: 2.288311004638672
46 4.401639e+03 1.773922e-03
* time: 2.304560899734497
47 4.401639e+03 1.773922e-03
* time: 2.353624105453491
48 4.401639e+03 1.773922e-03
* time: 2.405687093734741FittedPumasModel
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.639During 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.277989672784And 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.1609This 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.
coef(iv1cmt_results)which gives a named tuple of the coefficient valuescoeftable(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.19817The 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.12113Again, 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_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
endPumasModel
Parameters: θcl, θCLCR, θvc, Ω, σ_add, σ_prop
Random effects: η
Covariates: CLCR, WEIGHT
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: CONC
Observed: CONCand 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_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×5 DataFrame
| 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.RSE = infer_table.se ./ infer_table.estimate .* 100
infer_table7×6 DataFrame
| 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.