5 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 Sections 2.2.16 and 2.2.17) 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.
5.1 Basics
5.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
mod <- modlib("popex", compile = FALSE)We can print the model object to the console and see the matrix structure
mod.
.
. ---------------- source: popex.cpp ----------------
.
. project: /Users/kyleb/Rli...gsolve/models
. shared object: popex-so-e41611940d1e <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: 0This 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.5Use 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.5so the length of the omegalist object is also 1
om <- omat(mod)
length(om). [1] 1Functions 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). ...
. 3The omegalist (sigmalist) object can be converted to a standard R list
. 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.55.1.2 Multiple matrix lists
Let’s look at an example where there is a more complicated $OMEGA structure.
code <- '
$OMEGA @name first @labels a b c
1 2 3
$OMEGA @name second @labels d e
4 5
$OMEGA @name third @labels f g h i
6 7 8 9
'
mod <- mcode("multiple-matrices", code, compile = FALSE)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: /private/var/fol.../T/RtmpxbKKp1
. shared object: multiple-matrices-so-e4161e84987b <not loaded>
.
. 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 foundwe only see the $OMEGA matrices as the model was coded. Now the length
of the omegalist object is 3
length(omat(mod)). [1] 3and 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 4We can also check dim()
dim(omat(mod)). $first
. [1] 3 3
.
. $second
. [1] 2 2
.
. $third
. [1] 4 4The 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 9The 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 95.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 foundand 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 9Suppose 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 9Now 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
5.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.
5.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.16). 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 3bmat() for block matrix
bmat(1,2,3). [,1] [,2]
. [1,] 1 2
. [2,] 2 3cmat() 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.3000000mrgsolve 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:
m <- list(OMEGA1.1 = 0.9, OMEGA2.1 = 0.3, OMEGA2.2 = 0.4)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.4The 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
. 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.9744We can use as_bmat() with this data frame to extract the $OMEGA matrices
omegas <- as_bmat(exBoot, "OMEGA")
length(omegas). [1] 100
dim(exBoot). [1] 100 13
omegas[[6]]. [,1] [,2] [,3]
. [1,] 0.13390 -0.01927 0.1074
. [2,] -0.01927 0.16400 -0.0117
. [3,] 0.10740 -0.01170 0.3412
omegas[[16]]. [,1] [,2] [,3]
. [1,] 0.08126 0.01252 0.1050
. [2,] 0.01252 0.16860 0.0149
. [3,] 0.10500 0.01490 0.4062The 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
sigmas <- as_bmat(exBoot,"SIGMA") 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.
5.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
mod <- modlib("popex", compile = FALSE)
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.5We can turn that matrix to all zeros with
mod %>% zero_re() %>% omat(). $...
. [,1] [,2] [,3]
. ECL: 0 0 0
. EV: 0 0 0
. EKA: 0 0 0And 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()5.3.3 Example: unnamed matrix
Here is a model with only a 3x3 $OMEGA matrix
code <- '
$OMEGA
1 2 3
'
mod <- mcode("matrix", code, compile=FALSE)Let’s check the values in the matrix using omat()
mod %>% omat. $...
. [,1] [,2] [,3]
. 1: 1 0 0
. 2: 0 2 0
. 3: 0 0 3We also use omat() to update the values in the matrix
mod %>% omat(dmat(4,5,6)) %>% omat. $...
. [,1] [,2] [,3]
. 1: 4 0 0
. 2: 0 5 0
. 3: 0 0 6To 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.
ans <- try(mod %>% omat(dmat(11,23))). Error : improper signature: omat
ans. [1] "Error : improper signature: omat\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: improper signature: omat>5.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).
code <- '
$OMEGA @name IIV
1 2 3
$OMEGA @name IOV
4 5
'
mod <- mcode("iov", code, compile=FALSE)
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 foundNow, 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 12Again, an error is generated if we try to assign a 3x3 matrix to the IOV
position
ans <- try(mod %>% omat(IIV = dmat(1, 2))). Error : improper dimension: omat
ans. [1] "Error : improper dimension: omat\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: improper dimension: omat>5.3.5 Example: unnamed matrices
If we do write the model with unnamed matrices, we can still update them
code <- '
$OMEGA
1 2 3
$OMEGA
4 5
'
mod <- mcode("multi", code, compile = FALSE)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
mod %>% omat(list(dmat(5, 6, 7), dmat(8, 9))) %>% omat(). $...
. [,1] [,2] [,3]
. 1: 5 0 0
. 2: 0 6 0
. 3: 0 0 7
.
. $...
. [,1] [,2]
. 4: 8 0
. 5: 0 9