# Bootstrap PCA with varimax rotation

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)
unrot$eigenvectors 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:

ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
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.

-
add comment

## 1 Answer

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

-
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? –  Björn Malmgren Jun 22 '12 at 13:40
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. –  Michael Chernick Jun 22 '12 at 13:50
add comment