Bootstrapping eigenvalues for nonlinear PCA in r

Question

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)
}

Welcome to SO! You can bootstrap by resampling `mydata` with replacement by using `sample`, then re-running your code. If you need more help than that, you'll need to provide a reproducible example. — csgroen, Apr 12 '18 at 13:03

Rui Barradas · Accepted Answer · 2018-04-14T15:49:30.920

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)

Thank you so much for your help! This seems to work, except I have 10 variables in my dataset, and because of that it doesn't work. Could you help me change the script so that it runs with 10 variables? — Eszter, Apr 14 '18 at 12:05
@Eszter You should edit the question with the output of `dput(head(mydata, 30))` for us to have sample data. Like this it is difficult to do much more. — Rui Barradas, Apr 14 '18 at 12:41
@Eszter Thanks, that is exactly what I meant. Also, function `rescale` comes from what package? — Rui Barradas, Apr 14 '18 at 12:58
I've updated the question with the function now included. I got this function from my professor. Sorry if I'm making this difficult — Eszter, Apr 14 '18 at 13:09

Bootstrapping eigenvalues for nonlinear PCA in r

1 Answers1