# R/ci_f_ncp.R In confintr: Confidence Intervals

#### Documented in ci_f_ncp

#' Confidence Interval for the Non-Centrality Parameter of the F Distribution
#'
#' Based on the inversion principle, parametric confidence intervals for the non-centrality parameter Delta of the F distribution are calculated. Note that we do not provide bootstrap confidence intervals here to keep the input interface simple. A positive lower (1-alpha)*100%-confidence limit for the ncp goes hand-in-hand with a significant F test at level alpha.
#'
#' Note that, according to \code{?pf}, the results might be unreliable for very large F values.
#' @importFrom stats lm pf uniroot
#' @param x The result of \code{lm} or the F test statistic.
#' @param df1 The numerator degree of freedom, e.g. the number of parameters (including the intercept) of a linear regression. Only used if \code{x} is a test statistic.
#' @param df2 The denominator degree of freedom, e.g. n - df1 - 1 in a linear regression. Only used if \code{x} is a test statistic.
#' @param probs Error probabilites. The default c(0.025, 0.975) gives a symmetric 95% confidence interval.
#' @return A list with class \code{cint} containing these components:
#' \itemize{
#'   \item \code{parameter}: The parameter in question.
#'   \item \code{interval}: The confidence interval for the parameter.
#'   \item \code{estimate}: The estimate for the parameter.
#'   \item \code{probs}: A vector of error probabilities.
#'   \item \code{type}: The type of the interval.
#'   \item \code{info}: An additional description text for the interval.
#' }
#' @export
#' @examples
#' fit <- lm(Sepal.Length ~ ., data = iris)
#' ci_f_ncp(fit)
#' ci_f_ncp(fit, probs = c(0.05, 1))
#' ci_f_ncp(fit, probs = c(0, 0.95))
#' ci_f_ncp(x = 188.251, df1 = 5, df2 = 144)
#' @references
#' Smithson, M. (2003). Confidence intervals. Series: Quantitative Applications in the Social Sciences. New York, NY: Sage Publications.
ci_f_ncp <- function(x, df1 = NULL, df2 = NULL, probs = c(0.025, 0.975)) {
# Input checks and initialization
check_probs(probs)
iprobs <- 1 - probs
limits <- c(0, Inf)
stopifnot(inherits(x, "lm") || is.numeric(x))

# Distinguish input
if (inherits(x, "lm")) {
sx <- summary(x)
stopifnot("fstatistic" %in% names(sx))
fstat <- sx[["fstatistic"]]
stat <- fstat[["value"]]
df1 <- fstat[["numdf"]]
df2 <- fstat[["dendf"]]
}
if (is.numeric(x)) {
stopifnot(length(x) == 1L,
!is.null(df1),
!is.null(df2))
stat <- x
}

# Estimate
estimate <- f_to_ncp(stat, df1, df2)

# Calculate limits
if (probs[1] == 0) {
lci <- limits[1]
} else {
lci <- try(uniroot(function(ncp) pf(stat, df1 = df1, df2 = df2, ncp = ncp) - iprobs[1],
interval = c(0, estimate))$root, silent = TRUE) if (inherits(lci, "try-error")) { lci <- limits[1] } } if (probs[2] == 1) { uci <- limits[2] } else { # Upper limit might be improved upper_limit <- pmax(4 * estimate, stat * df1 * 4, df1 * 100) uci <- try(uniroot(function(ncp) pf(stat, df1 = df1, df2 = df2, ncp = ncp) - iprobs[2], interval = c(estimate, upper_limit))$root,
silent = TRUE)
if (inherits(uci, "try-error")) {
warning("Upper limit outside search range. Set to the maximum of the parameter range.")
uci <- limits[2]
}
}

# Organize output
cint <- check_output(c(lci, uci), probs, limits)
out <- list(parameter = "non-centrality parameter of the F-distribution",
interval = cint, estimate = estimate,
probs = probs, type = "F", info = "")
class(out) <- "cint"
out
}


## Try the confintr package in your browser

Any scripts or data that you put into this service are public.

confintr documentation built on Jan. 29, 2022, 1:08 a.m.