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

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/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

om <- omat(mod)
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 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 <- 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 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:

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.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

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.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

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.

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

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.5

We 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    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
 
$OMEGA
1 2 3
mod <- 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    3

We 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    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.

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>

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 IIV
1 2 3
$OMEGA @name IOV
4 5
mod <- 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 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

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>

6.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 5
mod <- 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