R/omega.R

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

"omegah" <-
function(m,nfactors=3,fm="minres",key=NULL,flip=TRUE, digits=2,title="Omega",sl=TRUE,labels=NULL, plot=TRUE,n.obs=NA,rotate="oblimin",Phi = NULL,option="equal",covar=FALSE,two.ok=FALSE,...) {
      #m is a correlation matrix, or if not, the correlation matrix is found
      #nfactors is the number of factors to extract
      #key allows items to be reversed scored  if desired
      #if Phi is not null, this implies that we have been given a factor matrix  -- added May 30, 2010
     if(!requireNamespace('GPArotation') && (rotate !="cluster")) {stop("I am sorry, you need to have the  GPArotation package installed")}
      cl <- match.call() 
      nvar <- dim(m)[2]
      raw.data <- NULL
      if(is.null(Phi)) {   #the normal case is to do the factor analysis of the raw data or the correlation matrix
      if(!isCorrelation(m)) {   #should also check for covariance matrix
                            n.obs <- dim(m)[1]
                            m <- as.matrix(m)
                            raw.data <- m #added 9/1/14
                            if(covar) {m <- cov(m,use="pairwise") } else {m <- cor(m,use="pairwise")}
                            } else {
                            if(!covar) m <- cov2cor(as.matrix(m))    #make sure it is a correlation matrix not a covariance or data matrix (if we change this, we will need to change the calculation for omega later)
                           }
                            
      if(is.null(colnames(m))) {  rownames(m) <- colnames(m) <- paste("V",1:nvar,sep="") }
       m.names <- colnames(m)
      
      if (!is.null(key)) { m <- diag(key) %*% m %*% diag(key)
                           colnames(m) <- m.names   #flip items if we choose to do so
                           flip <- FALSE   #we do this if we specify the key
                           } else {key <- rep(1,nvar) }
             signkey <- strtrim(key,1)
             signkey[signkey=="1"] <- ""
             m.names <- paste(m.names,signkey,sep="")
             colnames(m) <- rownames(m) <- m.names
      if ((nvar < 6) && (fm =="mle") ) {message(paste("In omega, 3 factors are too many for ",nvar," variables using mle.  Using minres instead",sep=""))
       fm <- "minres"} 
       } else { m.names <- rownames(m) }   #add the names if we have a factor input 
       
      gf <-schmid(m,nfactors,fm,digits,rotate=rotate,n.obs=n.obs,Phi=Phi,option=option,covar=covar,two.ok=two.ok, ...)
      if(!is.null(Phi)) { model <- m
                          nfactors <- dim(model)[2]
                          m <- factor.model(model,Phi=Phi,U2=FALSE)  #estimate the correlation matrix from the factor model
                          nvar <- dim(m)[2]
                          if(is.null(rownames(m))) {colnames(m) <- rownames(m) <- paste("V",1:nvar)}
                          }
                         
      gload <- gf$sl[,1]
      if (flip) {       #should we think about flipping items ?
       			key <- sign(gload)
      		 	key[key==0] <- 1     # a rare and weird case where the gloading is 0 and thus needs not be flipped
      			 if (sum(key) < nvar) {  #some items have negative g loadings and should be flipped  
            		 m <- diag(key) %*% m %*% diag(key)  #this is just flipping the correlation matrix so we can calculate alpha
            		 gf$sl[,1:(nfactors+1)] <-  diag(key) %*% gf$sl[,1:(nfactors+1)] 
            		 signkey <- strtrim(key,1)
            		 signkey[signkey=="1"] <- ""
            		 m.names <- paste(m.names,signkey,sep="")
            		 colnames(m) <- rownames(m) <- m.names
            		 rownames(gf$sl) <- m.names
          				
     		 
            }
          }
          Vt <- sum(m)   #find the total variance in the scale
     	  Vitem <- sum(diag(m)) 
     	  gload <- gf$sl[,1]	
       
     
      gsq <- (sum(gload))^2
      uniq <- sum(gf$sl[,"u2"])
      
      if((nfactors == 1) && (fm=="pc")) {gsq <- Vt - uniq
                                        warning("omega_h is not meaningful for a principal components analysis with one component")}   #weird condition when using fm=pc and 1 factor
      om.tot <- (Vt-uniq)/Vt
      om.limit <- gsq/(Vt-uniq)  
      alpha <- ((Vt-Vitem)/Vt)*(nvar/(nvar-1))
      sum.smc <- sum(smc(m,covar=covar))
      lambda.6 <- (Vt +sum.smc-sum(diag(m)))/Vt
      if (!is.null(digits)) {omega <-list(omega_h= gsq/Vt,alpha=alpha,lambda.6 = lambda.6,omega.tot =om.tot,schmid=gf ,key = key,title=title)
            dg <-max(digits-1,1)} else {
              omega <- list(omega_h= gsq/Vt,alpha=alpha,omega.tot=om.tot,schmid=gf,key=key,title=title)
      dg <- 1}
      ev <- colSums(gf$sl[,1:(nfactors+1)]^2)
      ECV <- ev[1]/sum(ev)
      omega.stats <- factor.stats(m,gf$sl[,1:(nfactors+1)],n.obs=n.obs)
      general.stats <- factor.stats(m,as.matrix(gf$sl[,1]),n.obs=n.obs)  #just get fit for the general factor
     if (nfactors<2) plot <- FALSE
    # if(require(Rgraphviz) && plot) {omega.model <-omega.graph(omega,title=title,sl=sl,labels=labels,digits=dg) } else {omega.model <- omega.sem(omega,sl=sl)}
   
     omega.model <- omega.sem(omega,sl=sl)
     
     #find the subset omegas
     omg <- omgo <- omt<-  rep(NA,nfactors+1)
     sub <- apply(gf$sl,1,function(x) which.max(abs(x[2:(nfactors+1)])))
     grs <- 0
     for(group in( 1:nfactors)) {
     groupi <- which(sub==group)
     if(length(groupi) > 0) {
      Vgr <- sum(m[groupi,groupi])
      gr <- sum(gf$sl[groupi,(group+1)])
      grs <- grs + gr^2
      omg[group+1] <- gr^2/Vgr
      omgo[group+1] <- sum(gf$sl[groupi,1])^2/Vgr
      omt[group+1] <- (gr^2+ sum(gf$sl[groupi,1])^2)/Vgr
     }
     omgo[1] <- sum(gf$sl[,1])^2/sum(m)  #omega h
     omg[1] <- grs/sum(m)  #omega of subscales
     omt[1] <- om.tot 
     om.group <- data.frame(total=omt,general=omgo,group=omg)
     rownames(om.group) <- colnames(gf$sl)[1:(nfactors+1)]

     }   #moved after tge bext line (6/21/18)
  #we should standardize the raw.data before doing the next step
  #change this to call factor.scores
  if(!is.null(raw.data)) {   #added mean imputation 4/22/23
       
        miss <- which(is.na(raw.data),arr.ind=TRUE)
    
       		item.means <- colMeans(raw.data,na.rm=TRUE)   #replace missing values with means
       		raw.data[miss]<- item.means[miss[,2]]
       		z.data <- scale(raw.data)
       		scores <- z.data %*% omega.stats$weights
       		} else {scores<- NULL} 
    # if(!is.null(raw.data)) {scores <- raw.data %*%  omega.stats$weights} else {scores<- NULL} 
    # }
     omega <- list(omega_h= gsq/Vt,omega.lim = om.limit,alpha=alpha,omega.tot=om.tot,G6=lambda.6,schmid=gf,key=key,stats = omega.stats,ECV=ECV,gstats = general.stats,call=cl,title=title,R = m,model=omega.model,omega.group=om.group,scores=scores,Call=cl)

      class(omega) <- c("psych","omega")
     if(plot)  omega.diagram(omega,main=title,sl=sl,labels=labels,digits=dg)
      return(omega)
      }
#April 4, 2011  added a check for fm=pc and nfactors == 1 to solve problem of omega_h < omega_t -- probably not a good idea.  removed
#January 9, 2014  added omega scores if the raw data are given


"omega" <- 
function(m,nfactors=3,fm="minres",n.iter=1,p=.05,poly=FALSE,key=NULL,flip=TRUE, digits=2,title="Omega",sl=TRUE,labels=NULL, plot=TRUE,n.obs=NA,rotate="oblimin",Phi = NULL,option="equal",covar=FALSE,...) {
 cl <- match.call()
  if(is.data.frame(m) || is.matrix(m)) {if((isCorrelation(m)) | (isCovariance(m) && covar)) {if(is.na(n.obs) && (n.iter > 1)) stop("You must specify the number of subjects if giving a correlation matrix")
                                # if(!require(MASS)) stop("You must have MASS installed to simulate data from a correlation matrix")
                                 }
                  }
 if(!is.data.frame(m) && !is.matrix(m)) {  n.obs=m$n.obs
                if(poly) {
                      pol <- list(rho=m$rho,tau = m$tau,n.obs=m$n.obs)
                        
                           m <- m$rho
                          } else { m <- m$R 
                          }
                          } else { #new data

          if(poly) { pol <- polychoric(m)
                     m <- pol$rho
                     n.obs <- pol$n.obs}  }       
            om <- omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,covar=covar,...) #call omega with the appropriate parameters
    

       
 if(is.na(n.obs) ) {n.obs <- om$stats$n.obs} 
 replicates <- list()
 if(n.iter > 1) {for (trials in 1:n.iter) {
 if(dim(m)[1] == dim(m)[2]) {#create data sampled from multivariate normal with correlation
                                      nvar <- dim(m)[1]
                                      #mu <- rep(0, nvar)
                                     # m <- mvrnorm(n = n.obs, mu, Sigma = m, tol = 1e-06, empirical = FALSE)
 #the next 3 lines replaces mvrnorm (taken from mvrnorm, but without the checks)
                                      eX <- eigen(m)
                                      m <- matrix(rnorm(nvar * n.obs),n.obs)
                                      m <-  t(eX$vectors %*% diag(sqrt(pmax(eX$values, 0)), nvar) %*%  t(m))
                            } else {m <- m[sample(n.obs,n.obs,replace=TRUE),]}
    if(poly) {pol <- polychoric(m)
            oms <- omegah(m=pol$rho,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=pol$n.obs,rotate=rotate,Phi = Phi,option=option,...) #call omega with the appropriate parameters
            } else {
 oms <- omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,...) #call omega with the appropriate parameters
 }                         
 # oms <-omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,...) #call fa with the appropriate parameters
  replicates[[trials]] <- list(omega=oms$omega_h,alpha=oms$alpha,omega.tot=oms$omega.tot,G6=oms$G6,omega.lim=oms$omega.lim)
  }
  
replicates <- matrix(unlist(replicates),ncol=5,byrow=TRUE)
z.replicates <- cbind(fisherz(replicates[,1:4]),replicates[,5])  #convert to z scores
means <- colMeans(z.replicates,na.rm=TRUE)
sds <-  apply(z.replicates,2,sd,na.rm=TRUE)
ci.lower <-  means + qnorm(p/2) * sds
ci.upper <- means + qnorm(1-p/2) * sds
ci <- data.frame(lower = ci.lower,upper=ci.upper)
ci <- rbind(fisherz2r(ci[1:4,]),ci[5,])
rownames(ci) <- c("omega_h","alpha","omega_tot","G6","omega_lim")
colnames(replicates) <- names(means) <- names(sds) <- rownames(ci) 
conf <- list(means = means,sds = sds,ci = ci,Call= cl,replicates=replicates)
om$Call=cl
results <- list(om = om,ci=conf) } else {om$Call=cl
                                         if(poly) {om$rho <- pol$rho
                                                   om$tau <- pol$tau
                                                   om$n.obs <- pol$n.obs
                                                   }
                                         results <- om}
class(results) <- c("psych","omega")
return(results)
 }
 #written April 25, 2011
 #adapted May 12, 2011 to be the primary version of omega