R/fa.parallel.poly.R

In psych: Procedures for Psychological, Psychometric, and Personality Research

 #parallel analysis of polychoric factor analysis
 "fa.parallel.poly" <- 
function(x,n.iter=10,SMC=TRUE,fm="minres",correct=TRUE,sim=FALSE,fa="both",global=TRUE) { 
p <- .05
 cl <- match.call()
.Deprecated("fa.parallel.poly", msg = "fa.parallel.poly is deprecated.  Please use the fa.parallel function with the cor='poly' option.")


n.obs <- dim(x)[1]
	tx <- table(as.matrix(x))
	nx <- dim(tx)[1] 
	if(dim(tx)[1] ==2) {tet <- tetrachoric(x,correct=correct)
	    typ = "tet"
	    if(sim) { tx.item <- matrix(apply(x,2,table),ncol=ncol(x))
	    px <- matrix(tx.item/colSums(tx.item),ncol=2,byrow=TRUE)
	    } } else {tet <- mixedCor(x,global=global)
	    typ = "poly"}
	    cat("\n")  #to try to clear the progress bar
        flush(stdout())
 r <- tet$rho
 f <- fa(r,n.obs=n.obs,SMC = SMC,covar=FALSE,fm=fm) #call fa with the appropriate parameters
 f$Call <- cl
 fl <- f$loadings  #this is the original
 nvar <- dim(fl)[1]

 e.values <- list(pc =list(),fa = list())
 replicates <- list()
 rep.rots <- list()
 for (trials in 1:n.iter) {
progressBar(trials,n.iter,"fa.parallel.poly")
 
 bad <- TRUE
 while(bad) {
 xs <- matrix(apply(x,2,function(y) sample(y,n.obs,replace=TRUE)),ncol=nvar) #do it column wise

tets <- polychoric(xs,progress=FALSE,global=global)
  r <- tets$rho
  bad <- any(is.na(r))
  }
  values.samp <- eigen(tets$rho)$values
  e.values[["pc"]][[trials]] <- values.samp
   					
   				
  fs <- fa(r,n.obs=n.obs,SMC = SMC,covar=FALSE,fm=fm) #call fa with the appropriate parameters
     e.values[["fa"]][[trials]] <- fs$values
     
if(sim && (typ=="tet")) {xsim <- matrix(apply(px,1,function(x) sample(2,n.obs,TRUE,prob=x)),ncol=nvar)
    rho.sim <- tetrachoric(xsim,correct=correct)
    values.sim <- eigen(rho.sim$rho)$values
    e.values[["pc.sim"]][[trials]] <- values.sim
    fsim <- fa(rho.sim$rho,n.obs=n.obs,SMC = SMC,covar=FALSE,fm=fm) #call fa with the appropriate parameters
     e.values[["fa.sim"]][[trials]] <- fsim$values

}
     
} 

cat("\n")  #to try to clear the progress bar

ei.pc <-describe(matrix(unlist(e.values$pc),ncol=nvar,byrow=TRUE))  #eigen values of pcs
ei.fa <- describe(matrix(unlist(e.values$fa),ncol=nvar,byrow=TRUE)) #eigen values of fa
fa.test <- which(!(f$values > ei.fa$mean))[1]-1
pc.test <- which(!(f$e.values > ei.pc$mean))[1] -1
#e.stats <- list(fa.values=f$values,pc.values=f$e.values,pc.sim=ei.pc,fa=ei.fa,nfact=fa.test,ncomp = pc.test)
if(sim && (typ=="tet")) {eis.pc <- describe(matrix(unlist(e.values$pc.sim),ncol=nvar,byrow=TRUE))  #eigen values of pcs                
       eis.fa <- describe(matrix(unlist(e.values$fa.sim),ncol=nvar,byrow=TRUE)) #eigen values of fa
} else {eis.pc <- NULL
        eis.fa <- NULL
}

#results <- list(fa = f,rho=tet$rho,tau=tet$tau,n.obs=n.obs,Call= cl,e.values=e.values,e.stats=e.stats)
results <- list(rho=tet$rho,tau=tet$tau,n.obs=n.obs,Call= cl,fa.values=f$values,pc.values=f$e.values,pc.sim=ei.pc,fa.sim=ei.fa,pcs.sim=eis.pc,fas.sim=eis.fa,nfact=fa.test,ncomp = pc.test)
class(results) <- c("psych","parallel")
cat('\n See the graphic output for a description of the results\n')
plot.poly.parallel(results,fa=fa)
return(results)
 }
 #written May 8 2011
 #corrected December 27, 2011 to pass the fm correctly (had always forced minres)  
 #modified Sept 16, 2013 to use mixed.cor instead of polychoric  (more robust) 
 #modified Oct 2, 2013 to add the ability to do random data as well 
 #modified Oct 22, 2013 to allow choice of what to plot
 #modified 1/15/14 to pass the global parameter
 #modified 3/24/14 to check for bad resamples
 "plot.poly.parallel" <-
 function(x,show.legend=TRUE,fa="both",...) {
 e.values <- x$pc.values
 values <-  x$fa.values
pcm <- x$pc.sim$mean
fam <- x$fa.sim$mean

switch(fa,

both = {  
 plot(e.values,type="b", main = "Eigen values of tetrachoric/polychoric matrix",ylab="Eigen values of original  and simulated factors and components",ylim=c(0,max(e.values)) ,xlab="Factor Number",pch=4,col="blue")
 points(values,type ="b",pch=2,col="blue")
 points(pcm,type="l",lty="dotted",pch=2,col="red")
 points(fam,type="l",lty="dashed",pch=2,col="red")
 if(!is.null(x$pcs.sim)) points(x$pcs.sim$mean,type="l",lty="dashed",pch=2,col="red")
 if(!is.null(x$fas.sim)) points(x$fas.sim$mean,type="l",lty="dashed",pch=2,col="red")
    },
    
 fa =  {plot(values,type="b", main = "Eigen values of tetrachoric/polychoric matrix",ylab="Eigen values of original  and simulated factors ",ylim=c(0,max(e.values)) ,xlab="Factor Number",pch=2,col="blue")
        points(fam,type="l",lty="dashed",pch=2,col="red")},
 
 pc = {plot(e.values,type="b", main = "Eigen values of tetrachoric/polychoric matrix",ylab="Eigen values of original  and simulated factors and components",ylim=c(0,max(e.values)) ,xlab="Factor Number",pch=4,col="blue")

 points(pcm,type="l",lty="dotted",pch=2,col="red")} )
 
if(show.legend) {
switch(fa, 
both ={
 if(!is.null(x$pcs.sim)) {
       legend("topright", c("PC  Actual Data", " PC  Resampled Data"," PC  Simulated Data","FA  Actual Data", " FA  Resampled Data", " FA  Simulated Data"), col = c("blue","red","red","blue","red","red"),pch=c(4,NA,NA,2,NA,NA),
       text.col = "green4", lty = c("solid","dotted","dashed","solid","dotted","dashed"),
       merge = TRUE, bg = 'gray90') 
       } else {legend("topright", c("PC  Actual Data", " PC  Resampled Data","FA  Actual Data", " FA  Resampled Data"), col = c("blue","red","blue","red"),pch=c(4,NA,2,NA),
       text.col = "green4", lty = c("solid","dotted","solid","dotted"),
       merge = TRUE, bg = 'gray90') 
 } 
   },

fa ={ legend("topright", c("FA  Actual Data", " FA  Resampled Data"), col = c("blue","red"),pch=c(2,NA),
       text.col = "green4", lty = c("solid","dotted"),
       merge = TRUE, bg = 'gray90')  },
pc =  { legend("topright", c("PC  Actual Data", " PC  Resampled Data"), col = c("blue","red"),pch=c(4,NA),
       text.col = "green4", lty = c("solid","dotted"),
       merge = TRUE, bg = 'gray90')  } )
 }  
fa.test <- which(!(values > fam))[1]-1
pc.test <- which(!(e.values > pcm))[1] -1
cat("Parallel analysis suggests that ")
cat("the number of factors = ",fa.test, " and the number of components = ",pc.test,"\n")
 }
 #modified Oct 2, 2013 to plot (if desired) the simulated as well as the resampled data
 #modified Oct 22, 2013 to allow the plot to choose beteween both, fa, and pc.