1"omega.bifactor" <- function(r,nfactors=3,rotate="bifactor",n.obs = NA,flip=TRUE,key=NULL,title="Omega from bifactor",...) {
2 cl <- match.call()
3 if(dim(r)[2] != dim(r)[1]) {n.obs <- dim(r)[1]
4                             r <- cor(r,use="pairwise")}
5 nvar <- dim(r)[2]
6  if(is.null(colnames(r))) {  rownames(r) <- colnames(r) <- paste("V",1:nvar,sep="") }
7       r.names <- colnames(r)
8
9 if (!is.null(key)) { r <- diag(key) %*% r %*% diag(key)
10                           colnames(r) <- r.names   #flip items if we choose to do so
11                           flip <- FALSE   #we do this if we specify the key
12                           } else {key <- rep(1,nvar) }
13
14 f <- fa(r,nfactors=nfactors,rotate=rotate,n.obs = n.obs)
15 if (flip) {       #should we think about flipping items ?
16       			key <- sign(f$loadings[,1])
17      		 	key[key==0] <- 1     # a rare and weird case where the gloading is 0 and thus needs not be flipped
18      			 if (sum(key) < nvar) {  #some items have negative g loadings and should be flipped
19            		 r <- diag(key) %*% r %*% diag(key)  #this is just flipping the correlation matrix so we can calculate alpha
20            		f$loadings <-  diag(key) %*% f$loadings
21            		 signkey <- strtrim(key,1)
22            		 signkey[signkey=="1"] <- ""
23            		 r.names <- paste(r.names,signkey,sep="")
24            		 colnames(r) <- rownames(r) <- r.names
25            		 rownames(f$loadings) <- r.names
26
27
28            }
29          }
30          Vt <- sum(r)   #find the total variance in the scale
31     	  Vitem <- sum(diag(r))
32     	  gload <- f$loadings[,1]
33
34
35       gsq <- (sum(gload))^2
36      uniq <- nvar - tr(f$loadings %*% t(f$loadings))
37
38      om.tot <- (Vt-uniq)/Vt
39      om.limit <- gsq/(Vt-uniq)
40      alpha <- ((Vt-Vitem)/Vt)*(nvar/(nvar-1))
41      sum.smc <- sum(smc(r))
42      lambda.6 <- (Vt +sum.smc-sum(diag(r)))/Vt
43      omega_h= gsq/Vt
44      omega <-list(omega_h= omega_h,alpha=alpha,G6 = lambda.6,omega.tot =om.tot ,key = key,title=title,f=f)
45      class(omega) <- c("psych","bifactor")
46      return(omega)
47      }
48
49
50
51
52