Iboot (version 0.1-1)

Iboot: Iterated bootstrap tests and confidence sets for a $p$-dimensional parameter

Description

This function provides both simple and re-calibrated bootstrap critical values for tests and confidence sets for a $p$-dimensional parameter, obtained by a single weighted bootstrap iteration. The considered statistic is the unstudentised version of the Rao statistic as proposed in Lunardon (2013). Details about the implemented algorithm can be found in Nankervis (2005) and Lunardon (2013).

Usage

Iboot(gradient, B=500, M=500, kB=0.01, alpha=c(0.1, 0.05), control.optim=list(), seed)

Arguments

gradient

An object of class data.matrix, or coercible to that class, of dimension n x p, where n is the sample size and p is the dimension of the parameter, containing the contributions of an estimating function

B

An integer indicating the number of outer level bootstrap replications.

M

An integer indicating the number of inner level bootstrap replications. If equal to 0 then single bootstrap is performed.

kB

The proportion of convex hull condition failures in the outer level that the user allows before stopping the algorithm. See ``Details''.

alpha

A vector specifying the desired nominal level(s) of test and confidence set.

control.optim

A list of optional control parameters as passed to the optimisation routine optim for method="L-BFGS-B". See ``Details''.

seed

A single value, interpreted as an integer, recommended to specify seeds.

Value

• A list of class ``Iboot'' containing:
• Call
  The matched call.
• oss
  The observed value of the statistic.
• boot
  The bootstrap distribution of the statistic (sorted into descending order).
• map
  The re-calibrated 1-alpha nominal levels.
• boot.quant
  The 1-alpha level bootstrap quantile(s).
• recalib.quant
  The re-calibrated 1-alpha level bootstrap quantile(s).
• fails.outer
  Actual proportion of failures in the outer level.
• failure
  An error code indicating wheter the algorithm has succeeded: [object Object],[object Object],[object Object]

Details

Iboot performs a single weighted bootstrap iteration in order to provide re-calibrated critical values for tests and confidence sets based on the unstudentised version of the Rao statistic for a parameter $\theta\in R^p,\,p\geq 1$.

Denoted with $y=(y_1,\dots,y_n)$ an i.i.d. sample of size $n$, where $y_i\in R^d,\,d\geq 1$, and with $U(\theta)=\sum_{i=1}^n u(\theta;y_i)$ an (unbiased) estimating function for $\theta$, the unstudentised version of the Rao statistic is defined as $$W_{us}(\theta)=n^{-1} U(\theta)^\top U(\theta).$$

The function Iboot implements the algorithm outlined in Lunardon (2013) that differs from the Nankervis's one for two sided confidence intervals because bootstrap is carried out in a weighted fashion. In the outer level bootstrap versions of $W_{us}(\theta)$, $W^{b}_{us}(\theta),\,b=1,\dots,B$, are obtained by resampling according to the $n$-dimensional vector of weights $w(\theta)=(w_1(\theta),\dots,w_n(\theta))$ associated to each contribution $u(\theta;y_i)$. The functional form of $w_i(\theta)$ is provided by Owen's empirical likelihood formulation, that is $$w_i(\theta)=n^{-1}(1+\lambda(\theta)^\top u(\theta;y_i))^{-1},$$ where $\lambda(\theta)\in R^p$ is the Lagrange multiplier (see Owen, 1990). Instead, in the inner level, bootstrap versions of $W^{b}_{us}(\theta)$ are obtained by resampling according to $w^b(\theta)$, computed as before by using the outer level contributions $u(\theta; y^b_i)$.

As weighted bootstrap might entail computational difficulties, i.e. the convex hull condition (see Owen, 1990) might not be satisfied in all B bootstrap replications, some precautions have been taken. In particular, the algorithm checks if the convex hull condition is satisfied for both the observed contributions $u(\theta;y_i)$ and for each of the B bootstrap resamplings in the outer level, i.e. $u(\theta;y^{b}_i)$.

The algorithm stops immediately if the convex hull condition fails for $u(\theta;y_i)$ whereas concludes after reaching B x kB convex hull condition failures in the outer level. In particular, in the former situation the observed value of the statistic is returned only, whereas in the latter both the observed value and the bootstrap distribution (based on a smaller number of resamplings than B) of the statistic are supplied.

Numerical optimisation of the objective function that supplies $\lambda(\theta)$ relies on the foreign function lbfgsb called from R by optim with method="L-BFGS-B". However, it is not possible to set bounds on the variables by specifying lower and upper, instead the optional parameters that can be set in control.optim are

[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

References

Lunardon, N. (2013). Prepivoting composite score statistics by weighted bootstrap iteration. E-print: arXiv/1301.7026.

Nankervis, J. (2005). Computational algorithms for double bootstrap confidence intervals. Computational statistics & data analysis, 45, 461--475.

Owen, A. (1990). Empirical likelihood ratio confidence regions. The Annals of Statistics, 18, 90--120.

See Also

one.boot, boot, optim

Examples

####Example 1: mean of a normal with known scale
n <- 20
mu <- 1

set.seed(1)
##contributions obtained from the score function
gr <- rnorm(n, mu) - mu

OBJ.Ib <- Iboot(gradient=gr, B=500, M=500, kB=0.01, alpha=c(0.1, 0.05, 0.01), seed=1)

##critical values for testing H0: mu=1, H1: mu!=1
OBJ.Ib
summary(OBJ.Ib)

#####exceeded number of convex hull condition failures in the outer level (kB=0)
OBJ.Ib <- Iboot(gradient=gr, B=500, M=500, kB=0, alpha=c(0.1, 0.05, 0.01), seed=1)
OBJ.Ib
summary(OBJ.Ib)

####Example 2: example 1 of Lunardon (2013)
n <- 20
q <- 10

##parameter
mu <- 0
sig2 <- 1
rho <- 0.5
theta <- c(mu, sig2, rho)

##function to compute the pairwise score contributions
pscore_theta <- function(theta, data) 
{
	n <- nrow(data)
	q <- ncol(data)
	mu <- theta[1]
	sig2 <- theta[2]
	rho <- theta[3]
	A <- matrix(rho, q, q)
	diag(A) <- -(q-1)
	A <- A/((1-rho^2)*sig2^2)
	B <- matrix(-(1+rho^2), q, q)
	diag(B) <- 2*rho*(q-1)
	B <- B/(sig2*(1-rho^2)^2)
	x_dot <- apply(data, 1, sum)
	p_mu <- ((q-1)/(sig2*(1+rho)))*(x_dot-q*mu)
	p_sig2 <- -0.5*apply(((data-mu)	p_rho <- q*(q-1)*rho*0.5/(1-rho^2)-0.5*apply(((data-mu)	RES <- cbind(p_mu, p_sig2, p_rho)
	colnames(RES) <- c("mu", "sig2", "rho")
	RES
}


###data simulation

##correlation matrix
S <- matrix(rho*sig2, q, q)
diag(S) <- sig2
##cholesky
cholS <- chol(S)

##generation from multivariate normal
set.seed(3)
Y <- matrix(rnorm(n*q), n, q)
##compute pairwise score conributions
gr <- pscore_theta(theta, Y)

OBJ.Ib <- Iboot(gradient=gr, B=500, M=500, kB=0.01, alpha=c(0.1, 0.05, 0.01), seed=3)

##critical values for testing H0: theta=(0, 1, 0.5), H1: theta!=(0, 1, 0.5)
OBJ.Ib
summary(OBJ.Ib)

##sampe size too small: convex hull failure at the beginning
n <- 10
set.seed(3)
Y <- matrix(rnorm(n*q), n, q)
##compute pairwise score conributions
gr <- pscore_theta(theta, Y)

OBJ.Ib <- Iboot(gradient=gr, B=500, M=500, kB=0.01, alpha=c(0.1, 0.05, 0.01), seed=3)
OBJ.Ib
summary(OBJ.Ib)