# R/cpciuupi2.R In ciuupi2: Kabaila and Giri (2009) Confidence Interval

#### Documented in cpciuupi2

#'Compute the coverage probability of the Kabaila & Giri (2009) CIUUPI
#'
#'Evaluate the coverage probability of the Kabaila & Giri (2009) confidence
#'interval that utilizes uncertain prior information (CIUUPI),
#'with minimum coverage 1 - \code{alpha}, at \code{gam}.
#'
#'@param gam A value of gamma or vector of gamma values at which the coverage
#'  probability function is evaluated
#'@param bsvec The vector (b(d/6),b(2d/6),...,b(5d/6),s(0),s(d/6),...,s(5d/6))
#'  computed using \code{bsciuupi2}
#'@param alpha The minimum coverage probability is 1 - \code{alpha}
#'@param m Degrees of freedom \code{n - p}
#'@param rho A known correlation
#'@param natural Equal to 1 (default) if the b and s functions are obtained by
#'  natural cubic spline interpolation or 0 if obtained by clamped cubic spline
#'  interpolation. This parameter must take the same value as that used in
#'  \code{bsciuupi2}
#'
#'@return The value(s) of the coverage probability of the Kabaila & Giri (2009)
#'  CIUUPI at \code{gam}.
#'
#'@details Suppose that \deqn{y = X \beta + \epsilon} where \eqn{y} is a random
#'  \eqn{n}-vector of responses, \eqn{X} is a known \eqn{n} by \eqn{p} matrix
#'  with linearly independent columns, \eqn{\beta} is an unknown parameter
#'  \eqn{p}-vector and \eqn{\epsilon} is a random \eqn{n}-vector with
#'  components that are independent and identically normally distributed with
#'  zero mean and unknown variance. The parameter of interest is \eqn{\theta = }
#'  \code{a}' \eqn{\beta}. The uncertain prior information is that \eqn{\tau = }
#'  \code{c}' \eqn{\beta} takes the value \code{t}, where \code{a} and \code{c}
#'  are specified linearly independent vectors and \code{t} is a specified
#'  number. \code{rho} is the known correlation between the least squares
#'  estimators of \eqn{\theta} and \eqn{\tau}. It is determined by the \eqn{n}
#'  by \eqn{p} design matrix X and the \eqn{p}-vectors a and c using
#'
#'  In the examples, we continue with the same 2 x 2 factorial example described
#'  in the documentation for \code{\link{find_rho}}.
#'
#'
#' @examples
#' alpha <- 0.05
#' m <- 8
#'
#' # Find the vector (b(d/6),...,b(5d/6),s(0),...,s(5d/6)) that specifies the
#' # Kabaila & Giri (2009) CIUUPI for the first definition of the
#' # scaled expected length (default) (takes about 30 mins to run):
#' \donttest{
#' bsvec <- bsciuupi2(alpha, m, rho = -0.7071068)
#' }
#'
#' # The result bsvec is (to 7 decimal places) the following:
#' bsvec <- c(-0.0287487, -0.2151595, -0.3430403, -0.3125889, -0.0852146,
#'             1.9795390,  2.0665414,  2.3984471,  2.6460159,  2.6170066,  2.3925494)
#'
#'
#' # Graph the coverage probability function
#' gam <- seq(0, 10, by = 0.1)
#' cp <- cpciuupi2(gam, bsvec, alpha, m, rho = -0.7071068)
#' plot(gam, cp, type = "l", lwd = 2, ylab = "", las = 1, xaxs = "i",
#' main = "Coverage Probability", col = "blue",
#' xlab = expression(paste("|", gamma, "|")), ylim = c(0.9490, 0.9510))
#' abline(h = 1-alpha, lty = 2)
#'
#'@references
#'
#'Kabaila, P. and Giri, K. (2009) Confidence intervals in regression utilizing
#'prior information.  Journal of Statistical Planning and Inference, 139, 3419 -
#'3429.
#'
#'@export

cpciuupi2 <- function(gam, bsvec, alpha, m, rho, natural = 1){

# Specify the values of the inputs to other functions
n.ints <- 6
N <- 33
eps <- 10^{-10}
n.nodes <- 20
d <- choice_d(m)$d # Set t.alpha for m and alpha t.alpha <- stats::qt(1 - alpha/2, m) # Find the nodes and weights of the legendre quadrature quad.info <- statmod::gauss.quad(n.nodes, kind="legendre") nodes <- quad.info$nodes
weights <- quad.info$weights # Specify where the knots for b and s are located # as inputs to other functions knots <- seq(0, d, by = d/n.ints) knots.all <- seq(-d, d, by = d/n.ints) # Specify the values of the inputs to compute the # outer integral in coverage probability h <- OuterPara(m, nu=m+1, N, eps)$h
wvec <- OuterPara(m, nu=m+1, N, eps)$wvec psinu.zvec <- OuterPara(m, nu=m+1, N, eps)$psinu.zvec
cons <- OuterPara(m, nu=m+1, N, eps)\$cons1

# Find the b and s functions
b.spl <- spline_b(bsvec, n.ints, knots.all, t.alpha, natural)
s.spl <- spline_s(bsvec, n.ints, knots.all, t.alpha, natural)

# Compute the coverage probability
res <- rep(0, length(gam))
for(i in 1:length(gam)){
res[i] <- compute_cov_trapez(gam[i], rho, knots, alpha,
t.alpha, nodes, weights, b.spl, s.spl,
wvec, psinu.zvec, h, cons)
}

out <- res

}


## Try the ciuupi2 package in your browser

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

ciuupi2 documentation built on March 11, 2021, 5:06 p.m.