# # Bootstrapping Eigenvalues and Screeplots (Second attempt): # # We use a simple percentile bootstrap and calculate 95% CI # for the eigenvalues and the scree plot. # # Create new function from existing package "psy"; # We are modifying function "scree.plot": scree.plot.boot.CI<-function (namefile, title = "Scree Plot", type = "R", use = "complete.obs", simu = "F") { mat <- namefile if (use == "complete.obs") mat <- na.omit(namefile) if (type == "R") eigenval <- eigen(cor(mat, use = "pairwise.complete.obs"), symmetric = TRUE)$values if (type == "V") eigenval <- eigen(cov(mat, use = "pairwise.complete.obs"), symmetric = TRUE)$values if (type == "E") eigenval <- namefile if (type == "M") eigenval <- eigen(namefile, symmetric = TRUE)$values nev <- length(eigenval) plot(eigenval, type = "b", pch = 16, bty = "o", main = title, xlab = "Dimension", ylab = "Eigenvalue") lines(c(1, nev), c(1, 1), lty = 2) if (is.numeric(simu) && (type == "R")) { n <- dim(mat)[1] p <- dim(mat)[2] matsimu <- matrix(nrow = n, ncol = p) int <- rep(1, n * p) attr(int, "dim") <- c(n, p) mat <- pmax(as.matrix(mat), int) eigen.temp.matrix<-matrix(rep(0,simu*p), ncol=p) for (i in 1:simu) { matnorm<-mat[sample(1:n, n, replace=T), ] attr(matnorm, "dim") <- c(n, p) matsimu <- (mat/mat) * matnorm eigenval <- eigen(cor(matsimu, use = "pairwise.complete.obs"))$values eigen.temp.matrix[i,]<-eigenval points(eigenval, type = "l") } } eigen.temp.matrix eigen.boot.mean<-apply(eigen.temp.matrix, 2, mean) eigen.boot.sd<-apply(eigen.temp.matrix, 2, sd) eigen.boot.quantile<-apply(eigen.temp.matrix, 2, quantile, probs=c(.025, .975)) return(list(boot.eigen.values=eigen.temp.matrix, bootstrap.means.eigenvalues=eigen.boot.mean, std.error.eigenvalues=eigen.boot.sd, CI.95.percent.eigenvalues=eigen.boot.quantile)) } # Load library to use data set from packages "psy" library(psy) data(expsy) expsy ##### Running the bootstrap scree plot and getting it's output: # # In the "$boot.eigen.values matrix we have a row representing a set of # bootstrapped eigenvalues. There are "simu" rows of eigenvalues # and "p" columns correpsonding to p eigenvalues for each bootstrap sample. # # In the $bootstrap.means.eigenvalues matrix we have the average boostrapped eigenvalues # for the pth variable # # In the $std.error.eigenvalues matrix we have the standard deviation of the # bootstrap eigenvalues for the pth variable # # In the $CI.95.percent.eigenvalues matrix we have the 2.5th and 97.5th # percentiles of the boostrapped eigenvalues - one for each variable. scree.plot.boot.CI(expsy[,1:10],simu=20,use="P")