#corrected estimate of communality, May 21, 2007 #removed "x" from factanal call, June 14, 2007 #added ability to do 2 factors by treating them with equal loadings Jan 2008 #added use of simplimax rotation June 2008 #corrected the sign of group factors to match the original factors "schmid" <- function (model, nfactors = 3, fm = "minres", digits=2,rotate="oblimin",n.obs=NA,option="equal",Phi=NULL,...) { #if Phi is not Null, then we have been given a factor matrix, otherwise #model is a correlation matrix, or if not, the correlation matrix is found #nfactors is the number of factors to extract if(!require(GPArotation)) {stop("I am sorry, you need to have the GPArotation package installed")} if(is.null(Phi)) { #the normal case normal.case <- TRUE nvar <-dim(model)[2] if(dim(model)[1] != dim(model)[2]) {n.obs <- dim(model)[1] model <- cor(model,use="pairwise")} if (fm =="pc") { fact <- principal(model, nfactors,n.obs=n.obs,...) } else {if ((fm == "pa") |(fm =="minres") | (fm =="wls") |(fm =="minres") |(fm =="ml") |(fm =="gls")) {fact <- fa(model, nfactors,n.obs=n.obs,rotate="varimax",fm=fm,...) } else { #fact <- factanal(covmat = model, factors = nfactors,n.obs=n.obs,...) stop("The method of factor extraction you specified is not avaialble") }} orth.load <- loadings(fact) } else {model <- as.matrix(model) Phi <- as.matrix(Phi) fact <- model %*% Phi #find the orthogonal matrix from the oblique pattern and the Phi matrix orth.load <- fact nfactors <- dim(fact)[2] normal.case <-FALSE} colnames(orth.load) <- paste("F",1:nfactors,sep="") if(nfactors == 1) { message("Omega_h for 1 factor is not meaningful, just omega_t") obminfact <-list(loadings= orth.load) factr <- 1 } else { #the normal case is nfactors > 2 if (rotate == "simplimax") {obminfact <- simplimax(orth.load)} else { if((rotate == "promax") | (rotate == "Promax") ) {obminfact <- Promax(orth.load) rotmat <- obminfact\$rotmat Phi <- obminfact\$Phi } else { if ((rotate == "cluster") | (rotate == "target")) {obminfact <- varimax(orth.load) obminfact <- target.rot(obminfact,...) loadings <- obminfact\$loadings Phi <- obminfact\$Phi} else { obminfact <- try(oblimin(orth.load)) if(class(obminfact)== as.character("try-error")) {obminfact <- Promax(orth.load) #special case for examples with exactly 2 orthogonal factors message("\nThe oblimin solution failed, Promax used instead.\n") #perhaps no longer necessary with patch to GPForth and GPFoblq in GPArotation rotmat <- obminfact\$rotmat Phi <- obminfact\$Phi} }} } } if(nfactors > 1) rownames(obminfact\$loadings) <- attr(model,"dimnames")[[1]] if(!normal.case) { fload <- model factr <- Phi} else { fload <- obminfact\$loadings #factr <- t(obminfact\$Th) %*% (obminfact\$Th) factr <- obminfact\$Phi} if (nfactors ==1) {gload <- c(1) warning("Omega_h and Omega_assymptotic are not meaningful with one factor") } else { colnames(factr) <- rownames(factr) <- paste("F",1:nfactors,sep="") #make it a vector if (nfactors>2) { gfactor <- fa(factr) #The first factor of the factor intercorrelation matrix gload <- loadings(gfactor) } else {gload<- c(NA,NA) #consider the case of two factors if(option=="equal") { gload[1] <- sqrt(abs(factr[1,2])) gload[2] <- sign(factr[1,2])*sqrt(abs(factr[1,2])) message("Three factors are required for identification -- general factor loadings set to be equal. Proceed with caution.")} else { if(option=="first") { gload[1] <- 1 # gload[2] <- abs(factr[1,2]) gload[2] <- (factr[1,2]) message("Three factors are required for identification -- general factor loading set to be 1 for group factor 1. Proceed with caution.")} else { gload[2] <- 1 # gload[1] <- abs(factr[1,2]) gload[1] <- (factr[1,2]) message("Three factors are required for identification -- general factor loadings are set to be 1 for group factor 2. Proceed with caution.")} } } } gprimaryload <- fload %*% gload colnames(gprimaryload) <- "g" u2 <- 1 - diag(orth.load %*% t(orth.load)) h2 <- 1 - u2 uniq <- 1 - fload^2 guniq <- 1 - gprimaryload^2 Ig <- matrix(0, ncol = nfactors, nrow = nfactors) diag(Ig) <- gload primeload <- fload %*% Ig g.percent <- gprimaryload^2/h2 colnames(g.percent) <- "p2" uniq2 <- 1 - uniq - primeload^2 uniq2[uniq2<0] <- 0 sm <- sign(fload) * sqrt(uniq2) #added June 1, 2010 to correctly identify sign of group factors colnames(sm) <- paste("F",1:nfactors,"*",sep="") if(!is.null(Phi)) { result <- list(sl = cbind(gprimaryload, sm,h2, u2,p =g.percent), orthog = orth.load, oblique=fload, phi =factr, gloading = gload)} else{ result <- list(sl = cbind(gprimaryload, sm,h2, u2,p=g.percent), orthog = orth.load, oblique=fload, phi =factr, gloading = gload,dof=fact\$dof,objective=fact\$criteria[1],STATISTIC=fact\$STATISTIC,PVAL=fact\$PVAL,RMSEA=fact\$RMSEA,BIC=fact\$BIC,rms = fact\$rms,crms=fact\$crms,n.obs=n.obs )} # class(result) <- c("psych" ,"fa") return(result) }