We use the epilepsy example, explained in more detail in the epil notebook, as example. Start by loading and preparing the data, as well as compling the TMB template:

data(Epil, package = "INLA")

center <- function(x) (x - mean(x))

Epil <- Epil %>%
  mutate(CTrt    = center(Trt),
         ClBase4 = center(log(Base/4)),
         CV4     = center(V4),
         ClAge   = center(log(Age)),
         CBT     = center(Trt * log(Base/4)))

N <- 59
J <- 4
K <- 6
X <- model.matrix(formula(~ 1 + CTrt + ClBase4 + CV4 + ClAge + CBT), data = Epil)
y <- Epil$y

make_epsilon_matrix <- function(N, J) {
  t(outer(1:N, 1:(N * J), function(r, c) as.numeric((J*(r - 1) < c) & (c <= J*r))))
}

E <- make_epsilon_matrix(N = 3, J = 2)

dat <- list(N = N, J = J, K = K, X = X, y = y, E = make_epsilon_matrix(N, J))

compile("epil.cpp")
## [1] 0
dyn.load(dynlib("epil"))

param <- list(
  beta = rep(0, K),
  epsilon = rep(0, N),
  nu = rep(0, N * J),
  l_tau_epsilon = 0,
  l_tau_nu = 0
)

obj <- MakeADFun(
  data = dat,
  parameters = param,
  random = c("epsilon", "nu"),
  DLL = "epil",
  silent = TRUE
)

aghq::marginal_laplace_tmb

We start by calling this function, to show where we’re trying to end up:

fit <- aghq::marginal_laplace_tmb(
  obj,
  k = 2,
  startingvalue = c(param$beta, param$l_tau_epsilon, param$l_tau_nu)
)
summary(fit)
## There are 295 random effects, but max_print = 30, so not computing their summary information.
## Set max_print higher than 295 if you would like to summarize the random effects.
## AGHQ on a 8 dimensional posterior with  2 2 2 2 2 2 2 2 quadrature points
## 
## The posterior mode is: 1.578628 -0.9485647 0.8803436 -0.1030889 0.4901293 0.3492073 1.555814 2.059541 
## 
## The log of the normalizing constant/marginal likelihood is: -679.2151 
## 
## The covariance matrix used for the quadrature is...
##               [,1]          [,2]          [,3]          [,4]          [,5]          [,6]          [,7]          [,8]
## [1,]  5.392878e-03  1.505446e-03 -5.983583e-04  2.944984e-05 -3.928243e-04 -5.832637e-04  1.885191e-03  0.0021018979
## [2,]  1.505446e-03  1.564623e-01  3.044970e-02  2.695215e-05 -2.584180e-02 -7.372289e-02  1.750390e-03  0.0009247506
## [3,] -5.983583e-04  3.044970e-02  1.671006e-02 -2.758497e-05 -1.656679e-03 -1.692942e-02  3.106606e-04 -0.0003139452
## [4,]  2.944984e-05  2.695215e-05 -2.758497e-05  7.426245e-03  4.152776e-05  1.798204e-05  5.293036e-05 -0.0011614479
## [5,] -3.928243e-04 -2.584180e-02 -1.656679e-03  4.152776e-05  1.171811e-01  1.738378e-02  3.604741e-03 -0.0008080433
## [6,] -5.832637e-04 -7.372289e-02 -1.692942e-02  1.798204e-05  1.738378e-02  4.015874e-02 -3.305936e-05 -0.0005011597
## [7,]  1.885191e-03  1.750390e-03  3.106606e-04  5.293035e-05  3.604741e-03 -3.305935e-05  7.475428e-02 -0.0087276816
## [8,]  2.101898e-03  9.247507e-04 -3.139452e-04 -1.161448e-03 -8.080433e-04 -5.011598e-04 -8.727682e-03  0.0569622640
## 
## Here are some moments and quantiles for the transformed parameter:
## 
##                     mean         sd        2.5%      median      97.5%
## beta           1.5749050 0.07334178  1.43851847  1.57583464  1.7255004
## beta1         -0.9532869 0.39549381 -1.70007626 -0.93193488 -0.1574583
## beta2          0.8812315 0.12936765  0.63526602  0.88991894  1.1388784
## beta3         -0.1036569 0.08617377 -0.26681406 -0.09877583  0.0692624
## beta4          0.4881112 0.34226987 -0.16023910  0.50863294  1.1747637
## beta5          0.3505870 0.20057825 -0.03072352  0.36485385  0.7500002
## l_tau_epsilon  1.4392481 0.24762163  1.02103120  1.33193847  2.0249114
## l_tau_nu       2.0437606 0.23713622  1.60322927  2.04090614  2.5368762

And now let’s go through line by line. These are all of the inputs to the function:

ff <- obj
k <- 2
startingvalue <- c(param$beta, param$l_tau_epsilon, param$l_tau_nu)
transformation <- aghq::default_transformation()
optresults <- NULL
basegrid <- NULL
control <- aghq::default_control_tmb()
#' My guess is that this function just makes sure that the control input is valid
validate_control(control, type = "tmb")
## [1] TRUE
#' And similarly, this probably makes sure the transformation input is valid
validate_transformation(transformation)
## [1] TRUE
transformation <- make_transformation(transformation)
transformation$totheta
## function (x) 
## x
## <bytecode: 0x7fac94108f08>
## <environment: namespace:base>
transformation$fromtheta
## function (x) 
## x
## <bytecode: 0x7fac94108f08>
## <environment: namespace:base>
transformation$jacobian
## function (theta) 
## {
##     out <- numeric(length(theta))
##     for (i in 1:length(theta)) {
##         out[i] <- det(abs(numDeriv::jacobian(totheta, fromtheta(theta[i]))))
##     }
##     out
## }
## <bytecode: 0x7fac56319698>
## <environment: 0x7fac46431af8>
#' Get names for the theta parameters from the TMB function template
#' Note that make.unique will make the names unique by apending sequence numbers to duplicates
thetanames <- NULL
if (exists('par', ff)) thetanames <- make.unique(names(ff$par), sep = "")
thetanames
## [1] "beta"          "beta1"         "beta2"         "beta3"         "beta4"         "beta5"         "l_tau_epsilon" "l_tau_nu"
#' Hessian
if (control$numhessian) {
  ff$he <- function(theta) numDeriv::jacobian(ff$gr, theta, method = "Richardson")
}
#' Now we're going to do the AGHQ
quad <- aghq(
  ff = ff,
  k = k,
  transformation = transformation,
  startingvalue = startingvalue,
  optresults = optresults,
  basegrid = basegrid,
  control = control
)
#' There is an option whereby quad will just return the normalising constant, and if so then tmb_marginal_laplace can also be used
#' just to return the normalising constant too. Note on this, why use tmb_marginal_laplace in this situation rather than aghq? Since
#' as far as I understand we're not doing the marginal laplace approximation if control$onlynormconst is TRUE
if (control$onlynormconst) return(quad)

distinctthetas <- quad$normalized_posterior$nodesandweights[, grep('theta', colnames(quad$normalized_posterior$nodesandweights))]
distinctthetas

Note that there are 256 “distinct thetas” corresponding to \(k\), here 2, raised to the power 8:

dim(distinctthetas)
## [1] 256   8
#' I think this will just be fixing a case where there is only one theta, and somehow a data.frame isn't created
if (!is.data.frame(distinctthetas)) distinctthetas <- data.frame(theta1 = distinctthetas)

modesandhessians <- distinctthetas

#' Replacing the theta, theta1, theta2, ... names with the names from TMB
if (is.null(thetanames)) {
  thetanames <- colnames(distinctthetas)
} else {
  colnames(modesandhessians) <- thetanames
  colnames(quad$normalized_posterior$nodesandweights)[grep('theta', colnames(quad$normalized_posterior$nodesandweights))] <- thetanames
}

modesandhessians$mode <- vector(mode = "list", length = nrow(distinctthetas))
modesandhessians$H <- vector(mode = "list", length = nrow(distinctthetas))

#' There is a loop coming up, which we will step through for i = 1, then let a loop do for i = 2, ...
i <- 1

#' Get the theta (currently modesandhessians, despite being called modesandhessians, only contains the theta, not any mode or hessian)
theta <- as.numeric(modesandhessians[i, thetanames])
theta
## [1]  1.5051918 -1.3640864  0.7078521 -0.1894402  0.1930595  0.4830202  1.2314444  1.8491570
#' Set up the mode and hessian of the random effects. This happens when you run the TMB objective with a particular theta
ff$fn(theta)
## [1] 675.2611
## attr(,"logarithm")
## [1] TRUE
#' Now pull the mode and hessian. Have to be careful about scoping

#' So this is the mode, for both random and fixed effects
mm <- ff$env$last.par

#' These are the indices of the random effects
ff$env$random
##   [1]   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
##  [43]  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
##  [85]  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
## [127] 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
## [169] 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## [211] 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
## [253] 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
## [295] 301
#' We can confirm that for the fixed effects, the "mode" is just the input values theta
#' (To get the indicies for the fixed effects, we just take the complement of the random effects)
mm[-ff$env$random]
##          beta          beta          beta          beta          beta          beta l_tau_epsilon      l_tau_nu 
##     1.5051918    -1.3640864     0.7078521    -0.1894402     0.1930595     0.4830202     1.2314444     1.8491570
theta
## [1]  1.5051918 -1.3640864  0.7078521 -0.1894402  0.1930595  0.4830202  1.2314444  1.8491570
#' Set the mode using the relevant indicies of mm
#' "mode" is a list column
modesandhessians[i, "mode"] <- list(list(mm[ff$env$random]))
modesandhessians[i, ]
#' This should be the Hessian evaluated at the mode
#' The argument random = TRUE set it so that...
H <- ff$env$spHess(mm, random = TRUE)

#' Alex has the following note on duplicating H here:
#' "Somehow, TMB puts all evaluations of spHess in the same memory location"
#' Don't totally understand this
H <- rlang::duplicate(H)

modesandhessians[i, "H"] <- list(list(H))

for (i in 2:nrow(distinctthetas)) {
  theta <- as.numeric(modesandhessians[i,thetanames])
  ff$fn(theta)
  mm <- ff$env$last.par
  modesandhessians[i,'mode'] <- list(list(mm[ff$env$random]))
  H <- ff$env$spHess(mm,random = TRUE)
  H <- rlang::duplicate(H)
  modesandhessians[i,'H'] <- list(list(H))
}

quad$modesandhessians <- modesandhessians

class(quad) <- c("marginallaplace", "aghq")

aghq::sample_marginal.marginallaplace

What we’d like to do is to sample from \[ \begin{align} \tilde p(x_i \, | \, y) &= \sum_z \tilde p_\text{G}(x_i \, | \, \theta(z), y) \tilde p_\text{LA}(\theta(z) \, | \, y) w_k(z), \\ &= \sum_z \tilde p_\text{G}(x_i \, | \, \theta(z), y) \lambda(z) \end{align} \] where \(\lambda(z) \in [0, 1]\) are weights such that \(\sum_z \lambda(z) = 1\). We can do this by

  1. Sample \(z_j \sim \text{Multinomial}(1, \lambda(z_1), \ldots, \lambda(z_J))\)
  2. Given \(z_j\) sample from \(\tilde p_\text{G}(x \, | \, \theta(z_j), y)\)
  3. Take the \(i\)th entry of \(\tilde p_\text{G}(x \, | \, \theta(z_j), y)\) as the sample from \(\tilde p_\text{G}(x_i \, | \, \theta(z_j), y)\)

Now that we have quad, it should be compatible with sample_marginal:

#' The number of samples to get
M <- 10

transformation <- quad$transformation
interpolation <- "auto"

#' Is this getting mc.cores, and if there isn't a value just setting it to 1 (integer)
numcores <- getOption("mc.cores", 1L)

#' Length of the random effects
d <- dim(quad$modesandhessians$H[[1]])[1]

simlist <- quad$modesandhessians

#' Avoid use of parallel computing on Windows
if (.Platform$OS.type == 'windows') numcores <- 1

The next step is to compute the Cholesky decomposition L of the Hessian H such that \[ H = LL^\top \] The reason we want the Cholesky is that we will be using it to same from the joint Gasusian approximation. The function chol from base will compute it for us. For example, for the first element:

h <- simlist$H[[1]]
h <- Matrix::forceSymmetric(h)
L <- chol(h, perm = FALSE)
L
## 295 x 295 sparse Matrix of class "dtCMatrix"
##                                                                                                                                                                                         
## [1,] 4.175814 .        .        .        .        .        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.9969627 0.8183811
## [2,] .        4.173932 .        .        .        .        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .         .        
## [3,] .        .        3.730687 .        .        .        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .         .        
## [4,] .        .        .        4.017652 .        .        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .         .        
## [5,] .        .        .        .        7.619454 .        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .         .        
## [6,] .        .        .        .        .        5.090901 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .         .        
##                                                                                                                                                                                    
## [1,] 0.8183811 0.721608 .         .         .         .         .         .         .         .         .         .         .         .         .        .        .        .       
## [2,] .         .        0.8178211 0.9963606 0.8178211 0.7210789 .         .         .         .         .         .         .         .         .        .        .        .       
## [3,] .         .        .         .         .         .         0.6608374 0.8230477 0.5230443 0.8053812 .         .         .         .         .        .        .        .       
## [4,] .         .        .         .         .         .         .         .         .         .         0.8769562 0.8769562 0.6366786 0.7742805 .        .        .        .       
## [5,] .         .        .         .         .         .         .         .         .         .         .         .         .         .         1.228496 2.190027 1.388735 2.362534
## [6,] .         .        .         .         .         .         .         .         .         .         .         .         .         .         .        .        .        .       
##                     
## [1,] .        ......
## [2,] .        ......
## [3,] .        ......
## [4,] .        ......
## [5,] .        ......
## [6,] 1.076779 ......
## 
##  ..............................
##  ........suppressing 215 columns and 283 rows in show(); maybe adjust 'options(max.print= *, width = *)'
##  ..............................
##                                                                                                                                                                              
## [290,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
## [291,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
## [292,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
## [293,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
## [294,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
## [295,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......

Doing this for the list using lapply:

simlist$L <- lapply(simlist$H,function(h) chol(Matrix::forceSymmetric(h), perm = FALSE))

There is another version, using the arguments perm and LDL (which I don’t understand) that may be quicker but ultimately will achieve the same thing:

# simlist$L <- lapply(simlist$H, function(h) Matrix::Cholesky(as(Matrix::forceSymmetric(h), "sparseMatrix"), perm = TRUE, LDL = FALSE))

Now let’s compute the \(\lambda(z)\) weights.

simlist$lambda <- exp(quad$normalized_posterior$nodesandweights$logpost_normalized) * quad$normalized_posterior$nodesandweights$weights
head(simlist$lambda)
## [1] 0.006359709 0.006024268 0.006375453 0.005857861 0.006285794 0.006311227

We can check that they do indeed sum to one:

sum(simlist$lambda)
## [1] 1

Now let’s sample from the multinomial:

j <- apply(stats::rmultinom(M, 1, simlist$lambda), 2, function(x) which(x == 1))
length(j) #' Should be M
## [1] 10
j
##  [1]  81   9 150 114 142 241   2 192 178 169
#' Values are number of samples to draw for each j
#' In this case they're likely to be all 1, since we're not doing many samples
tt <- table(j)

Now we construct a big Gaussian mixture matrix. The first step is to create a matrix with d columns, one for each random effect, and M rows, one for each sample, where each value is a standard normal \(\mathcal{N}(0, 1)\):

A <- matrix(stats::rnorm(M * d), nrow = M)

#' Now make a list where each element is a row of A, where the index is determined by j
A_split <- split(A, j)

#' Big Gaussian mixture matrix

#' For each element of the list, create a matrix with d rows (and 1 column? always?)
Z <- lapply(A_split, matrix , nrow = d)
str(Z)
## List of 10
##  $ 2  : num [1:295, 1] 0.0714 1.2809 0.1427 0.5559 -0.8705 ...
##  $ 9  : num [1:295, 1] 1.643 -0.481 -0.987 -0.91 -1.241 ...
##  $ 81 : num [1:295, 1] 0.00409 0.72422 -0.4331 -0.05492 -0.41522 ...
##  $ 114: num [1:295, 1] 2.1462 -0.7045 0.0772 0.4656 0.6458 ...
##  $ 142: num [1:295, 1] 0.6 -2.316 -0.201 2.28 0.238 ...
##  $ 150: num [1:295, 1] 1.228 -0.336 -0.711 0.137 -2.491 ...
##  $ 169: num [1:295, 1] -0.2433 -0.0836 2.1579 -2.1105 -1.0353 ...
##  $ 178: num [1:295, 1] -0.423 0.224 -0.73 -1.061 1.236 ...
##  $ 192: num [1:295, 1] -0.0852 1.3955 -0.0511 0.3434 1.0585 ...
##  $ 241: num [1:295, 1] -0.918 0.442 -1.587 -0.861 0.319 ...

In the code there is now a big call to mapply, which we are going to break down step-by-step. First set the inital values that function(.x, .y) will iterate over:

.x <- Z[[1]]
str(.x)
##  num [1:295, 1] 0.0714 1.2809 0.1427 0.5559 -0.8705 ...
.y <- names(Z)[[1]]
str(.y)
##  chr "2"

Now solve the equation L %*% .z = .x for .z

.z <- solve(simlist$L[[as.numeric(.y)]], .x)
str(.z)
## Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
##   ..@ Dim     : int [1:2] 295 1
##   ..@ Dimnames:List of 2
##   .. ..$ : NULL
##   .. ..$ : NULL
##   ..@ x       : num [1:295] 0.00481 0.34322 0.19887 0.00238 -0.24044 ...
##   ..@ factors : list()

Now add the mode (where the additional manipulation is so that if there are multiple samples from the same j then there will be the corresonding number of columns):

.z <- .z + do.call(cbind, rep(list(simlist$mode[[as.numeric(.y)]]), ncol(.x)))
str(.z)
## Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
##   ..@ Dim     : int [1:2] 295 1
##   ..@ Dimnames:List of 2
##   .. ..$ : chr [1:295] "epsilon" "epsilon" "epsilon" "epsilon" ...
##   .. ..$ : NULL
##   ..@ x       : num [1:295] -0.1073 0.2354 0.2356 -0.0245 -0.2359 ...
##   ..@ factors : list()

Do this using a mapply call:

samps <- mapply(
  function(.x, .y) as.numeric(solve(simlist$L[[as.numeric(.y)]], .x)) + do.call(cbind, rep(list(simlist$mode[[as.numeric(.y)]]), ncol(.x))),
  Z,
  names(Z),
  SIMPLIFY = FALSE
)

#' Order them properly
#' My guess is that you'd want them to be ordered according to how j was drawn, rather
#' than being grouped up by the value of lambda
ord <- numeric(length(j))
cumtab <- cumsum(c(0, tt))
cumtab <- cumtab[-length(cumtab)]
cnt <- numeric(length(unique(j)))
names(cnt) <- sort(unique(j))

for (i in 1:length(j)) {
  wc <- which(names(cnt) == j[i])
  cnt[wc] <- cnt[wc] + 1
  ord[i] <- cumtab[wc] + cnt[wc]
}

samps <- Reduce(cbind, samps)
samps <- samps[, ord]
samps
##                  [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]          [,8]          [,9]         [,10]
## epsilon -0.1713546458  0.389022874  0.152521897  0.470493882 -0.012219115 -0.190824045 -0.107323691  0.0625516298 -0.2096090324  2.531406e-03
## epsilon  0.1106289542 -0.134366457 -0.218479076 -0.246456540 -0.710196069  0.153568727  0.235403601  0.3643596032  0.1204384257  2.692746e-02
## epsilon -0.1415398653 -0.207014556  0.125563770  0.234674444  0.107799869 -0.144041573  0.235568434  0.3861041574  0.0598286223  8.605439e-01
## epsilon -0.0337958995 -0.187575488  0.006493643  0.188932119  0.607924446 -0.121668262 -0.024515312  0.2876272866 -0.2983284799 -2.256549e-01
## epsilon  0.2427837542  0.026336466 -0.401174563  0.054495035 -0.274370995  0.182292076 -0.235917928  0.1830314893  0.4560876982  1.149050e-01
## epsilon -0.3636780044 -0.492459636 -0.532600451  0.008650655 -0.057790355 -0.265716834  0.026031497  0.0427448359 -0.6981300053 -9.876891e-02
## epsilon -0.1192593624 -0.152029974 -0.888738469 -0.148845153 -0.537155021 -0.583634995 -0.174099305 -0.0760930647 -0.3448809335 -4.626149e-01
## epsilon  0.5141940902  0.618918547  0.044459612  0.330201237  0.600198671  0.729959896  0.362642867 -0.0272280724  0.6996757232  5.146388e-01
## epsilon -0.2992740943  0.125322850 -0.536013925 -0.244755115 -0.338146379 -0.367205071 -0.302650747 -0.3688465408  0.4079418657  2.000648e-01
## epsilon  0.6100716485  0.877791542  0.437093754  0.675182963  0.464290993  0.967151690  0.250326604  1.1803860794  0.9548306214  5.175403e-01
## epsilon  0.3170539621 -0.204037524 -0.159280044 -0.087551836  0.158874238  0.475545947 -0.142104221 -0.1241084525 -0.1714607058  2.675668e-01
## epsilon -0.1345684886  0.143542540 -0.105146467  0.106638817  0.054786285  0.256599323 -0.050338761 -0.0904872974 -0.2732322190  2.089260e-01
## epsilon  0.4163677825 -0.274166873 -0.419639450  0.064901928 -0.415158987  0.593051028  0.083589478  0.2249247077  0.4309855293  1.451465e-01
## epsilon -0.2794955174 -0.626992916 -0.426185614 -0.099286934  0.077849553 -0.099794574  0.160543994 -0.0431068500 -0.2481180613  7.681694e-02
## epsilon  0.1289643266 -0.079147711 -0.153194357  0.090546635 -0.490694213 -0.182010913 -0.071902232 -0.4174213243  0.5058955968  3.430933e-01
## epsilon -0.4779585908 -1.291315020 -0.741434172 -0.230251285 -1.055966455 -0.683036854 -0.839589374 -0.8641537154 -0.6335677020 -8.035414e-01
## epsilon -0.4361246393 -0.905275335 -0.880134554 -0.792305167 -1.066979038 -0.315350411 -0.795679650 -0.1790474236 -0.6389997986 -8.464472e-01
## epsilon  0.6719475384  0.748661208  0.357182301  0.413891425  0.073438343  0.734644855  0.159858901 -0.2369279987  0.3303313657  8.345306e-01
## epsilon -0.0224822378 -0.208816367 -0.807054326 -0.081635662 -0.059696979  0.134605215  0.291923358  0.0005960436 -0.6879147507  1.283643e-01
## epsilon  0.2104559039  0.033724972  0.040330114  0.272373939 -0.558507965 -0.204242815 -0.198663063 -0.2936892108  0.3053510247  1.667583e-02
## epsilon -0.0425284987 -0.262586222  0.274338109  0.021511396  0.116504034 -0.032198945 -0.479299515  0.1322380123 -0.0149582834 -7.597557e-03
## epsilon -0.0827547440  0.762701540  0.402077975  0.412628275 -0.321204332  0.503764004  0.498346405  0.7127311677 -0.1468700075  5.566111e-01
## epsilon -0.2297205319 -0.687466379 -0.341566581 -0.367233813  0.201369823  0.142197476 -0.599743472 -0.4377792483 -0.3627712317 -2.364402e-01
## epsilon  0.1767682574  0.193941469  0.206484100  0.123718513 -0.399749897  0.207172124 -0.065773158  0.0245163300  0.3229227805 -2.891482e-02
## epsilon  0.6433692374  0.945250818  0.733369145  0.466764384  0.711070591  0.787127136  0.599824093  0.5884106719  1.0456510011  1.144872e+00
## epsilon -0.6705814539 -0.440947870 -0.733169455 -0.689248080 -0.340755468 -0.498188228 -0.412828740 -0.9631872318 -0.5129935637 -8.726338e-01
## epsilon  0.3526064381 -0.827975194  0.103350033 -0.523501464  0.007336369  0.169079855 -0.196548326  0.1923590177  0.4417033226 -6.638183e-02
## epsilon  0.4527610254  0.500382151  0.002825913  0.393752780 -0.173584383  0.405490986  0.497901599  0.3600116111  0.3166010298  2.410882e-01
## epsilon -0.2731182940 -0.667748441 -0.406103038 -0.496286248 -0.235052219 -0.128292562 -0.556089546  0.0126651245 -0.3433591408 -5.242921e-01
## epsilon -0.2046312413  0.011935573 -0.017829558 -0.057968326  0.122000009 -0.002577661  0.012825782 -0.6992717620 -0.3687621702 -5.289882e-01
## epsilon -0.5063165873 -0.420796146 -0.105628764 -0.298233536 -0.563144435  0.459860959 -0.471047915 -0.0018500227 -0.6132876618  1.223206e-01
## epsilon  0.9589053304  0.633527893  0.612017323  0.162219890  0.926190363  0.793774725  0.801448920  0.3359575090  0.6496937760  1.530952e-01
## epsilon  0.4077795164  0.972430371  0.644180908  0.290946239 -0.148907762  0.545022616  0.266463466  0.4716389445  0.7298179093  3.101816e-01
## epsilon -0.1386571754 -0.260473389 -0.204674848 -0.092625493 -0.609398864 -0.527691490 -0.461168254 -0.1431688083 -0.0280018160 -6.440878e-02
## epsilon  1.0481600204  1.161569840  1.038893297  0.713273202  0.868860488  0.888736211  0.618722449  0.6592737572  0.9627734597  1.089023e+00
## epsilon  0.5068166807  0.452509140  0.771661536 -0.127475594  0.665975868  0.592740822  0.650495854  0.1377869335  0.6197787545  2.930961e-01
## epsilon  0.7570210118  0.076020204  0.126556693  0.132664612  0.376242568  0.758345335  0.298496867  0.3638329015  0.7894181977  3.675516e-01
## epsilon -0.4056937648 -0.373546103 -0.734241273 -0.091244277 -0.493388115 -0.675965788 -0.996822155 -0.5729928584 -0.6275519635 -9.410832e-01
## epsilon  0.2435515548 -0.047445660  0.279883546  0.293081348 -0.128256121 -0.035194933 -0.387767835 -0.2321350705  0.2509620370 -3.297987e-02
## epsilon -0.0603480078 -0.284049524 -0.744553883 -0.506629340  0.144201868 -0.047230574  0.457434060 -0.2593714222  0.2092039144  3.049175e-01
## epsilon -0.2766836911 -0.463143523 -0.767873092 -0.566848524 -0.824718450 -0.259266437 -0.632075487 -0.2137496728 -0.6334997654 -3.166307e-01
## epsilon  0.2833580938  1.121579013 -0.251831506  0.226484686  0.552783482  0.175873702 -0.328847921 -0.0013157359 -0.0982114316  7.156709e-01
## epsilon  0.5338103234  0.196756615  0.202715782  0.401468163  0.550180158  0.625155310  0.808310952  0.2693192151 -0.0882534628  5.349954e-01
## epsilon  0.4688574613 -0.046394560  0.023356376  0.163901839 -0.085270347  0.212976991 -0.200486973 -0.1487863116  0.0181865245  3.183495e-01
## epsilon  0.0571076726  0.786809063  0.009330139  0.059058403  0.413006624 -0.306588366  0.341416881 -0.2378411200  0.0885915502  3.158964e-01
## epsilon  0.2536015687  0.581621279  0.595248063  0.155224466  0.260895684  0.671388415  0.361848349 -0.1121106252  0.5100441796  5.229212e-01
## epsilon -0.0501982516  0.066493773  0.323828399  0.176781177  0.007695045  0.065976343 -0.122497975  0.0452288395 -0.1867048284 -3.769494e-02
## epsilon  0.4244526355 -0.311505006  0.132584368 -0.053064263 -0.098772716 -0.707523671 -0.079081070 -0.7355968744 -0.4028237614 -5.359148e-01
## epsilon  0.9678283077  0.685904520  0.281479259  0.332600321  0.811796760  0.698081762  0.360160759  0.6518412561  0.6429873989  2.253820e-01
## epsilon  0.0791330571 -0.785975903 -0.260103904 -0.226646810 -0.458823224  0.129800503 -0.388062173 -0.3651886440 -0.9176275707 -1.460296e-01
## epsilon  0.0213157106  0.005599186 -0.256911138 -0.148092269 -0.297839408 -0.167647015 -0.149687357  0.0911167918 -0.4835436682 -6.389711e-02
## epsilon -0.9962765184 -0.210664076 -0.686028973 -0.833694755 -0.988848264 -0.845192461 -0.285118522 -0.8563795776 -0.5195876616 -4.583191e-01
## epsilon  0.7324812916  0.579031128  0.199707637  0.280437275  0.392305852  0.377774816  0.296977561  0.3824091157  0.7021166558  3.613398e-01
## epsilon -0.2411070429  0.023810207 -0.312451386 -0.641808343 -0.608549535 -0.407695964  0.717260024 -0.5335108209 -0.5332883297  1.126761e-01
## epsilon  0.5078836281  0.358150501  0.419460116 -0.132270827  0.417917283  0.225899453  0.669846075  0.3553300463  0.0510229073  1.082207e-01
## epsilon  0.7316830746  0.806729666  1.141717744  0.683805028  1.383744387  1.121409970  1.377280223  0.8848713292  0.7843548550  9.448439e-01
## epsilon -0.6184744425 -0.469815751  0.061357163 -0.637677147 -0.795017004 -0.786767261  0.021097478 -0.3357223170 -0.5798175699 -6.205624e-01
## epsilon -0.5644790093 -1.237913445 -0.793970401 -1.009678447 -1.621150688 -0.565016160 -0.803669376 -1.2315783408 -1.6631103814 -1.093177e+00
## epsilon  0.4225471703  0.428528566 -0.360798361  0.169838180  0.070310749  0.224948953 -0.016641142 -0.0045639656  0.0242673403  4.047505e-01
## nu       0.1804105015  0.080234543  0.332988713  0.267667786  0.309153031  0.366419846  0.179549516  0.4951978720  0.3161076476  1.938132e-01
## nu       0.4884076328  0.071043603  0.242188660 -0.464991394 -0.056796960 -0.148294070 -0.380022625  0.3757094588 -0.0448491360  3.508063e-01
## nu      -0.0725982776 -0.044640657 -0.164516751 -0.592682537  0.361796728  0.175016964 -0.199134362 -0.0265040810 -0.4917389566 -3.936555e-01
## nu       0.2387622679 -0.074081781 -0.153294315  0.535306458 -0.045848268 -0.225912022  0.440938400 -0.3977049479  0.5377174163  4.689660e-02
## nu      -0.5630993493 -0.051079991  0.379751258  0.275613187 -0.161972940  0.009768852 -0.363260347 -0.1387325782 -0.2019722658  1.929397e-03
## nu       0.2268859142  0.171384068  0.271895286  0.155403802  0.114829553  0.071287114  0.385456290  0.2354492548 -0.0447158109  1.050749e-01
## nu       0.3955353340  0.059719586 -0.219517709 -0.257112701  0.603192821 -0.268363247 -0.487680580  0.4134374157 -0.1031027154 -9.039765e-02
## nu       0.1607591786 -0.064600765 -0.017061968  0.214532962  0.107696424  0.369062168  0.160627619  0.3941621439 -0.1617892930  2.988746e-01
## nu       0.2626793501  0.435876935  0.092842186 -0.372594208 -0.123670753  0.073463856 -0.115239659 -0.1181881618 -0.3716757209 -7.641015e-02
## nu       0.2689215136  0.222719778  0.242112380  0.530006343  0.287217119 -0.473470051 -0.171366423  0.4170348188  0.5190155767 -1.106151e-01
## nu      -0.0070056774 -0.271480018 -0.500985878 -0.478921721  0.202406810  0.037142960 -0.529163893  0.0040065974 -0.4608183467  9.511429e-02
## nu       0.9569463192  0.214860026  0.003227114 -0.036139120  0.046138123  0.626251976  0.005810970  0.5933334251 -0.0034100626  1.112437e-01
## nu       0.3560633325  0.098604433  0.165448492 -0.023018057 -0.227275429  0.262004605 -0.031231491 -0.3976011930 -0.0651563806 -4.150524e-01
## nu       0.5099054634  0.126209288 -0.146386800 -0.262717321 -0.042314389 -0.147450849  0.565274896  0.3768984053  0.2532250378 -1.234966e-01
## nu      -0.6440846266 -0.139839618 -0.146517707 -0.091275950  0.312166545 -0.442749598 -0.126290254 -0.0137272743 -0.2044594564 -3.016988e-03
## nu      -0.3693145506  0.193253657  0.203168545 -0.269857053  0.312808419  0.259238484  0.240947672  0.4825669333  0.0057232619  4.152826e-02
## nu      -0.3798436488 -0.709412264 -0.218051130 -0.111348139 -0.157300125 -0.387167380  0.195801194 -0.2605272171 -0.5103177496 -3.432773e-01
## nu       0.4791925253  0.146655426  0.361229691  0.608016563  0.542484963  0.543308615  0.056381039  0.1519687715  0.0296084437  3.547094e-01
## nu      -0.5697488629 -0.319182426 -0.238685373 -0.002265835 -0.109784604  0.203302293 -0.223163382 -0.0809142865  0.0620408005 -1.350070e-01
## nu       0.3262942572  0.434369731  0.266537040  0.836947469  0.137749356  0.737209244  0.659613466 -0.0484149985  0.4908889366  1.391962e-01
## nu       0.6836801114 -0.323437334  0.000487639  0.056335676 -0.298347788 -0.057590181  0.134101581  0.0148735471 -0.2494310680 -2.075757e-01
## nu      -0.6658015732  0.233258196  0.201944901 -0.606604800 -0.370589643  0.052016201 -1.052521147 -0.4685896492 -0.0330956190  1.409153e-02
## nu       0.2655154030  0.609949346  0.177545735  0.797782037  0.007545498 -0.090832098  0.136602008 -0.1806750449  0.2746319322  1.128102e-01
## nu       0.1116517278  0.389192154 -0.217948635 -0.110570171 -0.055261157 -0.115328192 -0.084098483 -0.0209824948  0.2990305781  1.256916e-01
## nu       0.5631139117  0.188025122  0.262110918  0.350714616  0.396831891 -0.034576426  0.111422210  0.2744818293  0.1244624102  6.979580e-01
## nu       0.1503445747 -0.100289843  0.295440139 -0.023861214  0.095833677  0.376338595 -0.145938130  0.3539166306  0.3368971355 -5.309718e-01
## nu      -0.4250973141 -0.755780340  0.328496578 -0.240682476 -0.115879549  0.226122329 -0.401210199 -0.4222324479 -0.6821443211 -2.699288e-01
## nu      -0.2279954722 -0.443409167 -0.272330665  0.216242760  0.178394857  0.587349221 -0.118470359  0.6723537254  0.1105925177  2.965567e-01
## nu       0.3924182084  0.594498596  0.606727934  0.928637799  0.439874220  0.298802424  0.623682454  0.6278193762  0.2088714740  6.739374e-01
## nu      -0.2499343862  0.027642086 -0.088699885 -0.060819733 -0.315383189  0.119799986  0.116825231 -0.1431023457 -0.4623471353 -1.028003e-01
## nu       0.0297291812  0.128414749 -0.343683681 -0.136103157 -0.208134241  0.071894178  0.273510791  0.2358591260 -0.2506064244  2.743557e-01
## nu      -0.4473859387 -0.498328406 -0.390494326 -0.547012703 -0.491122382 -0.615402932  0.256661443 -0.0730016968 -0.6870597706 -4.225929e-01
## nu      -0.2654299616 -0.224921393 -0.374108157 -0.035915933  0.447380773  0.239586186 -0.364664168  0.3617907179 -0.3629740691 -5.447794e-01
## nu       0.0371831004 -0.466639419  0.288213973  0.322809306  0.071798336  0.118786280 -0.290214258 -0.2504929225 -0.5307975233 -1.704196e-01
## nu       0.0329903079  0.274711768 -0.067579233 -0.176215556 -0.277817893 -0.378446972  0.093135190  0.0222631376 -0.6592716381 -2.071882e-01
## nu      -0.0001053712  0.396183917  0.103549531 -0.519677359  0.385715838 -0.615938386  0.062111637  0.4784877548 -0.1358845384 -1.754886e-02
## nu      -0.0766468473  0.649995113  0.590881534  0.477891598  0.368559738 -0.020339040  0.687616974  0.0588091196  0.3502844284  2.490478e-01
## nu       0.8126120732  0.371638837  0.207099280  0.441129922  0.668612353  0.133380465  0.504346647  0.3514197122 -0.0304403657  5.485302e-01
## nu      -0.3291068260 -0.107651892 -0.257152091 -0.329834602  0.263330573 -0.431831515  0.425155595  0.1625728472 -0.2752503267 -1.521985e-01
## nu      -0.5328483458 -0.539237211 -0.132511633 -0.359786668 -0.514055440 -0.195763249 -0.855226383 -0.7375405563 -0.5739580400 -1.296158e-01
## nu       0.3469401025  0.944797057  0.195636997  0.769927277  0.469188149  0.001364065  0.620100646  0.2124977777  0.5882258970  7.110342e-01
##  [ reached getOption("max.print") -- omitted 195 rows ]
#' This is now the length of the fixed effects
d <- length(quad$optresults$mode)

#' Just so that we can get the names of the theta (fixed effects)
thetanames <- colnames(quad$normalized_posterior$nodesandweights)[1:d]
theta <- simlist[j, thetanames]

#' In one dimension, R's indexing is not type consistent
if (!is.matrix(samps)) {
  samps <- rbind(samps)
  rownames(samps) <- NULL
}

if (!inherits(theta, "data.frame")) theta <- data.frame(theta1 = theta)

out <- list(
  samps = samps,
  theta = theta
)

out$theta
# Add the marginals for theta, with possible transformation
class(quad) <- "aghq" #' This is preventing sample_marginal.marginallaplace being used
out$thetasamples <- sample_marginal(quad, M, transformation, interpolation)
str(out)
## List of 3
##  $ samps       : num [1:295, 1:10] -0.1714 0.1106 -0.1415 -0.0338 0.2428 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:295] "epsilon" "epsilon" "epsilon" "epsilon" ...
##   .. ..$ : NULL
##  $ theta       :'data.frame':    10 obs. of  8 variables:
##   ..$ beta         : num [1:10] 1.51 1.51 1.65 1.65 1.65 ...
##   ..$ beta1        : num [1:10] -1.36 -1.36 -1.32 -1.32 -1.32 ...
##   ..$ beta2        : num [1:10] 0.708 0.708 0.898 0.692 0.898 ...
##   ..$ beta3        : num [1:10] -0.1894 -0.0171 -0.1892 -0.1886 -0.0168 ...
##   ..$ beta4        : num [1:10] 0.862 0.194 0.916 0.851 0.249 ...
##   ..$ beta5        : num [1:10] 0.519 0.484 0.453 0.637 0.418 ...
##   ..$ l_tau_epsilon: num [1:10] 1.8 1.23 1.31 1.87 1.29 ...
##   ..$ l_tau_nu     : num [1:10] 1.78 1.82 2.37 1.83 2.34 ...
##  $ thetasamples:List of 8
##   ..$ : num [1:10] 1.48 1.46 1.56 1.56 1.53 ...
##   ..$ : num [1:10] -0.405 -1.708 -1.193 -0.751 -0.381 ...
##   ..$ : num [1:10] 0.865 1.078 0.765 0.862 0.781 ...
##   ..$ : num [1:10] -0.04771 -0.04978 0.04304 -0.10947 0.00232 ...
##   ..$ : num [1:10] 0.306 1.017 1.14 -0.116 1.047 ...
##   ..$ : num [1:10] 0.1153 0.0511 0.417 0.2293 0.3849 ...
##   ..$ : num [1:10] 1.82 1.58 1.17 1.77 1.33 ...
##   ..$ : num [1:10] 2.54 2.03 2.29 2.06 2.26 ...

aghq::sample_marginal.aghq

Let’s have a look at what that call to sample_marginal.aghq is doing:

class(quad)
## [1] "aghq"
out <- list()

#' If the quadrature object somehow doesn't contain any information about the 
if (is.null(quad$marginals)) return(out)

str(quad$marginals)
## List of 8
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta1     : num [1:2] 1.51 1.65
##   ..$ logmargpost: num [1:2] 1.24 1.14
##   ..$ w          : num [1:2] 0.152 0.152
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta2     : num [1:2] -0.553 -1.344
##   ..$ logmargpost: num [1:2] -0.503 -0.48
##   ..$ w          : num [1:2] 0.817 0.817
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta3     : num [1:2] 0.751 1.01
##   ..$ logmargpost: num [1:2] 0.621 0.633
##   ..$ w          : num [1:2] 0.267 0.267
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta4     : num [1:2] -0.0169 -0.1893
##   ..$ logmargpost: num [1:2] 1.03 1.04
##   ..$ w          : num [1:2] 0.178 0.178
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta5     : num [1:2] 0.148 0.832
##   ..$ logmargpost: num [1:2] -0.341 -0.353
##   ..$ w          : num [1:2] 0.707 0.707
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta6     : num [1:2] 0.149 0.55
##   ..$ logmargpost: num [1:2] 0.182 0.195
##   ..$ w          : num [1:2] 0.414 0.414
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta7     : num [1:2] 1.28 1.83
##   ..$ logmargpost: num [1:2] 0.235 -0.684
##   ..$ w          : num [1:2] 0.565 0.565
##  $ :'data.frame':    2 obs. of  3 variables:
##   ..$ theta8     : num [1:2] 1.82 2.3
##   ..$ logmargpost: num [1:2] 0.0814 -0.0588
##   ..$ w          : num [1:2] 0.493 0.493

As usual, we’re going to start by doing the computation for the first element of theta explicitly:

i <- 1

We already have information about each marginal, contained in:

quad$marginals[[i]]

So we can use the function aghq::compute_quantiles to, my guess, use the inverse CDF method to sample from the marginal, given input of some \(\mathcal{U}(0, 1)\) draws:

unif <- stats::runif(M)
unif
##  [1] 0.50467636 0.73696872 0.05496506 0.80764065 0.14720269 0.04936790 0.59864599 0.77080723 0.15874619 0.14437648
compute_quantiles(quad$marginals[[i]], unif, interpolation = interpolation, transformation = transformation)
## 50.4676361335441% 73.6968720331788% 5.49650606699288% 80.7640645187348% 14.7202692227438%  4.9367897445336% 59.8645987687632% 77.0807233639061%  15.874619083479% 
##          1.577305          1.649638          1.446752          1.672279          1.472627          1.445281          1.606121          1.660518          1.475861 
##  14.437648258172% 
##          1.471745
for (i in 1:length(quad$marginals)) {
  out[[i]] <- unname(compute_quantiles(quad$marginals[[i]], stats::runif(M), interpolation = interpolation, transformation = transformation))
} 

#' The function doesn't name these yet, but in future it could (and also return a data.frame) via
d <- length(quad$optresults$mode)
thetanames <- colnames(quad$normalized_posterior$nodesandweights)[1:d]
names(out) <- thetanames
as.data.frame(out)

aghq::compute_quantiles.aghq

Does the compute_quantiles function work as we expect it to?

class(quad$marginals[[i]])
## [1] "data.frame"

The PDF and CDF are computed by compute_pdf_and_cdf (as the name suggests, and another layer deep that we can look). Code comments note that transformations are not applied at this point, and that it’s better to transform the quantiles later:

q <- stats::runif(M)

pdf_and_cdf <- compute_pdf_and_cdf(quad$marginals[[i]], interpolation = interpolation)

ggplot(pdf_and_cdf, aes(x = theta, y = cdf)) +
  geom_line() +
  theme_minimal()

out <- numeric(length(q))
increasing <- TRUE

validate_transformation(transformation)
## [1] TRUE
transformation <- make_transformation(transformation)

For each sample, take the maximum theta which still has CDF value less than the quantile desired:

for (i in 1:length(q)) out[i] <- pdf_and_cdf$theta[max(which(pdf_and_cdf$cdf < q[i]))]

#' Have to check if it's increasing or decreasing so can reverse order if necessary
increasing <- transformation$fromtheta(min(out)) <= transformation$fromtheta(max(out))
for (i in 1:length(out)) out[i] <- transformation$fromtheta(out[i])

if (!increasing) out <- rev(out)
names(out) <- paste0(as.character(100 * q), "%")
out
## 30.6657020933926% 58.0539908260107%  19.928448763676% 10.3568355785683% 27.6170151075348% 14.8865669965744% 31.1338842846453% 88.6961564188823% 92.8852251498029% 
##          1.856470          2.121179          1.757085          1.672034          1.827801          1.712171          1.860293          2.445136          2.491962 
## 82.9638364957646% 
##          2.382065

compute_pdf_and_cdf

Going to leave this as a to-do.