pedmod_profile_nleq: Computes Profile Likelihood Based Confidence Intervals for a...
In pedmod: Pedigree Models

pedmod_profile_nleq

R Documentation

Computes Profile Likelihood Based Confidence Intervals for a Nonlinear Transformation of the Variables

Description

Computes Profile Likelihood Based Confidence Intervals for a Nonlinear Transformation of the Variables

Usage

pedmod_profile_nleq(
  ptr,
  par,
  maxvls,
  minvls = -1L,
  alpha = 0.05,
  abs_eps,
  rel_eps,
  heq,
  heq_bounds = c(-Inf, Inf),
  delta,
  indices = NULL,
  maxvls_start = max(100L, as.integer(ceiling(maxvls/5))),
  minvls_start = if (minvls < 0) minvls else minvls/5,
  do_reorder = TRUE,
  use_aprx = FALSE,
  n_threads = 1L,
  cluster_weights = NULL,
  method = 0L,
  seed = 1L,
  verbose = FALSE,
  max_step = 15L,
  use_tilting = FALSE,
  vls_scales = NULL,
  control.outer = list(itmax = 100L, method = "BFGS", kkt2.check = FALSE, trace =
    FALSE),
  control.optim = list(fnscale = get_n_terms(ptr)),
  ...
)

Arguments

`ptr`	object from `pedigree_ll_terms` or `pedigree_ll_terms_loadings`.
`par`	numeric vector with the maximum likelihood estimator e.g. from `pedmod_opt`.
`maxvls`	maximum number of samples in the approximation for each marginal likelihood term.
`minvls`	minimum number of samples for each marginal likelihood term. Negative values provides a default which depends on the dimension of the integration.
`alpha`	numeric scalar with the confidence level required.
`abs_eps`	absolute convergence threshold for `eval_pedigree_ll` and `eval_pedigree_grad`.
`rel_eps`	rel_eps convergence threshold for `eval_pedigree_ll` and `eval_pedigree_grad`.
`heq`	function that returns a one dimensional numerical vector which should be profiled. It does not need to evaluate to zero at the maximum likelihood estimator.
`heq_bounds`	two dimensional numerical vector with bounds for `heq`.
`delta`	numeric scalar with an initial step to take. Subsequent steps are taken by `2^(<iteration number> - 1) * delta`. Two times the standard error is a good value or a guess thereof. Hessian approximations are not implemented as of this writing and therefore the user needs to provide some guess.
`indices`	zero-based vector with indices of which log marginal likelihood terms to include. Use `NULL` if all indices should be used.
`maxvls_start, minvls_start`	number of samples to use when finding the initial values for the optimization.
`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.
`method`	integer with the method to use. Zero yields randomized Korobov lattice rules while one yields scrambled Sobol sequences.
`seed`	seed to pass to `set.seed` before each gradient and function evaluation. Use `NULL` if the seed should not be fixed.
`verbose`	logical for whether output should be printed to the console during the estimation of the profile likelihood curve.
`max_step`	integer scalar with the maximum number of steps to take in either directions.
`use_tilting`	`TRUE` if the minimax tilting method suggested by Botev (2017) should be used. See 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.
`control.outer, control.optim, ...`	arguments passed to `auglag`

Examples


# similar examples to that in help("pedmod_profile_prop")
K <- matrix(c(
  0.5  , 0    , 0.25 , 0   , 0.25 , 0   , 0.125 , 0.125 , 0.125 , 0.125 ,
  0    , 0.5  , 0.25 , 0   , 0.25 , 0   , 0.125 , 0.125 , 0.125 , 0.125 ,
  0.25 , 0.25 , 0.5  , 0   , 0.25 , 0   , 0.25  , 0.25  , 0.125 , 0.125 ,
  0    , 0    , 0    , 0.5 , 0    , 0   , 0.25  , 0.25  , 0     , 0     ,
  0.25 , 0.25 , 0.25 , 0   , 0.5  , 0   , 0.125 , 0.125 , 0.25  , 0.25  ,
  0    , 0    , 0    , 0   , 0    , 0.5 , 0     , 0     , 0.25  , 0.25  ,
  0.125, 0.125, 0.25 , 0.25, 0.125, 0   , 0.5   , 0.25  , 0.0625, 0.0625,
  0.125, 0.125, 0.25 , 0.25, 0.125, 0   , 0.25  , 0.5   , 0.0625, 0.0625,
  0.125, 0.125, 0.125, 0   , 0.25 , 0.25, 0.0625, 0.0625, 0.5   , 0.25  ,
  0.125, 0.125, 0.125, 0   , 0.25 , 0.25, 0.0625, 0.0625, 0.25  , 0.5
), 10)

C <- matrix(c(
  1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
  0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
  0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
  0, 0, 0, 0, 0, 0, 0, 0, 1, 1
), 10L)

# simulates a data set.
#
# Args:
#   n_fams: number of families.
#   beta: the fixed effect coefficients.
#   sig_sq: the scale parameters.
sim_dat <- function(n_fams, beta = c(-1, 1, 2), sig_sq = c(3, 1)){
  # setup before the simulations
  Cmat <- 2 * K
  n_obs <- NROW(K)
  Sig <- diag(n_obs) + sig_sq[1] * Cmat + sig_sq[2] * C
  Sig_chol <- chol(Sig)

  # simulate the data
  out <- replicate(
    n_fams, {
      # simulate covariates
      X <- cbind(`(Intercept)` = 1, Continuous = rnorm(n_obs),
                 Binary = runif(n_obs) > .5)

      # assign the linear predictor + noise
      eta <- drop(X %*% beta) + drop(rnorm(n_obs) %*% Sig_chol)

      # return the list in the format needed for the package
      list(y = as.numeric(eta > 0), X = X,
           scale_mats = list(genetic = Cmat, environment = C))
    }, simplify = FALSE)

  # add attributes with the true values and return
  attributes(out) <- list(beta = beta, sig_sq = sig_sq)
  out
}

# simulate the data
set.seed(1)
dat <- sim_dat(200L)

# fit the model
ptr <- pedigree_ll_terms(dat, max_threads = 2L)
start <- pedmod_start(ptr = ptr, data = dat, n_threads = 2L)
fit <- pedmod_opt(ptr = ptr, par = start$par, use_aprx = TRUE, n_threads = 2L,
                  maxvls = 5000L, minvls = 1000L, abs_eps = 0, rel_eps = 1e-3)
fit$par # the estimate

# 90% likelihood ratio based confidence interval for the proportion of variance
# of the genetic effect
heq <- function(par){
 vars <- exp(tail(par, 2))
 vars[1] / (1 + sum(vars))
}
heq(fit$par)
prof_out_nleq <- pedmod_profile_nleq(
  ptr = ptr, fit$par, maxvls = 2500L, minvls = 500L, alpha = .1,
  abs_eps = 0, rel_eps = 1e-3, verbose = TRUE, use_aprx = TRUE,
  heq = heq, heq_bounds = c(0, Inf), delta = .2, n_threads = 2L)
prof_out_nleq$confs # the confidence interval for the proportion

# plot the log profile likelihood
plot(prof_out_nleq$xs, prof_out_nleq$p_log_Lik, pch = 16,
     xlab = "proportion of variance", ylab = "log profile likelihood")
abline(v = prof_out_nleq$confs, lty = 2)

pedmod documentation built on Sept. 11, 2022, 5:05 p.m.