1"factor.wls" <- 2function(r,nfactors=1,residuals=FALSE,rotate="varimax",n.obs = NA,scores=FALSE,SMC=TRUE,missing=FALSE,impute="median", min.err = .001,digits=2,max.iter=50,symmetric=TRUE,warnings=TRUE,fm="wls") { 3 cl <- match.call() 4 .Deprecated("fa",msg="factor.wls is deprecated. Please use the fa function with fm=wls.") 5 #this does the WLS or ULS fitting depending upon fm 6 "fit.residuals" <- function(Psi,S,nf,S.inv,fm) { 7 diag(S) <- 1- Psi 8 if(fm=="wls") {if(!is.null(S.inv)) sd.inv <- diag(1/diag(S.inv)) } else {if(!is.null(S.inv)) sd.inv <- ((S.inv))} #gls 9 10 eigens <- eigen(S) 11 eigens$values[eigens$values < .Machine$double.eps] <- 100 * .Machine$double.eps 12 13 if(nf >1 ) {loadings <- eigens$vectors[,1:nf] %*% diag(sqrt(eigens$values[1:nf])) } else {loadings <- eigens$vectors[,1] * sqrt(eigens$values[1] ) } 14 model <- loadings %*% t(loadings) 15 #weighted least squares weights by the importance of each variable 16 if(fm=="wls" ) {residual <- sd.inv %*% (S- model)^2 %*% sd.inv} else {if(fm=="gls") {residual <- (sd.inv %*% (S- model))^2 } else {residual <- (S-model)^2 } } # the uls solution usually seems better? 17 diag(residual) <- 0 18 error <- sum(residual) 19 } 20 21 #this code is taken (with minor modification to make ULS or WLS) from factanal 22 #it does the iterative calls to fit.residuals 23 "fit" <- function(S,nf,fm) { 24 S.smc <- smc(S) 25 if((fm=="wls") | (fm =="gls") ) {S.inv <- solve(S)} else {S.inv <- NULL} 26 if(sum(S.smc) == nf) {start <- rep(.5,nf)} else {start <- 1- S.smc} 27 res <- optim(start, fit.residuals, method = "L-BFGS-B", lower = .005, 28 upper = 1, control = c(list(fnscale = 1, parscale = rep(0.01, 29 length(start)))), nf= nf, S=S, S.inv=S.inv,fm=fm ) 30 31 if((fm=="wls") | (fm=="gls")) {Lambda <- FAout.wls(res$par, S, nf)} else { Lambda <- FAout(res$par, S, nf)} 32 result <- list(loadings=Lambda,res=res) 33 } 34 35 #these were also taken from factanal 36 FAout <- function(Psi, S, q) { 37 sc <- diag(1/sqrt(Psi)) 38 Sstar <- sc %*% S %*% sc 39 E <- eigen(Sstar, symmetric = TRUE) 40 L <- E$vectors[, 1L:q, drop = FALSE] 41 load <- L %*% diag(sqrt(pmax(E$values[1L:q] - 1, 0)), 42 q) 43 diag(sqrt(Psi)) %*% load 44 } 45 46 FAout.wls <- function(Psi, S, q) { 47 diag(S) <- 1- Psi 48 E <- eigen(S,symmetric = TRUE) 49 L <- E$vectors[,1:q,drop=FALSE] %*% diag(sqrt(E$values[1:q,drop=FALSE]),q) 50 return(L) 51 } ## now start the main function 52 53 54 ## now start the main function 55 56 if((fm !="pa") & (fm !="minres" )& (fm != "gls") & (fm != "wls")) {message("factor method not specified correctly, weighted least squares used") 57 fm <- "wls" } 58 59 n <- dim(r)[2] 60 if (n!=dim(r)[1]) { n.obs <- dim(r)[1] 61 if(scores) {x.matrix <- r 62 if(missing) { #impute values 63 miss <- which(is.na(x.matrix),arr.ind=TRUE) 64 if(impute=="mean") { 65 item.means <- colMeans(x.matrix,na.rm=TRUE) #replace missing values with means 66 x.matrix[miss]<- item.means[miss[,2]]} else { 67 item.med <- apply(x.matrix,2,median,na.rm=TRUE) #replace missing with medians 68 x.matrix[miss]<- item.med[miss[,2]]} 69 }} 70 r <- cor(r,use="pairwise") # if given a rectangular matrix, then find the correlations first 71 } else { 72 if(!is.matrix(r)) { r <- as.matrix(r)} 73 sds <- sqrt(diag(r)) #convert covariance matrices to correlation matrices 74 r <- r/(sds %o% sds) } #added June 9, 2008 75 if (!residuals) { result <- list(values=c(rep(0,n)),rotation=rotate,n.obs=n.obs,communality=c(rep(0,n)),loadings=matrix(rep(0,n*n),ncol=n),fit=0)} else { result <- list(values=c(rep(0,n)),rotation=rotate,n.obs=n.obs,communality=c(rep(0,n)),loadings=matrix(rep(0,n*n),ncol=n),residual=matrix(rep(0,n*n),ncol=n),fit=0)} 76 77 78 79 r.mat <- r 80 Phi <- NULL 81 colnames(r.mat) <- rownames(r.mat) <- colnames(r) 82 if(SMC) { 83 if(nfactors < n/2) {diag(r.mat) <- smc(r) } else {if (warnings) message("too many factors requested for this number of variables to use SMC, 1s used instead")} } 84 orig <- diag(r) 85 86 87 comm <- sum(diag(r.mat)) 88 err <- comm 89 i <- 1 90 comm.list <- list() 91 if(fm=="pa") { 92 while(err > min.err) #iteratively replace the diagonal with our revised communality estimate 93 { 94 eigens <- eigen(r.mat,symmetric=symmetric) 95 #loadings <- eigen.loadings(eigens)[,1:nfactors] 96 if(nfactors >1 ) {loadings <- eigens$vectors[,1:nfactors] %*% diag(sqrt(eigens$values[1:nfactors])) } else {loadings <- eigens$vectors[,1] * sqrt(eigens$values[1] ) } 97 model <- loadings %*% t(loadings) 98 99 new <- diag(model) 100 101 comm1 <- sum(new) 102 diag(r.mat) <- new 103 err <- abs(comm-comm1) 104 if(is.na(err)) {warning("imaginary eigen value condition encountered in factor.pa,\n Try again with SMC=FALSE \n exiting factor.pa") 105 break} 106 comm <- comm1 107 comm.list[[i]] <- comm1 108 i <- i + 1 109 if(i > max.iter) {if(warnings) {message("maximum iteration exceeded")} 110 err <-0 } 111 } 112 113 } 114 115 if((fm == "wls")| (fm=="gls")) { #added May 25, 2009 to do WLS fits 116 uls <- fit(r,nfactors,fm) 117 eigens <- eigen(r) #used for the summary stats 118 result$par <- uls$res 119 loadings <- uls$loadings 120 } 121 122 # a weird condition that happens with the Eysenck data 123 #making the matrix symmetric solves this problem 124 if(!is.real(loadings)) {warning('the matrix has produced imaginary results -- proceed with caution') 125 loadings <- matrix(as.real(loadings),ncol=nfactors) } 126 #make each vector signed so that the maximum loading is positive - probably should do after rotation 127 #Alternatively, flip to make the colSums of loading positive 128 if (FALSE) { 129 if (nfactors >1) { 130 maxabs <- apply(apply(loadings,2,abs),2,which.max) 131 sign.max <- vector(mode="numeric",length=nfactors) 132 for (i in 1: nfactors) {sign.max[i] <- sign(loadings[maxabs[i],i])} 133 loadings <- loadings %*% diag(sign.max) 134 135 } else { 136 mini <- min(loadings) 137 maxi <- max(loadings) 138 if (abs(mini) > maxi) {loadings <- -loadings } 139 loadings <- as.matrix(loadings) 140 if(fm == "minres") {colnames(loadings) <- "MR1"} else {colnames(loadings) <- "PA1"} 141 } #sign of largest loading is positive 142 } 143 144 #added January 5, 2009 to flip based upon colSums of loadings 145 if (nfactors >1) {sign.tot <- vector(mode="numeric",length=nfactors) 146 sign.tot <- sign(colSums(loadings)) 147 loadings <- loadings %*% diag(sign.tot) 148 } else { if (sum(loadings) <0) {loadings <- -as.matrix(loadings)} else {loadings <- as.matrix(loadings)} 149 colnames(loadings) <- "MR1" } 150 151 152 #end addition 153 if(fm == "wls") {colnames(loadings) <- paste("WLS",1:nfactors,sep='') } else {if(fm == "gls") {colnames(loadings) <- paste("GLS",1:nfactors,sep='')} else {colnames(loadings) <- paste("MR",1:nfactors,sep='')}} 154 rownames(loadings) <- rownames(r) 155 loadings[loadings==0.0] <- 10^-15 #added to stop a problem with varimax if loadings are exactly 0 156 157 model <- loadings %*% t(loadings) 158 159 f.loadings <- loadings #used to pass them to factor.stats 160 if(rotate != "none") {if (nfactors > 1) { 161 162 163 if (rotate=="varimax" | rotate=="quartimax") { 164 rotated <- do.call(rotate,list(loadings)) 165 loadings <- rotated$loadings 166 Phi <- NULL} else { 167 if ((rotate=="promax")|(rotate=="Promax")) {pro <- Promax(loadings) 168 loadings <- pro$loadings 169 Phi <- pro$Phi} else { 170 if (rotate == "cluster") {loadings <- varimax(loadings)$loadings 171 pro <- target.rot(loadings) 172 loadings <- pro$loadings 173 Phi <- pro$Phi} else { 174 175 if (rotate =="oblimin"| rotate=="quartimin" | rotate== "simplimax") { 176 if (!require(GPArotation)) {warning("I am sorry, to do these rotations requires the GPArotation package to be installed") 177 Phi <- NULL} else { ob <- do.call(rotate,list(loadings) ) 178 loadings <- ob$loadings 179 Phi <- ob$Phi} 180 } 181 }}} 182 183 }} 184 #just in case the rotation changes the order of the factors, sort them 185 #added October 30, 2008 186 187 if(nfactors >1) { 188 ev.rotated <- diag(t(loadings) %*% loadings) 189 ev.order <- order(ev.rotated,decreasing=TRUE) 190 loadings <- loadings[,ev.order]} 191 rownames(loadings) <- colnames(r) 192 if(!is.null(Phi)) {Phi <- Phi[ev.order,ev.order] } #January 20, 2009 but, then, we also need to change the order of the rotation matrix! 193 class(loadings) <- "loadings" 194 if(nfactors < 1) nfactors <- n 195 196 result <- factor.stats(r,loadings,Phi,n.obs) #do stats as a subroutine common to several functions 197 198 result$communality <- round(diag(model),digits) 199 result$uniquenesses <- round(diag(r-model),digits) 200 result$values <- round(eigens$values,digits) 201 result$loadings <- loadings 202 203 if(!is.null(Phi)) {result$Phi <- Phi} 204 if(fm == "pa") result$communality.iterations <- round(unlist(comm.list),digits) 205 206 if(scores) {result$scores <- factor.scores(x.matrix,loadings) } 207 result$factors <- nfactors 208 result$fn <- "factor.pa" 209 result$Call <- cl 210 class(result) <- c("psych", "fa") 211 return(result) } 212 213 #modified October 30, 2008 to sort the rotated loadings matrix by the eigen values. 214 215