mod <- modlib("popex", compile = FALSE)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 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.
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
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: 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 3nrow(om). ...
. 3The 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.55.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 c
1 2 3
$OMEGA @name second @labels d e
4 5
$OMEGA @name third @labels f g h i
6 7 8 9mod <- mread("inline/multiple-matrices.cpp", quiet = TRUE)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 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 95.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 = TRUEor
$SET collapse_sigma = TRUEThis will produce a model object with collapsed matrices when the model is loaded via mread().
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.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 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
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.9744We can use as_bmat() with this data frame to extract the $OMEGA matrices
omegas <- as_bmat(exBoot, "OMEGA")
length(omegas). [1] 100dim(exBoot). [1] 100 13omegas[[6]]. [,1] [,2] [,3]
. [1,] 0.13390 -0.01927 0.1074
. [2,] -0.01927 0.16400 -0.0117
. [3,] 0.10740 -0.01170 0.3412omegas[[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
// inline/matrix.cpp
$OMEGA
1 2 3mod <- mread("inline/matrix.cpp", compile = FALSE, quiet = TRUE)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: omatans. [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).
// inline/iov.cpp
$OMEGA @name IIV
1 2 3
$OMEGA @name IOV
4 5mod <- mread("inline/iov.cpp", compile = FALSE, quiet = TRUE)
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: omatans. [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
// inline/multi.cpp
$OMEGA
1 2 3
$OMEGA
4 5mod <- mread("inline/multi.cpp", compile = FALSE, quiet = TRUE)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