<- modlib("popex", compile = FALSE) mod
6 Model Matrices
Model matrices include $OMEGA
(for subject-level variability) and $SIGMA
(for residual unexplained variability). These matrices are coded into the model file (see Section 2.2.17 and Section 2.2.18) and can be manipulated in different ways via the model object. Because, $OMEGA
and $SIGMA
matrices are handled using identical approach (only the names of the functions change), we will focus on working with $OMEGA
in the following examples with references to the equivalent functions that can be used to work on $SIGMA
. Also note that, for simplicity, we will not compile the examples presented in this chapter.
6.1 Basics
6.1.1 Simple matrix lists
We can look at the popex
model in the internal library for a starting example to show how model matrices can be seen from the model object.
Once the model is loaded
We can print the model object to the console and see the matrix structure
mod
.
.
. ---------------- source: popex.cpp ----------------
.
. project: /Users/kyleb/ren...gsolve/models
. shared object: popex-so-898d509d9a68 <not loaded>
.
. time: start: 0 end: 240 delta: 0.5
. add: <none>
.
. compartments: GUT CENT [2]
. parameters: TVKA TVCL TVV WT [4]
. captures: CL V ECL IPRED DV [5]
. omega: 3x3
. sigma: 1x1
.
. solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
. ------------------------------------------------------
In this output, we see that omega
is a 3 by 3 matrix and sigma
is 1 by 1. We can view both matrices by calling revar()
on the model object
revar(mod)
. $omega
. $...
. [,1] [,2] [,3]
. ECL: 0.3 0.0 0.0
. EV: 0.0 0.1 0.0
. EKA: 0.0 0.0 0.5
.
.
. $sigma
. $...
. [,1]
. 1: 0
This shows the 3x3 $OMEGA
matrix with all off-diagonals set to zero and the 1x1 $SIGMA
which is currently fixed to 0.
The $OMEGA
matrix can be extracted with the omat()
function
omat(mod)
. $...
. [,1] [,2] [,3]
. ECL: 0.3 0.0 0.0
. EV: 0.0 0.1 0.0
. EKA: 0.0 0.0 0.5
Use the smat()
function to extract the $SIGMA
matrix. The result of these calls are matlist
objects; for $OMEGA
the class is omegalist
(which inherits from matlist
)
omat(mod) %>% class()
. [1] "omegalist"
. attr(,"package")
. [1] "mrgsolve"
and for $SIGMA
it is sigmalist
. These are lists of matrices. In this example, there is just one $OMEGA
block in the code
blocks(mod, OMEGA)
.
. Model file: popex.cpp
.
. $OMEGA
. @labels ECL EV EKA
. 0.3 0.1 0.5
so the length of the omegalist
object is also 1
<- omat(mod)
om length(om)
. [1] 1
Functions are provided to check the names
names(om)
. [1] "..."
and the labels
labels(om)
. [[1]]
. [1] "ECL" "EV" "EKA"
as well as getting the dimensions or number of rows
dim(om)
. $...
. [1] 3 3
nrow(om)
. ...
. 3
The omegalist
(sigmalist
) object can be converted to a standard R list
as.list(om) %>% str()
. List of 1
. $ ...: num [1:3, 1:3] 0.3 0 0 0 0.1 0 0 0 0.5
. ..- attr(*, "dimnames")=List of 2
. .. ..$ : chr [1:3] "ECL" "EV" "EKA"
. .. ..$ : chr [1:3] "ECL" "EV" "EKA"
or it can be rendered as a matrix
as.matrix(om)
. [,1] [,2] [,3]
. [1,] 0.3 0.0 0.0
. [2,] 0.0 0.1 0.0
. [3,] 0.0 0.0 0.5
6.1.2 Multiple matrix lists
Let’s look at an example where there is a more complicated $OMEGA
structure.
// inline/multiple-matrices.cpp
$OMEGA @name first @labels a b c1 2 3
$OMEGA @name second @labels d e4 5
$OMEGA @name third @labels f g h i6 7 8 9
<- mread("inline/multiple-matrices.cpp", quiet = TRUE) mod
Each $OMEGA
block codes a diagonal matrix; the interesting feature is that there are 3 different $OMEGA
blocks.
Now, when we look at the model
mod
.
.
. ---------- source: multiple-matrices.cpp ----------
.
. project: /Users/kyleb/git...-guide/inline
. shared object: multiple-matrices.cpp-so-898d3b868a4
.
. time: start: 0 end: 24 delta: 1
. add: <none>
.
. compartments: <none>
. parameters: <none>
. captures: <none>
. omega: 3x3,2x2,4x4
. sigma: 0x0
.
. solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
. ------------------------------------------------------
We see that $OMEGA
is one 3x3 matrix, one 2x2 matrix and on 4x4 matrix, in that order. Calling revar()
on this model object
revar(mod)
. $omega
. $first
. [,1] [,2] [,3]
. a: 1 0 0
. b: 0 2 0
. c: 0 0 3
.
. $second
. [,1] [,2]
. d: 4 0
. e: 0 5
.
. $third
. [,1] [,2] [,3] [,4]
. f: 6 0 0 0
. g: 0 7 0 0
. h: 0 0 8 0
. i: 0 0 0 9
.
.
. $sigma
. No matrices found
we only see the $OMEGA
matrices as the model was coded. Now the length of the omegalist
object is 3
length(omat(mod))
. [1] 3
and the number of rows is 3 in the first matrix, 2 in the second, and 4 in the third
nrow(omat(mod))
. first second third
. 3 2 4
We can also check dim()
dim(omat(mod))
. $first
. [1] 3 3
.
. $second
. [1] 2 2
.
. $third
. [1] 4 4
The omega matrix can be converted from this segmented list into a single block matrix
as.matrix(omat(mod))
. [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
. [1,] 1 0 0 0 0 0 0 0 0
. [2,] 0 2 0 0 0 0 0 0 0
. [3,] 0 0 3 0 0 0 0 0 0
. [4,] 0 0 0 4 0 0 0 0 0
. [5,] 0 0 0 0 5 0 0 0 0
. [6,] 0 0 0 0 0 6 0 0 0
. [7,] 0 0 0 0 0 0 7 0 0
. [8,] 0 0 0 0 0 0 0 8 0
. [9,] 0 0 0 0 0 0 0 0 9
The result is 9x9 with all off diagonals between the different list positions set to zero.
Otherwise, we might work with this object as a list
as.list(omat(mod))
. $first
. a b c
. a 1 0 0
. b 0 2 0
. c 0 0 3
.
. $second
. d e
. d 4 0
. e 0 5
.
. $third
. f g h i
. f 6 0 0 0
. g 0 7 0 0
. h 0 0 8 0
. i 0 0 0 9
6.2 Collapsing matrices
The functionality described in this subsection is new with mrgsolve 1.0.0.
The structure of the matrices and the order in which they appear in the list matter when updating each matrix (see below). We saw how to create one big matrix out of all the smaller matrices using the as.matrix()
. This section describes how to combine matrices within the confines of the matlist
object.
mrgsolve provides functions to collapse (or combine) matrices in a matlist
object. We can call collapse_omega()
, passing in the model object
collapse_omega(mod) %>% revar()
. $omega
. $...
. [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
. a: 1 0 0 0 0 0 0 0 0
. b: 0 2 0 0 0 0 0 0 0
. c: 0 0 3 0 0 0 0 0 0
. d: 0 0 0 4 0 0 0 0 0
. e: 0 0 0 0 5 0 0 0 0
. f: 0 0 0 0 0 6 0 0 0
. g: 0 0 0 0 0 0 7 0 0
. h: 0 0 0 0 0 0 0 8 0
. i: 0 0 0 0 0 0 0 0 9
.
.
. $sigma
. No matrices found
and now we have an omegalist
object inside a model object and the matrix is a single 9x9 $OMEGA
matrix. The row names have been retained, but there is now no name for the matrix; this can be provided when collapsing
collapse_omega(mod, name = "only") %>% omat()
. $only
. [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
. a: 1 0 0 0 0 0 0 0 0
. b: 0 2 0 0 0 0 0 0 0
. c: 0 0 3 0 0 0 0 0 0
. d: 0 0 0 4 0 0 0 0 0
. e: 0 0 0 0 5 0 0 0 0
. f: 0 0 0 0 0 6 0 0 0
. g: 0 0 0 0 0 0 7 0 0
. h: 0 0 0 0 0 0 0 8 0
. i: 0 0 0 0 0 0 0 0 9
Suppose that we only want to combine the first two matrices, leaving the third matrix alone. In that case, call collapse_omega()
with the range argument
collapse_omega(mod, range = c(1,2), name = "first_second") %>% omat()
. $first_second
. [,1] [,2] [,3] [,4] [,5]
. a: 1 0 0 0 0
. b: 0 2 0 0 0
. c: 0 0 3 0 0
. d: 0 0 0 4 0
. e: 0 0 0 0 5
.
. $third
. f g h i
. f: 6 0 0 0
. g: 0 7 0 0
. h: 0 0 8 0
. i: 0 0 0 9
Now the matlist topology has changed; there are still 9 (total) rows and columns, but the matlist object is length 2 with 5x5 in the first position (newly named first_second
) an the old 4x4 matrix in the second position. Collapsing matrices is an irreversible process; at this time there is no mechanism to cut matrices back into smaller chunks. But collapsing matrices can be very helpful when they need to be updated.
Use collapse_sigma()
to collapse $SIGMA
matrices if needed.
Also, the function collapse_matrix()
can be called on a omegalist
or sigmalist
object to collapse
omat(mod) %>% collapse_matrix(range = c(2,NA))
. $first
. a b c
. a: 1 0 0
. b: 0 2 0
. c: 0 0 3
.
. $...
. [,1] [,2] [,3] [,4] [,5] [,6]
. d: 4 0 0 0 0 0
. e: 0 5 0 0 0 0
. f: 0 0 6 0 0 0
. g: 0 0 0 7 0 0
. h: 0 0 0 0 8 0
. i: 0 0 0 0 0 9
6.2.1 Collapse on load
You can also collapse $OMEGA
or $SIGMA
from within the model specification file so that the matrices are collapsed on load.
In the $SET
set block (Section 2.2.19), pass either collapse_omega
or collapse_sigma
set to TRUE
$SET collapse_omega = TRUE
or
$SET collapse_sigma = TRUE
This will produce a model object with collapsed matrices when the model is loaded via mread()
.
6.3 Updating $OMEGA
and $SIGMA
Like the values of parameters in the parameter list, we may want to update the values in $OMEGA
and $SIGMA
matrices. We can do so without re-compiling the model.
6.3.1 Matrix helper functions
mrgsolve
keeps $OMEGA
and $SIGMA
in block matrices (regardless of whether the off-diagonal elements are zeros or not). Recall that in the model specification file we can enter data for $OMEGA
and $SIGMA
as the lower triangle of the matrix (see Section 2.2.17). In R
, we need to provide a matrix (as an R
object). mrgsolve
provides some convenience functions to help … allowing the user to enter lower diagonals instead of the full matrix.
dmat()
for diagonal matrix
dmat(1,2,3)
. [,1] [,2] [,3]
. [1,] 1 0 0
. [2,] 0 2 0
. [3,] 0 0 3
bmat()
for block matrix
bmat(1,2,3)
. [,1] [,2]
. [1,] 1 2
. [2,] 2 3
cmat()
for a block matrix where the diagonal elements are variances and the off-diagonals are taken to be correlations, not covariances
cmat(0.1, 0.87,0.3)
. [,1] [,2]
. [1,] 0.1000000 0.1506884
. [2,] 0.1506884 0.3000000
mrgsolve
will convert the correlations to covariances.
mrgsolve
also provides as_bmat()
and as_dmat()
for converting other R
objects to matrices or lists of matrices.
Consider this list with named elements holding the data for a matrix:
<- list(OMEGA1.1 = 0.9, OMEGA2.1 = 0.3, OMEGA2.2 = 0.4) m
These data could form either a 3x3 diagonal matrix or a 2x2 block matrix. But the names suggest a 2x2 form. as_bmat()
can make the matrix like this
as_bmat(m, "OMEGA")
. [,1] [,2]
. [1,] 0.9 0.3
. [2,] 0.3 0.4
The second argument is a regular expression that mrgsolve
uses to find elements in the list to use for building the matrix.
Frequently, we have estimates in a data frame like this
data(exBoot)
head(exBoot)
. run THETA1 THETA2 THETA3 OMEGA11 OMEGA21 OMEGA22 OMEGA31 OMEGA32 OMEGA33
. 1 1 -0.7634 2.280 0.8472 0.12860 0.046130 0.2874 0.13820 -0.02164 0.3933
. 2 2 -0.4816 2.076 0.5355 0.12000 0.051000 0.2409 0.06754 -0.07759 0.3342
. 3 3 -0.5865 2.334 -0.4597 0.11460 0.097150 0.2130 0.16650 0.18100 0.4699
. 4 4 -0.6881 1.824 0.7736 0.14990 0.000003 0.2738 0.24700 -0.05466 0.5536
. 5 5 0.2909 1.519 -1.2440 0.07308 0.003842 0.2989 0.06475 0.05078 0.2500
. 6 6 0.1135 2.144 -1.0040 0.13390 -0.019270 0.1640 0.10740 -0.01170 0.3412
. SIGMA11 SIGMA21 SIGMA22
. 1 0.002579 0 1.0300
. 2 0.002228 0 1.0050
. 3 0.002418 0 1.0890
. 4 0.002177 0 0.8684
. 5 0.001606 0 0.8996
. 6 0.002134 0 0.9744
We can use as_bmat()
with this data frame to extract the $OMEGA
matrices
<- as_bmat(exBoot, "OMEGA")
omegas length(omegas)
. [1] 100
dim(exBoot)
. [1] 100 13
6]] omegas[[
. [,1] [,2] [,3]
. [1,] 0.13390 -0.01927 0.1074
. [2,] -0.01927 0.16400 -0.0117
. [3,] 0.10740 -0.01170 0.3412
16]] omegas[[
. [,1] [,2] [,3]
. [1,] 0.08126 0.01252 0.1050
. [2,] 0.01252 0.16860 0.0149
. [3,] 0.10500 0.01490 0.4062
The result of calling as_bmat
or as_dmat
is a list of matrices, one for each row in the data frame.
Note in this example, we could have called
<- as_bmat(exBoot,"SIGMA") sigmas
to grab the $SIGMA
matrices.
For help on these helper functions, see ?dmat
, ?bmat
, ?cmat
, ?as_bmat
, ?as_dmat
in the R
help system after loading mrgsolve
.
6.3.2 Fill a matrix with zeros
Sometimes we write a population model that includes random effects, but we would like to simulate from that same model without the random effects implemented. For example, we want to simulate some typical PK profiles from a population PK model that includes IIV on some parameters and / or RUV on the simulated outputs.
To do this, pass the model through the zero_re()
function. By default, this will convert all $OMEGA
and $SIGMA
matrix elements to zeros. See the R help file (?zero_re
) to see some options for selectively zeroing out only one or the other.
For example we have this population PK model
<- modlib("popex", compile = FALSE)
mod
omat(mod)
. $...
. [,1] [,2] [,3]
. ECL: 0.3 0.0 0.0
. EV: 0.0 0.1 0.0
. EKA: 0.0 0.0 0.5
We can turn that matrix to all zeros with
%>% zero_re() %>% omat() mod
. $...
. [,1] [,2] [,3]
. ECL: 0 0 0
. EV: 0 0 0
. EKA: 0 0 0
And when we simulate right after that, all ETA(n)
will be zero as well and you’ll get your fixed-effects simulation (the following is for example only and is not evaluated)
%>%
mod zero_re() %>%
ev(amt = 100) %>%
mrgsim() %>%
plot()
6.3.3 Example: unnamed matrix
Here is a model with only a 3x3 $OMEGA
matrix
// inline/matrix.cpp
$OMEGA1 2 3
<- mread("inline/matrix.cpp", compile = FALSE, quiet = TRUE) mod
Let’s check the values in the matrix using omat()
%>% omat mod
. $...
. [,1] [,2] [,3]
. 1: 1 0 0
. 2: 0 2 0
. 3: 0 0 3
We also use omat()
to update the values in the matrix
%>% omat(dmat(4,5,6)) %>% omat mod
. $...
. [,1] [,2] [,3]
. 1: 4 0 0
. 2: 0 5 0
. 3: 0 0 6
To update $OMEGA
, we must provide a matrix of the same dimension, in this case 3x3. An error is generated if we provide a matrix with the wrong dimension.
<- try(mod %>% omat(dmat(11,23))) ans
. Error : improper signature: omat
ans
. [1] "Error : improper signature: omat\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: improper signature: omat>
6.3.4 Example: named matrices
When there are multiple $OMEGA
matrices, it can be helpful to assign them names. Here, there are two matrices: one for interindividual variability (IIV
) and one for interoccasion variability (IOV
).
// inline/iov.cpp
$OMEGA @name IIV1 2 3
$OMEGA @name IOV4 5
<- mread("inline/iov.cpp", compile = FALSE, quiet = TRUE)
mod
revar(mod)
. $omega
. $IIV
. [,1] [,2] [,3]
. 1: 1 0 0
. 2: 0 2 0
. 3: 0 0 3
.
. $IOV
. [,1] [,2]
. 4: 4 0
. 5: 0 5
.
.
. $sigma
. No matrices found
Now, we can update either IIV
or IOV
(or both) by name
%>%
mod omat(IOV = dmat(11,12), IIV = dmat(13, 14, 15)) %>%
omat()
. $IIV
. [,1] [,2] [,3]
. 1: 13 0 0
. 2: 0 14 0
. 3: 0 0 15
.
. $IOV
. [,1] [,2]
. 4: 11 0
. 5: 0 12
Again, an error is generated if we try to assign a 3x3 matrix to the IOV
position
<- try(mod %>% omat(IIV = dmat(1, 2))) ans
. Error : improper dimension: omat
ans
. [1] "Error : improper dimension: omat\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: improper dimension: omat>
6.3.5 Example: unnamed matrices
If we do write the model with unnamed matrices, we can still update them
// inline/multi.cpp
$OMEGA1 2 3
$OMEGA 4 5
<- mread("inline/multi.cpp", compile = FALSE, quiet = TRUE) mod
In this case, the only way to update is to pass in a list of matrices, where (in this example) the first matrix is 3x3 and the second is 2x2
%>% omat(list(dmat(5, 6, 7), dmat(8, 9))) %>% omat() mod
. $...
. [,1] [,2] [,3]
. 1: 5 0 0
. 2: 0 6 0
. 3: 0 0 7
.
. $...
. [,1] [,2]
. 4: 8 0
. 5: 0 9