---
title: "Get started"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Get started}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## Installation
```{r, eval=FALSE}
# Stable release
install.packages("SimuRg")
# Development version:
remotes::install_github("ms-decisions/SimuRg")
```
```{r setup}
library(SimuRg)
```
This package is designed for pharmacometricians and provides a complete workflow for building, fitting, and evaluating population PK/PD models. Model fitting, the first step of the workflow, is performed using the `sg_fit()` function. Currently, two fitting engines are supported: Monolix 2023 and Simurg cybernetic core. Both engines must be installed separately, as they are not distributed with the package. The fitting engine is selected using the opt_name argument.
## Model calibration
To fit a model, several mandatory arguments must be provided. The first is `model`, which specifies the path to the model file. When using the Monolix fitting engine (`opt_name = "Monolix"`, the default), the model must be written in the mlxtran format. When using the SimuRg Core engine (`opt_name = "Simurg"`), both rxode2 and mlxtran model formats are supported. In the following example, mlxtran format is used.
```{r model}
library(tibble)
library(dplyr)
library(stringr)
library(readr)
model <- system.file("extdata", "models", "model_PK_1c.txt", package = "SimuRg")
read_lines(model)
```
### Data specification
Next, the `data` argument must be used to specify the path to the input dataset. The dataset should be provided in an ADPPK-like format.
The structure of the dataset must be described using the headers argument, which is a list of column specifications. Each column specification is itself a list with the following elements:
- `name` - character string specifying the name of the dataset column;
- `use` - character string defining the role of the column according to the Monolix data format specification (e.g., `"id"`, `"time"`, `"observation"`, `"covariate"`);
- `type` - character string specifying the covariate type. For columns with `use = "covariate"`, this must be either "continuous" or "categorical". For all other column types, this field should be `NULL`.
```{r data}
data <- system.file("extdata", "datasets", "dspk-warf.csv", package = "SimuRg")
data_content <- read.csv(data)
head(data)
headers <- list(list(name = "ID", use = "identifier", type = NULL),
list(name = "TIME", use = "time", type = NULL),
list(name = "DV", use = "observation", type = "continuous"),
list(name = "DVID", use = "observationtype", type = NULL),
list(name = "ADM", use = "administration", type = NULL),
list(name = "AMT", use = "amount", type = NULL),
list(name = "EVID", use = "eventidentifier", type = NULL),
list(name = "MDV", use = "missingdependentvariable", type = NULL),
list(name = "AGE", use = "covariate", type = "continuous"),
list(name = "AGE_centered", use = "covariate", type = "continuous"),
list(name = "SEX", use = "covariate", type = "categorical"),
list(name = "WEIGHT", use = "covariate", type = "continuous"),
list(name = "BMI", use = "covariate", type = "continuous"),
list(name = "CLCR", use = "covariate", type = "continuous"),
list(name = "CYP2C9_gentyp", use = "covariate", type = "categorical"),
list(name = "VKORC1_gentyp", use = "covariate", type = "categorical"),
list(name = "G1_1", use = "ignore", type = NULL),
list(name = "G1_2", use = "ignore", type = NULL),
list(name = "G1_3", use = "ignore", type = NULL),
list(name = "G2_2", use = "ignore", type = NULL),
list(name = "G2_3", use = "ignore", type = NULL),
list(name = "G3_3", use = "ignore", type = NULL),
list(name = "GG", use = "ignore", type = NULL),
list(name = "AG", use = "ignore", type = NULL),
list(name = "AA", use = "ignore", type = NULL))
```
### Statistical components
After specifying the model and dataset, the statistical components of the model must be defined. The first is the `theta` argument, which is a data frame describing the model parameters and their estimation settings. It must contain the following columns:
- `NAME` – character string specifying the parameter name;
- `TRANS` – character string specifying the parameter distribution. Supported values are `"normal"`, `"logNormal"`, and `"logitNormal"`;
- `INIT` – numeric value specifying the initial estimate or, for fixed parameters, the fixed value;
- `LB` – numeric value specifying the lower bound for parameters with a `"logitNormal"` distribution. For all other distributions, this value should be `NA`;
- `UB` – numeric value specifying the upper bound for parameters with a `"logitNormal"` distribution. For all other distributions, this value should be `NA`.
- `EST` – logical value indicating whether the parameter should be estimated (`TRUE`) or fixed (`FALSE`).
In our example, all parameters will be estimated and will have the lognormal distribution
```{r theta}
theta <- tribble(~NAME, ~TRANS, ~INIT, ~LB, ~UB, ~EST,
"Cl", "logNormal", 0.2, NA, NA, TRUE,
"V", "logNormal", 20, NA, NA, TRUE,
"ka", "logNormal", 0.2, NA, NA, TRUE
)
```
The second statistical component defines the random effects and is specified using the `re` argument. This object consists of two square matrices with dimensions equal to the number of model parameters:
- `init` – specifies the initial values of the random-effects covariance matrix.
- `est` – specifies how each element of the covariance matrix is treated during estimation:
- `TRUE` – the corresponding element is estimated.
- `FALSE` – the corresponding element is fixed at its initial value.
- `NA` – the corresponding random effect is omitted from the model.
The rows and columns of both matrices correspond to the model parameters defined in the `theta` data frame.
Between-occasion variability, defined by `occ` parameter, is specified in the same way, as the `re` parameter.
In our case, we will add the between subjects variability to `Cl` and `ka` parameters. No between-occasion variability will be added.
```{r re}
re <- list(init = tribble(~Cl, ~V, ~ka,
1, 0, 0,
0, 0, 0,
0, 0, 1) %>% as.matrix(),
est = tribble(~Cl, ~V, ~ka,
TRUE, NA, NA,
NA, NA, NA,
NA, NA, TRUE) %>% as.matrix())
occ <- list(init = tribble(~Cl, ~V, ~ka,
0, 0, 0,
0, 0, 0,
0, 0, 0) %>% as.matrix(),
est = tribble(~Cl, ~V, ~ka,
NA, NA, NA,
NA, NA, NA,
NA, NA, NA) %>% as.matrix())
```
The last statistical component specifies the residual unexplained variability (RUV) model and is provided through the ruv argument. This object is a list of observation-specific specifications, where each element is itself a list with the following fields:
- `YNAME` – character string specifying the observation name (typically "y1", "y2", etc.);
- `DVID` – numeric identifier of the observation type, corresponding to the values in the DVID column of the dataset;
- `TRANS` – character string specifying the residual error distribution. Supported values are "normal", "logNormal", and "logitNormal";
- `PRED` – character string specifying the name of the prediction variable defined in the model;
- `ERR` – character string specifying the residual error model. Supported values are "constant" (additive error), "proportional" (proportional error), and "combined1" (combined additive and proportional error);
- `INIT` – numeric vector containing the initial values of the residual error parameters. The required length depends on the selected error model;
- `EST` – logical vector indicating whether each residual error parameter should be estimated (`TRUE`) or fixed (`FALSE`). This vector must have the same length as INIT;
- `BLQM` – below-limit-of-quantification (BLQ) handling method. If no BLQ method is used, this field should be NULL.
In our example, we have one DVID, therefore ruv object will have the following structure:
```{r ruv}
ruv <- list(YNAME = "y1", DVID = 1, TRANS = "normal", PRED = "Cc",
ERR = "combined1", INIT = c(1, 1), EST = c(TRUE, TRUE), BLQM = NULL)
```
### Covariates
Finally, covariate effects can be specified using the `cov` argument. This argument is a list of covariate specifications, where each element is itself a list with the following fields:
- `PAR` – character string specifying the name of the model parameter to which the covariate effect is applied.
- `COVNAME` – character string specifying the name of the covariate.
- `FUNC` – character string specifying the functional form of the covariate effect. Supported values are `"linear"` for continuous covariates and `"categorical"` for categorical covariates.
- `TRANS` – character string specifying the covariate transformation. Supported values are `"median"` for continuous covariates and `"reference"` for categorical covariates.
- `INIT` – numeric value specifying the initial estimate of the covariate effect.
- `EST` – logical value indicating whether the covariate effect should be estimated (`TRUE`) or fixed (`FALSE`).
```{r cov}
covs <- list(list(PAR = "V", COVNAME = "AGE", FUNC = "linear",
TRANS = "median", INIT = 1, EST = TRUE),
list(PAR = "ka", COVNAME = "SEX", REF = 0, INIT = 1, EST = TRUE))
```
### Model calibration
The `sg_fit()` function also provides several optional arguments:
- `project_name` – character string specifying the project name. This name is used for the generated `.mlxtran` file and the output directory containing the fitting results.
- `fit` – logical value indicating whether model fitting should be performed. If `FALSE`, only the control object (either a `.mlxtran` project or a `GCO` object, depending on the selected engine) is generated.
- `path_to_save_output` – path to the directory where the output files will be saved.
- `path_to_fitter` – path to the fitting engine executable.
- `max_wait_time` – maximum time (in seconds) to wait for the fitting process to complete.
The model can now be fitted by calling `sg_fit()`. Since neither Monolix nor Simurg core is available in the environment used to build this vignette, the fitting step is skipped by setting `fit = FALSE`. In this case, only the control object is generated. The resulting control object can subsequently be submitted to the corresponding fitting engine in a separate environment where the required software is installed.
```{r fit}
output_path <- str_c(tempdir(), "/")
task_opt <- paste("populationParameters()", "individualParameters()",
"logLikelihood()", sep = "\n")
result <- sg_fit(model, data, headers, theta, ruv, re, occ, covs,
project_name = "my_project", fit = FALSE, # set fit = TRUE for fit
path_to_save_output = output_path)
```
### Reading model calibration results
Another option of SimuRg package is not to fit a model inside the package, but just read Monolix/Simurg fit results. For this goal, `sg_converter()` function can be used. Its use and specification is simpler, than the use of `sg_fit()`, as it has only three arguments:
- folder_path - character string with path to the directory with Monolix project files
- proj_name - character string with the name of the project
```{r read}
test_folder <- system.file("extdata", "Monolix_objects", package = "SimuRg")
if (substr(test_folder, nchar(test_folder), nchar(test_folder)) != "/")
test_folder <- str_c(test_folder, "/")
pro_name <- "proj-solo"
message("Resolved folder: ", test_folder)
message("Folder exists: ", dir.exists(test_folder))
result <- sg_converter(folder_path = test_folder, proj_name = pro_name)
```
`sg_converter()` and `sg_fit()` objects return `GCO` and `GFO` objects, which are inputs for other functions of the package. Please check the documentation to see the.
## Goodness-of-fit
### Basic GoF plots
Now, when the model was calibrated, one can do model diagnostics. SimuRg package include following functions for model goodness of fit diagnostics:
* `sg_gof_obpr()` - observed versus predicted plot;
* `sg_gof_tp()` - time profiles visualization;
* `sg_gof_par_cov()` - visualization of random effect/individual parameters vs covariated;
* `sg_gof_par_dist()` - plot the distribution of random effects/individual parameters;
* `sg_gof_res_dist()` - plot the distribution of the residual
* `sg_gof_res()` - create residual diagnostic plots
```{r gof}
sg_gof_obpr(result$GFO)
sg_gof_tp(result$GFO)
sg_gof_par_dist(result$GFO)
sg_gof_res_dist(result$GFO)
sg_gof_res(result$GFO)
sg_gof_res(result$GFO, vs_time = F)
```