[WIP] Refactor simblock etc

[shenc-metrum]

Summary

Example

library(simpar)
library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

Create some simpar() inputs

theta <- c(0.125, 0.333, 1.51)
covar <- diag(c(0.01, 0.02, 0.8))
omega <- list(bmat(0.1, 0.05, 0.2), dmat(0.1, 0.2, 0.3))
sigma <- matrix(0.1)

Check omega

omega
#> [[1]]
#>      [,1] [,2]
#> [1,] 0.10 0.05
#> [2,] 0.05 0.20
#> 
#> [[2]]
#>      [,1] [,2] [,3]
#> [1,]  0.1  0.0  0.0
#> [2,]  0.0  0.2  0.0
#> [3,]  0.0  0.0  0.3

Traditional format

By default, simpar() still gives you the same output it always has … as a data.frame

set.seed(1357)
x <- simpar(
  nsim = 2,
  theta = theta,
  covar = covar,
  omega = omega,
  sigma = sigma,
  odf = list(50, 50),
  sdf = 1000
)

x
#>       TH.1   TH.2  TH.3  OM1.1   OM2.1 OM2.2   OM3.3    OM4.3  OM4.4     OM5.3
#> 1 -0.05954 0.1942 1.837 0.1086 0.07823 0.273 0.10800 -0.02187 0.1957 -0.009684
#> 2  0.19850 0.3855 3.090 0.1186 0.07770 0.272 0.08704  0.07370 0.3623 -0.000002
#>      OM5.4  OM5.5   SG1.1
#> 1 -0.03524 0.2679 0.10320
#> 2  0.09462 0.2445 0.09873

Updated format for mrgsolve

This output is a list of lists, with one draw from the THETA, SIGMA, OMEGA distributions in each slot of the list.

We also show here new, optional behavior where simpar() will fix off-diagonal elements in simulated output when the input $OMEGA matrix has zeros on the off-diagonals.

set.seed(6428)
x <- simpar(
  nsim = 2,
  theta = theta,
  covar = covar,
  omega = omega,
  sigma = sigma,
  odf = list(50, 50),
  sdf = 1000,
  mrgsolve_style = TRUE,
  omega_diag = TRUE
)

x[[1]]
#> $param
#>   THETA1 THETA2 THETA3
#> 1 0.1182 0.5727 0.8184
#> 
#> $omega
#> $omega[[1]]
#>         [,1]    [,2]
#> [1,] 0.11520 0.04663
#> [2,] 0.04663 0.16250
#> 
#> $omega[[2]]
#>         [,1]   [,2]   [,3]
#> [1,] 0.09439 0.0000 0.0000
#> [2,] 0.00000 0.2008 0.0000
#> [3,] 0.00000 0.0000 0.2243
#> 
#> 
#> $sigma
#> $sigma[[1]]
#>        [,1]
#> [1,] 0.1011

The idea is to be able to send this right into mrgsolve

mod <- modlib("1005", capture = "THETA1") %>% ev(amt = 100)
#> Loading required namespace: xml2
#> Building 1005 ... done.

set.seed(0998)
x <- simpar(
  nsim = 2,
  theta = theta,
  covar = covar,
  omega = as.matrix(omat(mod)),
  sigma = as.matrix(smat(mod)),
  odf = 50,
  sdf = 1000,
  mrgsolve_style = TRUE,
  omega_diag = TRUE
)

out <- lapply(x, function(draw) {
  mod <- update(mod, data = draw)
  mrgsim(mod, output = "df")
})

^{Created on 2023-03-10 with reprex v2.0.2}