I have microfossil data from across the Cretaceous-Tertiary boundary, 20 species (=variables) and 27 samples (=observations). I have run a PCA with varimax rotation on these data using the rgr package. I now want to determine the significance of the PC loadings using bootstrap BCa confidence intervals from package boot, but I can't get it to work. I have found using Krzanowski cross-validation that 5 PCs appear to be significant (NRETAIN in the script below is thus equal to 5). The script I have been attempting looks like this (the data matrix is called “x”):

BS.vmax <- function(d, f)
{    n <- nrow(d)
     unrot <-gx.mva(d)
         vmax <- gx.rotate(unrot, nrot=NRETAIN)
         vload <- vmax$rload[,1:NRETAIN]

BS.vmax.boot <- boot(x, BS.vmax, R = 99)
boot.ci(BS.vmax.boot, type = "bca")

This script returns a BS.vmax.boot value that looks like this:

boot(data = x, statistic = BS.vmax, R = 99)
Bootstrap Statistics :
          original  bias    std. error
t1*   -0.745596283       0           0
t2*   -0.902044503       0           0
lines removed here
t99*  -0.144587913       0           0
t100*  0.076592734       0           0

which doesn't make any sense (because all biases and SEs are zero)! I will later on increase the number of bootstrap replicates to 2000; this is just for testing.


I don't think the increase in the number of bootstrap replications will help much. The problem is probably due to the small sample size of 27. Even second order accurate bootstrap confidence intervals such as BCa are approximate and the coverage can be way off in small samples. Schenker illustrated this problem with the BC confidence intervals for variance parameters when the population distribution is a highly skewed chi square distribution. Chernick and LaBudde also showed that this would occur for BCa with highly skewed distributions such as lognormal. I have mentioned this previously on this site. Links to the papers are as follows:

Schenker (1) http://www.tandfonline.com/doi/abs/10.1080/01621459.1985.10478123#preview My Bootstrap Methods StatProb article with several references including the Chernick-LaBudde article and our book on the bootstrap (2) http://statprob.com/encyclopedia/ResamplingMethods.html

| cite | improve this answer | |
  • $\begingroup$ Michael, thanks for your response. My post may not have been entirely clear unfortunately, since it was supposed to be concerned with the R script itself and not primarily with the statistics. The dataset provided was just an example, and I want to create an R script that would work in the general case. The problem is that I get a bias of 0 and a std. error also of 0 in the bootstrap statistics for all t1* to t100* as I show in my question, which means that the boot.ci function won't work. How can the script be rewritten to function? $\endgroup$ – Björn Malmgren Jun 22 '12 at 13:40
  • $\begingroup$ Oh I see. I can't help you with that but if the result doesn't make sense at least my answer tells you why yhat might be the case. $\endgroup$ – Michael R. Chernick Jun 22 '12 at 13:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.