pedmod_profile_nleq {pedmod} | 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 |
par |
numeric vector with the maximum likelihood estimator e.g. from
|
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
|
rel_eps |
rel_eps convergence threshold for
|
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
|
delta |
numeric scalar with an initial step to take. Subsequent steps
are taken by |
indices |
zero-based vector with indices of which log marginal
likelihood terms to include. Use |
maxvls_start , minvls_start |
number of samples to use when finding the initial values for the optimization. |
do_reorder |
|
use_aprx |
|
n_threads |
number of threads to use. |
cluster_weights |
numeric vector with weights for each cluster. Use
|
method |
integer with the method to use. Zero yields randomized Korobov lattice rules while one yields scrambled Sobol sequences. |
seed |
seed to pass to |
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 |
|
vls_scales |
can be a numeric vector with a positive scalar for each
cluster. Then |
control.outer , control.optim , ... |
arguments passed to
|
See Also
pedmod_opt
, pedmod_sqn
,
pedmod_profile
, and pedmod_profile_prop
.
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)