Confidence intervals of loadings in principal components in R

Question

I am using following code for principal component analysis of first 4 columns of iris data set using prcomp function in R:

> prcomp(iris[1:4])
Standard deviations:
[1] 2.0562689 0.4926162 0.2796596 0.1543862

Rotation:
                     PC1         PC2         PC3        PC4
Sepal.Length  0.36138659 -0.65658877  0.58202985  0.3154872
Sepal.Width  -0.08452251 -0.73016143 -0.59791083 -0.3197231
Petal.Length  0.85667061  0.17337266 -0.07623608 -0.4798390
Petal.Width   0.35828920  0.07548102 -0.54583143  0.7536574

How can I get confidence intervals of these values in R? Is there any package that can do it? Thanks for your help.

Answer 1

You could use bootstrapping on this. Simply re-sample your data with the bootstrapping package and record the principal components computed every time. Use the resulting empirical distribution to get your confidence intervals.

The boot package makes this pretty easy.

Here is an example calculating the Confidence Interval at 95% for the first PCA component with respect to Sepal.Length:

library(boot)

getPrcStat <- function (samdf,vname,pcnum){
  prcs <- prcomp(samdf[1:4]) # returns matrix
  return(prcs$rotation[ vname,pcnum ])   # pick out the thing we need
}

bootEst <- function(df,d){
   sampledDf <- df[ d, ]  # resample dataframe 
   return(getPrcStat(sampledDf,"Sepal.Length",1))
}

bootOut <- boot(iris,bootEst,R=10000)
boot.ci(bootOut,type=c("basic"))

The output is:

  BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
  Based on 10000 bootstrap replicates

  CALL : 
  boot.ci(boot.out = bootOut, type = c("basic"))

  Intervals : 
  Level      Basic         
  95%   ( 0.3364,  1.1086 )  
  Calculations and Intervals on Original Scale

So using the usual basic bootstrap method we get a 95 percent confidence interval of between 0.3364 and 1.1086. There are plenty of other more advanced statistical methods that can be used too, but you need to know what you are doing.