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