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).
