Bootstrapping eigenvalues for nonlinear PCA in r

I am running nonlinear PCA in r, using the homals package. Here is a chunk of the code I am using as an example:

res1 <- homals(data = mydata, rank = 1, ndim = 9, level = "nominal")
res1 <- rescale(res1)

I want to generate 1000 bootstrap estimates of the eigenvalues in this analysis (with replacement), but I can't figure out the code. Does anyone have any suggestions?

Sample data:

dput(head(mydata, 30))

structure(list(`W age` = c(45L, 43L, 42L, 36L, 19L, 38L, 21L, 
27L, 45L, 38L, 42L, 44L, 42L, 38L, 26L, 48L, 39L, 37L, 39L, 26L, 
24L, 46L, 39L, 48L, 40L, 38L, 29L, 24L, 43L, 31L), `W education` = c(1L, 
2L, 3L, 3L, 4L, 2L, 3L, 2L, 1L, 1L, 1L, 4L, 2L, 3L, 2L, 1L, 2L, 
2L, 2L, 3L, 3L, 4L, 4L, 4L, 2L, 4L, 4L, 4L, 1L, 3L), `H education` = c(3L, 
3L, 2L, 3L, 4L, 3L, 3L, 3L, 1L, 3L, 4L, 4L, 4L, 4L, 4L, 1L, 2L, 
2L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 4L), `N children` = c(10L, 
7L, 9L, 8L, 0L, 6L, 1L, 3L, 8L, 2L, 4L, 1L, 1L, 2L, 0L, 7L, 6L, 
8L, 5L, 1L, 0L, 1L, 1L, 5L, 8L, 1L, 0L, 0L, 8L, 2L), `W religion` = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), `W employment` = c(1L, 
1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L), `H occupation` = c(3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 3L, 1L, 1L, 3L, 2L, 4L, 2L, 2L, 
2L, 2L, 4L, 3L, 1L, 1L, 1L, 3L, 1L, 1L, 2L, 2L, 1L), `Standard of living` = 
c(4L, 
4L, 3L, 2L, 3L, 2L, 2L, 4L, 2L, 3L, 3L, 4L, 3L, 3L, 1L, 4L, 4L, 
3L, 1L, 1L, 1L, 4L, 4L, 4L, 3L, 4L, 4L, 2L, 4L, 4L), Media = c(0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Contraceptive = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L)), .Names = c("W age", 
"W education", "H education", "N children", "W religion", "W employment", 
"H occupation", "Standard of living", "Media", "Contraceptive"
), row.names = c(NA, 30L), class = "data.frame")
>

I was given the rescale function to use with the homals package, to do optimal scaling. Here is the function:

rescale <- function(res) {
    # Rescale homals results to proper scaling
    n <- nrow(res$objscores)
    m <- length(res$catscores)
    res$objscores <- (n * m)^0.5 * res$objscores
    res$scoremat <- (n * m)^0.5 * res$scoremat
    res$catscores <- lapply(res$catscores, FUN = function(x) (n * m)^0.5 * x)
    res$cat.centroids <- lapply(res$cat.centroids, FUN = function(x) (n * m)^0.5 * x)
    res$low.rank <- lapply(res$low.rank, FUN = function(x) n^0.5 * x)
    res$loadings <- lapply(res$loadings, FUN = function(x) m^0.5 * x)
    res$discrim <- lapply(res$discrim, FUN = function(x) (n * m)^0.5 * x)
    res$eigenvalues <- n * res$eigenvalues
    return(res)
}

Solution

The standard way to bootstrap in R is to use base package boot.
I am not very satistied with the code that follows because it is throwing lots of warnings. But maybe this is due to the dataset I have tested it with. I have used the dataset and 3rd example in help("homals").

I have run 10 bootstrap replicates only.

library(homals)
library(boot)

boot_eigen <- function(data, indices){
    d <- data[indices, ]
    res <- homals(d, active = c(rep(TRUE, 4), FALSE), sets = list(c(1,3,4),2,5))
    res$eigenvalues
}

data(galo)

set.seed(7578)    # Make the results reproducible
eig <- boot(galo, boot_eigen, R = 10)

eig
#
#ORDINARY NONPARAMETRIC BOOTSTRAP
#
#
#Call:
#boot(data = galo, statistic = boot_eigen, R = 10)
#
#
#Bootstrap Statistics :
#     original      bias    std. error
#t1* 0.1874958  0.03547116 0.005511776
#t2* 0.2210821 -0.02478596 0.005741331

colMeans(eig$t)
#[1] 0.2229669 0.1962961

If this also doesn't run properly in your case, please say so and I will delete the answer.

EDIT.

In order to answer to the discussion in the comments, I have changed the function boot_eigen, the call to homals now follows the question code and rescale is called before returning.

boot_eigen <- function(data, indices){
    d <- data[indices, ]
    res <- homals(data = d, rank = 1, ndim = 9, level = "nominal")
    res <- rescale(res)
    res$eigenvalues
}

set.seed(7578)    # Make the results reproducible
eig <- boot(mydata, boot_eigen, R = 10)