Working with SO files in R with libsoc

Navigating the SO structure

The SO R object closely follows the tree structure of the SO XML. We can extract the population parameter estimates like this:

param <- so$SOBlock[[1]]$Estimation$PopulationEstimates$MLE
param

##         TVCL     TVV     IVCL      IVV SIGMA_1_1_
## 1 0.00534888 1.42034 0.130798 0.177443  0.0188202

The result is a data.frame with the population estimates and the names of the parameters as column names. Note the way we go deeper in the SO tree with successive use of the $-opterator. As an SO can contain more than one SOBlock the SOBlock-structure is a list. In the example we simply select the first (and only) SOBlock.

libsoc also has a special function for extracting the parameter estimates across all SOBlocks

so$all_population_estimates()

##         TVCL     TVV     IVCL      IVV SIGMA_1_1_
## 1 0.00534888 1.42034 0.130798 0.177443  0.0188202

Getting information on parameters

If the SO is linked to a PharmML (for example if it was generated with “nmoutput2so -generate_pharmml”) parameter information can be retrieved from the ParameterModel.

cols <- colnames(param)
so$variability_type(cols)

##                   TVCL                    TVV                   IVCL 
##      "structParameter"      "structParameter" "parameterVariability" 
##                    IVV             SIGMA_1_1_ 
## "parameterVariability"        "residualError"

It is also possible to get information on which parameters represent correlations.

so$correlation_parameters(cols)

##       TVCL        TVV       IVCL        IVV SIGMA_1_1_ 
##      FALSE      FALSE      FALSE      FALSE      FALSE

This example has no correlation parameters.

The names of the random variables corresponding to the variability parameters can be retrieved:

so$random_variable_from_variability_parameter(cols)

##       TVCL        TVV       IVCL        IVV SIGMA_1_1_ 
##         NA         NA     "ETA1"     "ETA2"     "EPS1"

Obtaining some parameter uncertainty measures

The covariance matrix is of standard R matrix type:

so$SOBlock[[1]]$Estimation$PrecisionPopulationEstimates$MLE$CovarianceMatrix

##                    TVCL          TVV         IVCL          IVV   SIGMA_1_1_
## TVCL        1.45861e-07  1.14987e-05 -1.95288e-05  7.77462e-07 -1.29140e-07
## TVV         1.14987e-05  9.64296e-03 -1.98572e-03  3.38779e-03  1.45110e-04
## IVCL       -1.95288e-05 -1.98572e-03  7.65663e-03 -4.79990e-04  6.53171e-05
## IVV         7.77462e-07  3.38779e-03 -4.79990e-04  2.56664e-03  1.06472e-04
## SIGMA_1_1_ -1.29140e-07  1.45110e-04  6.53171e-05  1.06472e-04  1.52394e-05

The standard error for the parameters can be extracted for each SOBlock separately:

so$SOBlock[[1]]$Estimation$PrecisionPopulationEstimates$MLE$StandardError

##    Parameter          SE
## 1       TVCL 0.000381917
## 2        TVV 0.098198600
## 3       IVCL 0.087502200
## 4        IVV 0.050662000
## 5 SIGMA_1_1_ 0.003903770

Or more conveniently for all SOBlocks together:

so$all_standard_errors()

##          TVCL       TVV      IVCL      IVV SIGMA_1_1_
## 1 0.000381917 0.0981986 0.0875022 0.050662 0.00390377

Relative standard errors can currently only be retrieved for each SOBlock separately:

so$SOBlock[[1]]$Estimation$PrecisionPopulationEstimates$MLE$RelativeStandardError

##    Parameter       RSE
## 1       TVCL  7.140130
## 2        TVV  6.913739
## 3       IVCL 66.898729
## 4        IVV 28.551140
## 5 SIGMA_1_1_ 20.742447

Getting the OFV

so$SOBlock[[1]]$Estimation$OFMeasures$Deviance

## [1] 740.278

Getting the Predictions

pred <- so$SOBlock[[1]]$Estimation$Predictions
head(pred, 20)

##    ID  TIME   PRED  IPRED
## 1   1   0.0 17.601 17.704
## 2   1   2.0 17.469 17.575
## 3   1  12.5 19.256 19.395
## 4   1  24.5 20.870 21.045
## 5   1  37.0 22.374 22.588
## 6   1  48.0 23.930 24.181
## 7   1  60.5 25.294 25.584
## 8   1  72.5 26.641 26.970
## 9   1  85.3 27.851 28.222
## 10  1  96.5 29.165 29.573
## 11  1 108.5 30.341 30.789
## 12  1 112.5 29.887 30.344
## 13  2   0.0 10.561 11.212
## 14  2   2.0 10.482 11.140
## 15  2   4.0 13.078 13.908
## 16  2  16.0 15.176 16.218
## 17  2  27.8 17.192 18.450
## 18  2  40.0 19.095 20.575
## 19  2  52.0 20.927 22.631
## 20  2  63.5 20.040 21.803

Getting information on table columns

Each table in the SO structure has a columnType attribute that gives an array of columnTypes taken directly from the xml.

attributes(pred)$columnType

## [[1]]
## [1] "id"
## 
## [[2]]
## [1] "idv"
## 
## [[3]]
## [1] "undefined"
## 
## [[4]]
## [1] "undefined"

The portable way of finding the column name/number for id, idv and dv columns is to use the provided functions.

libsoc::id_column(pred)

## [1] 1

libsoc::id_column_name(pred)

## [1] "ID"

libsoc::idv_column(pred)

## [1] 2

libsoc::idv_column_name(pred)

## [1] "TIME"

Getting the Residuals

res <- so$SOBlock[[1]]$Estimation$Residuals$ResidualTable
head(res, 20)

##    ID  TIME       RES     IRES      WRES      IWRES     CWRES
## 1   1   2.0  -0.16929 -0.27518 -0.140800 -0.0156570 -0.142740
## 2   1 112.5   1.11300  0.65596  0.177620  0.0216180  0.176380
## 3   2   2.0  -0.78157 -1.43970 -1.175400 -0.1292400 -1.205800
## 4   2  63.5   4.56020  2.79730  0.973350  0.1283000  0.943370
## 5   2 135.5   3.94870  0.80411  0.335060  0.0249750  0.303240
## 6   3   1.5  -3.00270  1.01180  0.271450  0.0595560  0.255090
## 7   3  83.5  -5.01830  0.45785 -0.058389  0.0196150 -0.105980
## 8   3 134.3 -11.33000 -4.57310 -1.178000 -0.1583900 -1.230700
## 9   4   1.8   7.79310  0.24606  0.791600  0.0119710  0.638650
## 10  4  59.3   7.61640 -0.64012  0.407420 -0.0260850  0.286270
## 11  4 130.8  10.91000  1.57290  1.325300  0.0522080  1.213500
## 12  5   2.0  -4.66680  0.38435 -0.063200  0.0278200 -0.093037
## 13  5  59.5  -5.57280  0.31364 -0.159420  0.0175350 -0.205840
## 14  5 132.0 -10.09100 -3.11030 -1.059500 -0.1328600 -1.117200
## 15  6   1.8   2.21680  1.63500  0.742970  0.0941580  0.743250
## 16  6  59.3  -3.78550 -4.56150 -1.541700 -0.2086500 -1.543600
## 17  6 142.8   4.66410  3.58240  0.982730  0.1238800  0.979120
## 18  7   2.0   4.62330 -0.11309  0.447580 -0.0062783  0.344150
## 19  7  73.8   5.92770  1.17140  0.879840  0.0526960  0.842960
## 20  7 165.3   3.82140 -0.64356  0.080479 -0.0243370  0.067617

Reading messages

Messages from a run is stored in the TaskInformation structure. Each message contains a Name a Severity from 1 to 10, the message content and the name of the tool that emitted the message.

messages <- so$SOBlock[[1]]$TaskInformation$Message
messages[[1]]$Name

## [1] "estimation_successful"

messages[[1]]$Content

## [1] "1"

messages[[1]]$Severity

## [1] 1

messages[[1]]$Toolname

## [1] "NONMEM"

messages[[13]]$Content

## [1] "This SOBlock was created with nmoutput2so version 4.7.14"

iOFV

OFV values for each individual separately is in the OFMeasures section

iofv <- so$SOBlock[[1]]$Estimation$OFMeasures$IndividualContribToLL
head(iofv)

##   ID    ICtoLL
## 1  1  7.653419
## 2  2 11.528829
## 3  3 12.559465
## 4  4 11.696525
## 5  5 11.163075
## 6  6 13.959221