pedmod (version 0.2.4)

eval_pedigree_ll: Approximate the Log Marginal Likelihood

Description

Approximate the log marginal likelihood and the derivatives with respect to the model parameters.

Usage

eval_pedigree_ll(
  ptr,
  par,
  maxvls,
  abs_eps,
  rel_eps,
  indices = NULL,
  minvls = -1L,
  do_reorder = TRUE,
  use_aprx = FALSE,
  n_threads = 1L,
  cluster_weights = NULL,
  standardized = FALSE,
  method = 0L,
  use_tilting = FALSE,
  vls_scales = NULL
)
eval_pedigree_grad(
  ptr,
  par,
  maxvls,
  abs_eps,
  rel_eps,
  indices = NULL,
  minvls = -1L,
  do_reorder = TRUE,
  use_aprx = FALSE,
  n_threads = 1L,
  cluster_weights = NULL,
  standardized = FALSE,
  method = 0L,
  use_tilting = FALSE,
  vls_scales = NULL
)
eval_pedigree_hess(
  ptr,
  par,
  maxvls,
  abs_eps,
  rel_eps,
  indices = NULL,
  minvls = -1L,
  do_reorder = TRUE,
  use_aprx = FALSE,
  n_threads = 1L,
  cluster_weights = NULL,
  standardized = FALSE,
  method = 0L,
  use_tilting = FALSE,
  vls_scales = NULL
)

Value

eval_pedigree_ll: a scalar with the log marginal likelihood approximation. It has an attribute called "n_fails" which shows the number of log marginal likelihood term approximations which do not satisfy the abs_eps and rel_eps criteria and an attribute called std with a standard error estimate based on the delta rule.

eval_pedigree_grad: a vector with the derivatives with respect to par. An attribute called "logLik" contains the log marginal likelihood approximation. There will also be "n_fails"

attribute like for eval_pedigree_ll and an attribute called "std" which first element is the standard error estimate of the log likelihood based on the delta method and the last elements are the standard error estimates of the gradient. The latter ignores the Monte Carlo error from the likelihood approximation.

eval_pedigree_hess: a matrix with the hessian with respect to par. An attribute called "logLik" contains the log marginal likelihood approximation and an attribute called "grad"

contains the gradient· The attribute "hess_org" contains the Hessian with the scale parameter on the identity scale rather than the log scale. "vcov" and "vcov_org" are the covariance matrices from the hessian and "hess_org".

Arguments

ptr

object from pedigree_ll_terms or pedigree_ll_terms_loadings.

par

numeric vector with parameters. For an object from pedigree_ll_terms these are the fixed effect coefficients and log scale parameters. The log scale parameters should be last. For an object from pedigree_ll_terms_loadings these are the fixed effects and the coefficients for scale parameters.

maxvls

maximum number of samples in the approximation for each marginal likelihood term.

abs_eps

absolute convergence threshold.

rel_eps

relative convergence threshold.

indices

zero-based vector with indices of which log marginal likelihood terms to include. Use NULL if all indices should be used.

minvls

minimum number of samples for each marginal likelihood term. Negative values provides a default which depends on the dimension of the integration.

do_reorder

TRUE if a heuristic variable reordering should be used. TRUE is likely the best value.

use_aprx

TRUE if a less precise approximation of pnorm and qnorm should be used. This may reduce the computation time while not affecting the result much.

n_threads

number of threads to use.

cluster_weights

numeric vector with weights for each cluster. Use NULL if all clusters have weight one.

standardized

logical for whether to use the standardized or direct parameterization. See standardized_to_direct and the vignette at vignette("pedmod", package = "pedmod").

method

integer with the method to use. Zero yields randomized Korobov lattice rules while one yields scrambled Sobol sequences.

use_tilting

TRUE if the minimax tilting method suggested by Botev (2017) should be used. See tools:::Rd_expr_doi("10.1111/rssb.12162").

vls_scales

can be a numeric vector with a positive scalar for each cluster. Then vls_scales[i] * minvls and vls_scales[i] * maxvls is used for cluster i rather than minvls and maxvls. Set vls_scales = NULL if the latter should be used.

Details

eval_pedigree_hess is only implemented for objects from pedigree_ll_terms.

Examples

# three families as an example
fam_dat <- list(
  list(
    y = c(FALSE, TRUE, FALSE, FALSE),
    X = structure(c(
      1, 1, 1, 1, 1.2922654151273, 0.358134905909256, -0.734963997107464,
      0.855235473516044, -1.16189500386223, -0.387298334620742,
      0.387298334620742, 1.16189500386223),
      .Dim = 4:3, .Dimnames = list( NULL, c("(Intercept)", "X1", ""))),
    rel_mat = structure(c(
      1, 0.5, 0.5, 0.125, 0.5, 1, 0.5, 0.125, 0.5, 0.5,
      1, 0.125, 0.125, 0.125, 0.125, 1), .Dim = c(4L, 4L)),
    met_mat = structure(c(1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1),
                        .Dim = c(4L, 4L))),
  list(
    y = c(FALSE, FALSE, FALSE),
    X = structure(c(
      1, 1, 1, -0.0388728997202442, -0.0913782435233639,
      -0.0801619722392612, -1, 0, 1), .Dim = c(3L, 3L)),
    rel_mat = structure(c(
      1, 0.5, 0.125, 0.5, 1, 0.125, 0.125, 0.125, 1), .Dim = c(3L, 3L)),
    met_mat = structure(c(
      1, 1, 0, 1, 1, 0, 0, 0, 1), .Dim = c(3L, 3L))),
  list(
    y = c(TRUE, FALSE),
    X = structure(c(
      1, 1, 0.305275750370738, -1.49482995913648,  -0.707106781186547,
      0.707106781186547),
      .Dim = 2:3, .Dimnames = list( NULL, c("(Intercept)", "X1", ""))),
    rel_mat = structure(c(1, 0.5,  0.5, 1), .Dim = c(2L, 2L)),
    met_mat = structure(c(1, 1, 1, 1), .Dim = c(2L,  2L))))

# get the data into the format needed for the package
dat_arg <- lapply(fam_dat, function(x){
  # we need the following for each family:
  #   y: the zero-one outcomes.
  #   X: the design matrix for the fixed effects.
  #   scale_mats: list with the scale matrices for each type of effect.
  list(y = as.numeric(x$y), X = x$X,
       scale_mats = list(x$rel_mat, x$met_mat))
})

# get a pointer to the C++ object
ptr <- pedigree_ll_terms(dat_arg, max_threads = 1L)

# approximate the log marginal likelihood
beta <- c(-1, 0.3, 0.2) # fixed effect coefficients
scs <- c(0.5, 0.33)     # scales parameters

set.seed(44492929)
system.time(ll1 <- eval_pedigree_ll(
  ptr = ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e5,
  rel_eps = 1e-5, minvls = 2000, use_aprx = FALSE))
ll1 # the approximation

# with the approximation of pnorm and qnorm
system.time(ll2 <- eval_pedigree_ll(
  ptr = ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e5,
  rel_eps = 1e-5, minvls = 2000, use_aprx = TRUE))
all.equal(ll1, ll2, tolerance = 1e-5)

# cluster weights can be used as follows to repeat the second family three
# times and remove the third
system.time(deriv_w_weight <- eval_pedigree_grad(
  ptr = ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e6,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE,
  cluster_weights = c(1, 3, 0)))

# the same as manually repeating second cluster and not including the third
dum_dat <- dat_arg[c(1, 2, 2, 2)]
dum_ptr <- pedigree_ll_terms(dum_dat, 1L)
system.time(deriv_dum <- eval_pedigree_grad(
  ptr = dum_ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e6,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE))
all.equal(deriv_dum, deriv_w_weight, tolerance = 1e-3)

# the hessian is computed on the scale parameter scale rather than on the
# log of the scale parameters
system.time(hess_w_weight <- eval_pedigree_hess(
  ptr = ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e6,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE,
  cluster_weights = c(1, 3, 0)))

system.time(hess_dum <- eval_pedigree_hess(
  ptr = dum_ptr, par = c(beta, log(scs)), abs_eps = -1, maxvls = 1e6,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE))
attr(hess_w_weight, "n_fails") <- attr(hess_dum, "n_fails") <- NULL
all.equal(hess_w_weight, hess_dum, tolerance = 1e-3)

# the results are consistent with the gradient output
all.equal(attr(deriv_dum, "logLik"), attr(hess_dum, "logLik"),
          tolerance = 1e-5)

hess_grad <- attr(hess_dum, "grad")
all.equal(hess_grad, deriv_dum, check.attributes = FALSE,
          tolerance = 1e-3)

# with loadings
dat_arg_loadings <- lapply(fam_dat, function(x){
  list(y = as.numeric(x$y), X = x$X, Z = x$X[, 1:2],
       scale_mats = list(x$rel_mat, x$met_mat))
})

ptr_loadings <-
  pedigree_ll_terms_loadings(dat_arg_loadings, max_threads = 1L)

scs <- c(log(0.5) / 2, 0.1, log(0.33) / 2, 0.2) # got more scales parameters
eval_pedigree_ll(
  ptr = ptr_loadings, par = c(beta, scs), abs_eps = -1, maxvls = 1e4,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE)
eval_pedigree_grad(
  ptr = ptr_loadings, par = c(beta, scs), abs_eps = -1, maxvls = 1e4,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE)

# can recover the result from before
scs <- c(log(0.5) / 2, 0, log(0.33) / 2, 0)
ll3 <- eval_pedigree_ll(
  ptr = ptr_loadings, par = c(beta, scs), abs_eps = -1, maxvls = 1e4,
  rel_eps = 1e-3, minvls = 2000, use_aprx = TRUE)
all.equal(ll1, ll3, tolerance = 1e-5)