1 2 3 4 5 6 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 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 91 92 93 94 95 96 97 98 99 100
|
#' @rdname p_value_wald
#'
#' @param ci Confidence Interval (CI) level. Default to 0.95 (95\%).
#' @param dof Degrees of Freedom. If not specified, for \code{ci_wald()}, defaults to model's residual degrees of freedom (i.e. \code{n-k}, where \code{n} is the number of observations and \code{k} is the number of parameters). For \code{p_value_wald()}, defaults to \code{Inf}.
#'
#' @inheritParams simulate_model
#' @inheritParams standard_error
#' @inheritParams model_parameters.default
#'
#' @importFrom stats qt coef
#' @export
ci_wald <- function(model, ci = .95, dof = NULL, effects = c("fixed", "random", "all"), component = c("all", "conditional", "zi", "zero_inflated", "dispersion", "precision", "scale", "smooth_terms", "full", "marginal"), robust = FALSE, ...) {
effects <- match.arg(effects)
component <- match.arg(component)
out <- lapply(ci, function(i) {
.ci_wald(model = model, ci = i, dof = dof, effects = effects, component = component, robust = robust, method = "wald", ...)
})
out <- do.call(rbind, out)
row.names(out) <- NULL
out
}
#' @importFrom insight get_parameters n_obs
#' @importFrom stats qt complete.cases
#' @keywords internal
.ci_wald <- function(model, ci, dof, effects, component, robust = FALSE, method = "wald", se = NULL, ...) {
if (inherits(model, "emmGrid")) {
params <- insight::get_parameters(model, effects = effects, component = component, merge_parameters = TRUE)
} else {
params <- insight::get_parameters(model, effects = effects, component = component)
}
# check if all estimates are non-NA
params <- .check_rank_deficiency(params, verbose = FALSE)
method <- tolower(method)
# if we have adjusted SE, e.g. from kenward-roger, don't recompute
# standard errors to save time...
if (is.null(se)) {
stderror <- if (isTRUE(robust)) {
standard_error_robust(model, ...)
} else {
switch(
method,
"wald" = standard_error(model, component = component),
"kenward" = ,
"kr" = se_kenward(model),
"ml1" = se_ml1(model),
"betwithin" = se_betwithin(model),
"satterthwaite" = se_satterthwaite(model),
standard_error(model, component = component)
)
}
if (.is_empty_object(stderror)) {
return(NULL)
}
# filter non-matching parameters
if (nrow(stderror) != nrow(params)) {
params <- stderror <- merge(stderror, params, sort = FALSE)
}
se <- stderror$SE
}
if (is.null(dof)) {
# residual df
dof <- degrees_of_freedom(model, method = "any")
# make sure we have a value for degrees of freedom
if (is.null(dof) || length(dof) == 0) {
dof <- Inf
} else if (length(dof) > nrow(params)) {
# filter non-matching parameters
dof <- dof[1:nrow(params)]
}
}
alpha <- (1 + ci) / 2
fac <- suppressWarnings(stats::qt(alpha, df = dof))
out <- cbind(
CI_low = params$Estimate - se * fac,
CI_high = params$Estimate + se * fac
)
out <- as.data.frame(out)
out$CI <- ci * 100
out$Parameter <- params$Parameter
out <- out[c("Parameter", "CI", "CI_low", "CI_high")]
if ("Component" %in% names(params)) out$Component <- params$Component
if ("Effects" %in% names(params) && effects != "fixed") out$Effects <- params$Effects
if (anyNA(params$Estimate)) {
out[stats::complete.cases(out), ]
} else {
out
}
}
|