R/factor.scores.R

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

"factor.scores" <- function(x,f,Phi=NULL,method=c("Thurstone","tenBerge","Anderson","Bartlett","Harman","components"),rho=NULL,missing=FALSE,impute="none") {
#the normal case is f is the structure matrix and Phi is not specified
#Note that the Grice formulas distinguish between Pattern and Structure matrices
#I need to confirm that I am doing this

    if(length(method) > 1) method <- "tenBerge"   #the default
    if(method=="regression") method <- "Thurstone"
    if(method %in% c("tenberge", "Tenberge","tenBerge","TenBerge")) method <- "tenBerge"
    if(length(class(f)) > 1) { if(inherits(f[2] ,"irt.fa" )) f <- f$fa  }
    
     if(!is.matrix(f)) {Phi <- f$Phi
     f <- loadings(f)
      if(ncol(f)==1) {method <- "Thurstone"}
      }
     nf <- dim(f)[2]
      if(is.null(Phi)) Phi <- diag(1,nf,nf)
     if(dim(x)[1] == dim(f)[1]) {r <- as.matrix(x)
         square <- TRUE} else { 
          square <- FALSE
         if(!is.null(rho)) {r <- rho } else {
          r <- cor(x,use="pairwise") #find the correlation matrix from the data
      }}
      
      S <- f %*% Phi   #the Structure matrix 
   switch(method,   
    "Thurstone" = { w <- try(solve(r,S),silent=TRUE )  #these are the factor weights (see Grice eq. 5)
     	if(inherits(w,"try-error")) {message("In factor.scores, the correlation matrix is singular, the pseudo inverse is  used")
              # r <- cor.smooth(r)
              w  <- Pinv(r) %*% S }
        
     # w <- try(solve(r,S),silent=TRUE)    #redudant,  this will fail again   (pointed out by Chanlder McClellan 4/4/23)
     # if(inherits(w,"try-error")) {message("I was unable to calculate the factor score weights, factor loadings used instead")
     #          w <- f}
      colnames(w) <- colnames(f)
      rownames(w) <- rownames(f)
       }, 
      
  "tenBerge" = { #Following Grice equation 8 to estimate scores for oblique solutions (with a correction to the second line where r should r.inv
        L <- f %*% matSqrt(Phi)
        r.5 <- invMatSqrt(r)
       
        r <- cor.smooth(r)
        inv.r <- try(solve(r),silent=TRUE)
        if(inherits(inv.r, as.character("try-error")))  {warning("The tenBerge based scoring could not invert the correlation matrix, regression scores found instead")
                                                      ev <- eigen(r)
      ev$values[ev$values < .Machine$double.eps] <- 100 * .Machine$double.eps
        r <- ev$vectors %*% diag(ev$values) %*% t(ev$vectors)
        diag(r)  <- 1
       w <- solve(r,f)}  else {
        C <- r.5 %*% L %*% invMatSqrt(t(L) %*% inv.r %*% L)    #note that this is the correct formula, per Grice personal communication
        w <- r.5 %*% C %*% matSqrt(Phi)}
        colnames(w) <- colnames(f)
        rownames(w) <- rownames(f)
        },

 
       
       "Harman" = { #Grice equation 10 -- 
     #   m <- t(f)  %*% f  #factor intercorrelations 
     m <- f %*% t(S)  #should be this  (the model matrix)  Revised August 31, 2017
     diag(m) <- 1  #Grice does not say this, but it is necessary to make it work!
       inv.m <- solve(m)
     #  w <- f %*%inv.m  
     w <- inv.m %*% f
       }, 
       
       
        
    "Anderson" =  { #scores for orthogonal factor solution will be orthogonal  Grice Eq 7 and 8
    I <- diag(1,nf,nf)
    h2 <-  diag( f %*% Phi %*% t(f))
    U2 <- 1 - h2
    inv.U2 <- diag(1/U2)
    w <- inv.U2 %*% f %*% invMatSqrt(t(f) %*% inv.U2 %*% r %*% inv.U2 %*% f)
    colnames(w) <- colnames(f)
    rownames(w) <- rownames(f)
    },
    
   "Bartlett" = {    #Grice eq 9  # f should be the pattern, not the structure 
    I <- diag(1,nf,nf)
    h2 <-  diag( f %*% Phi %*% t(f))
    U2 <- 1 - h2
    inv.U2 <- diag(1/U2)
    w <- inv.U2 %*% f %*% (solve(t(f) %*% inv.U2 %*% f))
    colnames(w) <- colnames(f)
    rownames(w) <- rownames(f)
    },
    "none" = {w <- NULL},
    
    "components" = {w <- try(solve(r,f),silent=TRUE )    #basically, just do the regression/Thurstone approach for components
                    w <- f }
    )
    
    
    #now find a few fit statistics
    if(is.null(w)) {results <- list(scores=NULL,weights=NULL)} else {
     R2 <- diag(t(w) %*% S)  #this had been   R2 <- diag(t(w) %*% f)   Corrected Sept 1, 2017
     if(any(R2 > 1) || (prod(!is.nan(R2)) <1) || (prod(R2) < 0) ) {#message("The matrix is probably singular -- Factor score estimate results are likely incorrect")
                      R2[abs(R2) > 1] <- NA
                      R2[R2 <= 0] <- NA
                     }
     #if ((max(R2,na.rm=TRUE) > (1 + .Machine$double.eps)) ) {message("The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.")}
      r.scores <- cov2cor(t(w) %*% r %*% w) #what actually is this?
     
    
  if(square) {  #that is, if given the correlation matrix
     class(w) <- NULL
     results <- list(scores=NULL,weights=w)
      results$r.scores <- r.scores 
   	  results$R2 <- R2   #this is the multiple R2 of the scores with the factors
     } else {
     

         #missing <- rowSums(is.na(x))
    if(missing && (impute !="none")) {
       x <- data.matrix(x)
        miss <- which(is.na(x),arr.ind=TRUE)
        if(impute=="mean") {
       		item.means <- colMeans(x,na.rm=TRUE)   #replace missing values with means
       		x[miss]<- item.means[miss[,2]]} else { 
       		item.med   <- apply(x,2,median,na.rm=TRUE) #replace missing with medians
        	x[miss]<- item.med[miss[,2]]}   #this only works if items is a matrix
     }
      
     if(method !="components") {if(impute!="none")  {scores <- factorScoresSapa(weights=w,  items = scale(x))} else {
    
     scores <- scale(x) %*% w }} else {  #standardize the data before doing the regression if using factors, 
        scores <- x %*% w}       # for components, the data have already been zero centered and, if appropriate, scaled
     results <- list(scores=scores,weights=w)
     results$r.scores <- r.scores
     results$missing <- missing 
   	  results$R2 <- R2   #this is the multiple R2 of the scores with the factors
     }
     }
   
     return(results) }
     #how to treat missing data?  see score.item
    
     
     
     
"matSqrt" <- function(x) {
   e <- eigen(x)
    e$values[e$values < 0] <- .Machine$double.eps
   sqrt.ev <- sqrt(e$values)   #need to put in a check here for postive semi definite
   result <- e$vectors %*% diag(sqrt.ev) %*% t(e$vectors)
   result}
   
   
"invMatSqrt" <- function(x) {
   e <- eigen(x)
   if(is.complex(e$values)) {warning("complex eigen values detected by invMatSqrt, results are suspect")
                 result <- x
      } else {
      
       e$values[e$values < .Machine$double.eps] <- 100 * .Machine$double.eps
   inv.sqrt.ev <- 1/sqrt(e$values)   #need to put in a check here for postive semi definite
   result <- e$vectors %*% diag(inv.sqrt.ev) %*% t(e$vectors) }
   result}