# How to actually use the Empirical Influence Function for BCa Bootstrap Intervals?

In the course of seeking reassurance for another part of a hobby analysis* I found a stack answer which mentioned The Jackknife, the Bootstrap and Other Resampling Plans (Efron, 1980). Having managed to find the paper online, the bit about non-parametric bias and skew adjustments to the bounds of bootstrap CIs caught my eye. The problem is that I don't really understand how to do it as described in that paper.

Efron (1980) as I am reading it suggests that the core ingredients of the BCa CI method are the $$U_i$$ and $$a$$ values. The derivation of the latter is simple given you know the $$U_i$$ with formula (pg. 21):

$$a \doteq \frac{1}{6}\frac{\sum\limits_{i=0}^n U_i^3}{\left( \sum\limits_{i=0}^n U_i^2 \right)^{\frac{3}{2}}}$$

The $$U_i$$ have a more complicated formula, coming from what's termed the empirical influence function (also pg. 21):

It is this formula for the $$U_i$$ that I don't know how to work with.

Clicking around in some more stack questions, I found a reference to Bootstrap Methods and their Application (Davison and Hinkley, 1997) but I couldn't find that online. Nevertheless, I did find a slideshow by Davison on ResearchGate based on it which proposes that the $$a$$ value can be jackknifed like so (I've changed the subscript $$j$$ to $$i$$ to match Efron, 1980):

$$l_i \approx l_{jack,i} = (n-1)(\hat{\theta} - \hat{\theta}_i)$$

where $$n$$ is the number of observations in the original sample.

I have tested my understanding of the jackknife estimate of $$a$$ against the $$U_i$$ values provided by Efron (1980) for a $$n=15$$ example (pg. 22) like so:

# first my function for the jackknifing

slide.jack = function(x, y){
theta.hat = cor(x, y)
L = length(x)

theta.i = numeric(L)

for(i in 1:L){
theta.i[i] = cor(x[-i],
y[-i])
}

(L - 1) * (theta.hat - theta.i)
}

# and now for the values as given in Efron (1980):

lsat = c(576, 635, 558, 578, 666, 580, 555, 661, 651, 605, 653, 575,
545, 572, 594)

gpa = c(3.39, 3.3, 2.81, 3.03, 3.44, 3.07, 3, 3.43, 3.36, 3.13, 3.12,
2.74, 2.76, 2.88, 2.96)

given.ui = c(-1.507, 0.168, 0.273, 0.004, 0.525, -0.049, -.1, 0.477,
0.31, 0.004, -0.526, -0.091, 0.434, 0.125, -0.048)


The values given by slide.jack do indeed generally approximate the values in given.ui so I can just use this process, but before I found the slideshow I searched for "empirical influence function R" in Google and was directed to empinf in package boot. However, I couldn't figure out how to obtain values close to given.ui using empinf:

wcor = function(data, w){
cor(w * data$$x, w * data$$y)
}
# a weights based function

empinf(data = data.frame(x = lsat, y = gpa, w = 1/15),
statistic = wcor, stype = "w", type = "inf")


While I would like to be able to use empinf on this data (presumably the whack answers I got are reflective of human error on my end), my main focus is learning how to put the& to work.

Incidentally, I achieve similarly good approximations with both of:

eco.paper = function(x, y){
L = length(x)
store = numeric(L)
for(i in 1:L){
store[i] = (cor(x, y) - cor(x[-i], y[-i])) / (1 / (L - 1))
}
store
}

# L-1 worked better than just L
# based on a formula in "Asymptotic and Bootstrap Inference for
# Inequality and Poverty Measures"
# (Davidson and Flachaire, 2005), pg. 8
# the division by (1 / (L-1)) was my own innovation; Davidson and
# Flachaire's formula just looked like it was a common ratio away from
# the actual U_i

# I present now a failed attempt to understand a different a formula
# for a that produced a value similar to 1 / (L-1) for a ... totally
# wrong for a but potentially useful for calculating the UI

accelerator.corr = function(x, y){
L = length(x)
theta.i = numeric(L)

for(i in 1:L){
theta.i[i] = cor(x[-i], y[-i])
}

theta.dot = mean(theta.i)

list(a = (1 / 6) * ( sum((theta.i - theta.dot)^3) /
(sum((theta.i - theta.dot)^2))^1.5),
jacknife_U_i = theta.i)
}

# and now the second stage

twostage = function(x, y){
L = length(x)
store = numeric(L)
for(i in 1:L){
store[i] = (cor(x, y) - cor(x[-i], y[-i])) /
accelerator.corr(x, y)\$a
}
store
}

(comps = data.frame(given.ui,
eco = eco.paper(lsat, gpa),
jacknife = slide.jack(lsat, gpa),
twostage = twostage(lsat, gpa)))
sum((comps[, 1] - comps[, 2])^2)
sum((comps[, 1] - comps[, 3])^2)
sum((comps[, 1] - comps[, 4])^2)


The thing is I have no idea whether eco or twostage actually generally work, or if they just happened to for the specific $$U_i$$ in this example.

*I have asked ChatGPT to produce "a value from 0 to 100, which represents your guess of what the average person would rate the film from 0-100" for 330 movies. Collecting this data is extremely tedious and after many days of mindlessly copying and pasting I have N = 30 guesses for each film. My concern is whether the values are "sensible" or "hallucinations". I believe that if they are sensible, the distribution of guesses ought to be symmetric (unless boundary distortions occur).

I have managed to figure out the answers to some of my questions but I am still at a loss with how to work with the empirical influence function itself. Luckily, I think I can now just rely on empinf. I'd still like to know, though.

# Regarding my problems with empinf

I have downloaded so many pdfs in the last few days and in re-reading the Package 'boot' pdf, I discovered an example function:

abc.ci(bigcity, corr, conf=c(0.90,0.95))


corr is what I often think cor is so this piqued my interest immediately. I plugged it in like so:

(nw.ui = empinf(data = data.frame(x = lsat, y = gpa),
statistic = corr, stype = "w", type = "inf"))


and this yielded $$U_i$$ values that very closely approximate the given.ui I copied from Efron. The problem was simply that in trying several schemes to weight cor I was doing several stupid things. What I needed was something more like:

rewrite.weight = function(d, w){
num = sum(d$$x * d$$y * w) - ((sum(d$$x * w)) * sum(d$$y * w))
denom = sqrt(sum(d$$x^2 * w) - (sum(d$$x *w))^2) *
sqrt(sum(d$$y^2 * w) - (sum(d$$y *w))^2)

num / denom
}


Basically how I was trying to incorporate the weights before was mathematically invalid.

It also turned out that I just needed to create a variable for the weights... it didn't need anything in it.

# Regarding my Approximating Functions

Again, I was making stupid maths mistakes -- in my defence, here I was very tired -- because eco.paper and slide.jack are just rearranged functions with different variable names. The key bits:

• theta.hat is the same as cor(x, y),
• theta.i = cor(x[-i], y[-i]) is literally just cor(x[-i], y[-i]) but stored as a variable, and
• (L - 1) * is equivalent to / (1 / (L - 1))

twostage would appear to have worked entirely coincidentally.