I can't find this type of bootstrap problem described anywhere, where the data is so unbalanced: 50 y-values vs 2500 values yielding 50 means that become the x's. Another thing: the 50 strata yielding the means are internally autocorrelated time series. The bootstrapped mean histograms looked mostly normal distributed.
Problem 1: I just can't intuit if the improvement in confidence interval is too good to be true when I bootstrap from all 2500 values instead of directly from the 50 x, y-pairs: A span of 0.39 improved to 0.039. I.e. 10 times narrower. It makes no sense when having only 50 y-values available - or does it? I started to doubt either the bootstrap in this context, or my own R-code provided below, or my theoretical understanding. Should I trust the result or is there a theoretical reason not to? (I have checked out a few papers and I know bootstrap is not always that precise and that better variants and methods exist for correlation CIs, e.g: Confidence intervals for correlations when data are not normal)
bootstrap directly on means: R = 0.608 ( 0.3878, 0.7790 )
bootstrap from underlying values: R = 0.608 ( 0.5916, 0.6290 )
Problem 2: Why is including all data better according to this experiment? I can reason it both ways:
(Wrong way) Say I bootstrap from 50 fix y-values, and 50 X's that are averages. Then I have fix x-values with no variation and should get a narrower CI than if I introduce an uncertainty - a spread in each x average - by using all 2500 values.
(correct by experiment) On the other hand: Say I use all 2500 values in a two-step bootstrap computing first the mean and from it the x,y-correlation. Then I use more information than is inherent in the 50 x-averages, and should get a narrower confidence interval for the Pearson correlation.
I expected the second method to yield narrower confidence interval for Pearson's R since more information is utilized, but I don't know how or why and I am curious.
This is not commercial in any way. I help out a friend that is about to publish a paper I believe could be beneficial to general health. I benefit from it myself mostly by learning R (yes I am an R and Tidyverse newbie), and by the thrill of doing some real science for a good cause. I am a system developer by profession.
Here is the R code for using all data, with the unnecessaries removed:
testst <-function(tb,i) {
# Generate correlation between yStrat and xStrat.
# arg: tibble with about 2500 rows with about 50 strata
tbi <- tb[i, ] # Bootstrap sample
yTbi <- tbi %>% group_by(StratumID) %>%
summarise(yStrat=last(YPrevalence)) # All Y identical per stratum
# Could have used unique as well
# Apply mean to each stratum of bootstrapped draws from stratum
xTbi <- tbi %>% group_by(StratumID) %>%
summarise(xStrat = mean(XConsumption, na.rm = TRUE))
# Apply cor to fixed y values against each stratum mean
pears <- cor(x = xTbi$xStrat,
y = yTbi$yStrat,
method = "pearson")
return(c(pears))
}
# Executing on tibble tbd with ~2500 rows, ~50 strata. One fix Y-value per stratum
tbd$StratumID <- as_factor(tbd$StratumID)
bt <- boot(tbd,testst,R=10000,strata=tbd$StratumID)
bc<-boot.ci(bt,type="bca",index=1)
