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

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

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

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

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()`.

## 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    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)
``````.   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`.

### 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.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()``````

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

### 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 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>``````

### 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 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``````