-
Notifications
You must be signed in to change notification settings - Fork 1
/
NMsim-ResidVar.Rmd
99 lines (83 loc) · 3.69 KB
/
NMsim-ResidVar.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
---
title: "NMsim - Simulation of residual variability"
output:
rmarkdown::html_vignette:
toc: true
code_folding: show
Suggests: markdown
VignetteBuilder: knitr
vignette: >
%\VignetteIndexEntry{Residual}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
header-includes:
- \usepackage{ae}
---
```{r,include = FALSE}
##knitr::opts_chunk$set(dev = "cairo_pdf")
knitr::opts_chunk$set(
collapse = TRUE
,comment = "#>"
,fig.width=7
,cache=FALSE
)
## this changes data.table syntax. I think we can do without.
## knitr::opts_chunk$set(tidy.opts=list(width.cutoff=60), tidy=TRUE)
```
```{r,setup,include=F}
## library(devtools)
## unloadNamespace("NMsim")
## unloadNamespace("NMdata")
## load_all("~/wdirs/NMdata")
## load_all()
library(NMsim)
library(data.table)
library(NMdata)
## library(dplyr)
library(tibble)
library(ggplot2)
library(patchwork)
library(tracee)
## library(tidyr)
library(fst)
library(knitr)
## NMdataConf(path.nonmem="/opt/NONMEM/nm75/run/nmfe75")
## NMdataConf(path.nonmem="/opt/nonmem/nm751/run/nmfe75")
## NMdataConf(dir.psn=NULL)
theme_set(theme_bw())
this.script <- "NMsim-simulate.Rmd"
writeOutput <- TRUE
file.project <- function(...)file.path(system.file("examples",package="NMsim"),...)
## file.project <- function(...)file.path("../inst/examples",...)
## file.project <- function(...)file.path("~/wdirs/NMsim/inst/examples",...)
found.files <- list.files(file.project("nonmem/NMsim"),pattern="noname\\.(lst|xml|ext|cov|cor|coi|phi|msf|msfi|msfo|tab)",full.names=TRUE)
unlink(found.files)
run.simuls <- FALSE
NMdataConf(as.fun="data.table")
```
## Add residual variability
The best way to simulate with residual variability is to include the it in the estimation control stream as described in this vignette. `NMsim` currently does not provide any automated way to add simulation of
residual variability with Nonmem. It does provide a method to
simulate residual variability in R, based on the Nonmem parameter
estimates. This should only be used in case one has an existing Nonmem without residual variability simulated, and it is not feasible to modify the model control stream for some reason. The function is called `addResVar()` and supports additive,
proportional, and combined (additive and proportional) error
models. It can also add the residual error on log scale (exponential
error model).
`addResVar` supports both estimation using `$SIGMA` and `$THETA` (in Nonmem). The user has to specify which of the two methods were used in the Nonmem model using the `par.type` argument. The other thing that must be specified is the parameter numbers for the standard deviations or variances. The model simulated in this vignette has this combined error model estimated using the `$SIGMA` matrix:
```
Y=F+F*ERR(1)+ERR(2)
```
We now specify for `addResVar` to find the variance for the proportional component in `$SIGMA[1,1]` and the one for the additive component in `$SIGMA[2,2]`. In this case where `SIGMA` is used, the off-diagonal (covariance) elements of the `$SIGMA` matrix are also used for the simulation.
```{r,sim-simplest,eval=FALSE}
file.mod <- file.project("nonmem/xgxr021.mod")
simres <- NMsim(file.mod=file.mod,
data=dat.sim)
```
```{r,eval=FALSE}
simres.with.resvar <- addResVar(simres,path.ext=fnExtension(file.mod,"ext"),par.type="SIGMA",prop=1,add=2)
```
If `par.type="THETA"` the default assumption is that the thetas
represent standard deviation (in contrast to when using
`par.type="SIGMA"`). This can be modified using the `scale.par`
argument. There are arguments to avoid negative observations and
several other features. But again, this should be the last resort.