search

metrumresearchgroup / simpar

Create Parameters for Simulation with Uncertainty
https://metrumresearchgroup.github.io/simpar/
3 stars 1 forks source link

Use independent rinvchisq draws for diagonal matrix #7

Open kylebaron opened 3 years ago

kylebaron commented 3 years ago

Summary

As a user, I want simpar to simulate OMEGA and SIGMA with all off-diagonals equal to zero when the input OMEGA or SIGMA are diagonal matrices.

Tests

kylebaron commented 3 years ago

Using simblock results in non-zero off-diagonal elements.

Discussed with @curtisKJ and @JimRogersMetrum

kylebaron commented 3 years ago

Current behavior as it stands in metrumrg

library(tidyverse)
library(simpar)
theta <- c(1,2,3,4,5)/10
covar <- diag(0.1, 5)/seq(1,5)
omega <- diag(c(1,2,3,4))
sigma <- diag(c(10,100))

omega
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    2    0    0
#> [3,]    0    0    3    0
#> [4,]    0    0    0    4

set.seed(12345)
a <- simpar(10, theta, covar, omega, sigma)

select(as_tibble(a), contains("OM"))
#> # A tibble: 10 x 10
#>    OM1.1   OM2.1 OM2.2   OM3.1   OM3.2 OM3.3   OM4.1  OM4.2   OM4.3 OM4.4
#>    <dbl>   <dbl> <dbl>   <dbl>   <dbl> <dbl>   <dbl>  <dbl>   <dbl> <dbl>
#>  1 1.50   0.359   2.06  0.0788  0.138   2.38 -0.674  -0.121  0.317   3.56
#>  2 1.42   0.865   2.62 -0.647  -0.888   2.72  0.0879  1.82  -0.0897  7.96
#>  3 1.70  -0.225   1.34  0.740  -0.222   3.94 -0.662  -0.271 -0.942   3.75
#>  4 1.35   0.438   3.77 -0.434  -0.801   2.43  0.725   1.11  -1.37    4.12
#>  5 1.31   0.378   5.84  0.0436  1.18    4.58 -0.361   6.41  -0.810  21.5 
#>  6 1.92   1.25    3.52 -1.31   -1.44    4.44  0.301   0.856  0.113   8.22
#>  7 0.907 -0.0927  1.32 -0.259  -0.163   3.87 -0.221  -0.714  0.476   3.04
#>  8 0.888 -0.486   4.73  0.401  -1.62    2.17  0.270   1.21  -0.465   3.83
#>  9 1.60   0.361   3.11 -0.562   0.0960  3.3   2.48   -2.41  -2.97   16.0 
#> 10 1.12   0.541   5.23  0.201   0.242   3.19  0.231  -1.45  -0.215   4.45

Created on 2021-07-14 by the reprex package (v2.0.0)

kylebaron commented 1 year ago

Updated behavior with source code

Handle diagonal matrices when simulating with metrumrg::simblock

Kyle Baron Tuesday, Nov 29, 2022

library(mrgsolve)
library(metrumrg)

Introduction

Proposed code to update metrumrg::simblock to be able to handle OMEGA and SIGMA matrices with zeros in off-diagonal elements.

The source code

This will be displayed at the bottom of the document

source("simblock-update.R")

Run some tests

OMEGA will be a 3x3, 2x2, 1x1, and 2x2 list of matrices

omega <- list(dmat(1,2,3), bmat(1,0.2,3), matrix(5), cmat(1, 0.94, 3) )
odf <- c(20, 20, 20, 20)

The first matrix is “diagonal”

omega[[1]]
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
is_diagonal(omega[[1]])
## [1] TRUE

The 2nd and 4th matrices are not

omega[[2]]
##      [,1] [,2]
## [1,]  1.0  0.2
## [2,]  0.2  3.0
omega[[4]]
##          [,1]     [,2]
## [1,] 1.000000 1.628128
## [2,] 1.628128 3.000000

simpar already does inverse chi-squared when there is just one row and column

omega[[3]]
##      [,1]
## [1,]    5

mrgsolve-style output

ans1 <- simulate_matrix(10, df = odf, cov = omega, style = "mrgsolve")
ans2 <- simulate_matrix(10, df = odf*1000, cov = omega, style = "mrgsolve")

Note that we have modified the outputs here to better match what mrgsolve is expecting. The output object is a list with length n and each element is a list of matrices that match the pattern of what we passed in as cov. So in this case, ans1 will be a list of length 10 and each position is a list of a 3x3, 2x2, 1x1, and 2x2 matrices. These should have some variability

ans1[[3]]
## [[1]]
##           [,1]     [,2]     [,3]
## [1,] 0.8679706 0.000000 0.000000
## [2,] 0.0000000 2.419742 0.000000
## [3,] 0.0000000 0.000000 3.972554
## 
## [[2]]
##         [,1]    [,2]
## [1,] 1.09199 0.66013
## [2,] 0.66013 2.72499
## 
## [[3]]
##          [,1]
## [1,] 4.261505
## 
## [[4]]
##         [,1]    [,2]
## [1,] 0.76871 1.04220
## [2,] 1.04220 1.94461

Increased degrees of freedom; results should be close to the input matrix

ans2[[3]]
## [[1]]
##          [,1]     [,2]     [,3]
## [1,] 1.009624 0.000000 0.000000
## [2,] 0.000000 1.981528 0.000000
## [3,] 0.000000 0.000000 2.987464
## 
## [[2]]
##         [,1]    [,2]
## [1,] 0.98844 0.20855
## [2,] 0.20855 2.99226
## 
## [[3]]
##          [,1]
## [1,] 5.003767
## 
## [[4]]
##         [,1]    [,2]
## [1,] 1.00828 1.64667
## [2,] 1.64667 3.03274

Check the correlation in the 4th matrix

cov2cor(omega[[4]])
##      [,1] [,2]
## [1,] 1.00 0.94
## [2,] 0.94 1.00
cov2cor(ans2[[3]][[4]])
##           [,1]      [,2]
## [1,] 1.0000000 0.9416693
## [2,] 0.9416693 1.0000000

simpar-style output

We can also get the simpar-style output to verify

ans3 <- simulate_matrix(10, df = odf, cov = omega, style = "simpar")
ans4 <- simulate_matrix(10, df = odf*1000, cov = omega, style = "simpar")

colMeans(ans4)
##     OM1.1     OM2.1     OM2.2     OM3.1     OM3.2     OM3.3     OM4.4     OM5.4 
## 0.9991717 0.0000000 2.0014055 0.0000000 0.0000000 2.9893315 0.9988010 0.1984040 
##     OM5.5     OM6.6     OM7.7     OM8.7     OM8.8 
## 3.0019560 4.9921640 0.9945240 1.6198240 2.9891100
colMeans(ans3)
##     OM1.1     OM2.1     OM2.2     OM3.1     OM3.2     OM3.3     OM4.4     OM5.4 
## 1.1381017 0.0000000 2.2927958 0.0000000 0.0000000 3.5027941 1.0871230 0.3582757 
##     OM5.5     OM6.6     OM7.7     OM8.7     OM8.8 
## 2.9673040 6.2954427 1.2100860 2.0347270 3.8260280

Source code

stopifnot(requireNamespace("metrumrg"))
#' Internal metrumrg function to generate names for outputs
#' We have to include this here because it is only available inside another 
#' metrumrg function
impliedNames <- function(x){
  vars <- sum(x)
  crit <- cumsum(x)-x+1 
  diag <- diag(rep(crit,x))
  good <- outer(colSums(diag),rowSums(diag),FUN='==')
  half <- half(good)
  names(half[half])
}

#' Check if this is a diagonal matrix
#' @param matrix the object to check
#' @param eps tolerance for checking that off diagonals are zero
is_diagonal <- function(matrix, eps = 1e-12) {
  # check if this is a diagonal matrix
  off <- matrix[upper.tri(matrix)]
  all(abs(off)  < eps)
}

#' Implement riwish for a diagonal matrix
#' When we find that all off diagonal elements are 0, we will simulate
#' as a series of independent invchisq simulations
#' 
#' @param n number of simulations to perform; passed to rinvchisq
#' @param df degrees of freedom; passed to rinvchisq
#' @param cov the matrix to simulate
#' 
riwish_diag <- function(n, df, cov) {
  # implement riwish for a diagonal matrix
  # this is a series of calls to metrumrg::rinvchisq
  x <- cov[upper.tri(cov, diag = TRUE)]
  lapply(x, function(xx) rinvchisq(n, df, xx))
}

#' metrumrg::simblock replacement that handles diagonal omega or sigma 
#' 
#' @param n number of simulations; scalar
#' @param df degrees of freedom; scalar
#' @param cov a single matrix to simulate
sblock <- function(n,df,cov) {
  if(df < nrow(cov)) stop('df is less than matrix length')
  if(length(cov)==1) return(rinvchisq(n,df,cov))
  # Handle diagonal omega or sigma
  if(is_diagonal(cov)) {
    res <- do.call(cbind, riwish_diag(n, df, cov))
    return(res)
  }
  # Handle omega or sigma with off diagonal elements
  s <- dim(cov)[1]
  ncols <- s*(s+1)/2
  res <- matrix(nrow=n, ncol=ncols)
  for(i in 1:n)res[i,] <- half(posmat(riwish(s,df-s+1,df*cov)))
  res
}

#' Simulate OMEGA or SIGMA matrices
#' 
#' @param n number of simulations
#' @param df degrees of freedom; must be a list of numeric values with same 
#' length as `cov`
#' @param cov a list of numeric matrices
#' @param prefix added to numbered outputs; only used when `style` is `"simpar"`
simulate_matrix <- function(n, df, cov, prefix = "OM", style = c("simpar", "mrgsolve")) {
  stopifnot(is.list(cov))
  stopifnot(all(sapply(cov, inherits, "matrix")))
  stopifnot(length(df)==length(cov))
  stopifnot(all(sapply(df, is.numeric)))
  stopifnot(length(n)==1)
  stopifnot(is.numeric(n))
  style <- match.arg(style)
  omg <- lapply(seq_along(cov), function(x) list(n = n, df = df[[x]], cov = cov[[x]]))

  # simpar style output
  if(style=="simpar") {
    omg <- do.call(cbind,lapply(omg, function(x) do.call(sblock, x)))
    labels <- paste0(prefix, impliedNames(sapply(cov, ord)))
    dimnames(omg) <- list(NULL,labels)
    return(omg)
  }

  # mrgsolve style output
  if(style=="mrgsolve") {
    requireNamespace("mrgsolve")
    omg <- lapply(omg, function(x) do.call(sblock,x))
    omg <- lapply(omg, function(x) mrgsolve::as_bmat(as.data.frame(x)))
    omg <- lapply(seq(n), function(nn) {
      lapply(omg, function(x) {
        x[[nn]]
      })
    })
    return(omg)
  }
  stop("something bad happened")
}
kylebaron commented 1 year ago

Direction