"omega.bifactor" <- function(r,nfactors=3,rotate="bifactor",n.obs = NA,flip=TRUE,key=NULL,title="Omega from bifactor",...) {
cl <- match.call()
if(dim(r)[2] != dim(r)[1]) {n.obs <- dim(r)[1]
r <- cor(r,use="pairwise")}
nvar <- dim(r)[2]
if(is.null(colnames(r))) { rownames(r) <- colnames(r) <- paste("V",1:nvar,sep="") }
r.names <- colnames(r)
if (!is.null(key)) { r <- diag(key) %*% r %*% diag(key)
colnames(r) <- r.names #flip items if we choose to do so
flip <- FALSE #we do this if we specify the key
} else {key <- rep(1,nvar) }
f <- fa(r,nfactors=nfactors,rotate=rotate,n.obs = n.obs)
if (flip) { #should we think about flipping items ?
key <- sign(f$loadings[,1])
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
r <- diag(key) %*% r %*% diag(key) #this is just flipping the correlation matrix so we can calculate alpha
f$loadings <- diag(key) %*% f$loadings
signkey <- strtrim(key,1)
signkey[signkey=="1"] <- ""
r.names <- paste(r.names,signkey,sep="")
colnames(r) <- rownames(r) <- r.names
rownames(f$loadings) <- r.names
}
}
Vt <- sum(r) #find the total variance in the scale
Vitem <- sum(diag(r))
gload <- f$loadings[,1]
gsq <- (sum(gload))^2
uniq <- nvar - tr(f$loadings %*% t(f$loadings))
om.tot <- (Vt-uniq)/Vt
om.limit <- gsq/(Vt-uniq)
alpha <- ((Vt-Vitem)/Vt)*(nvar/(nvar-1))
sum.smc <- sum(smc(r))
lambda.6 <- (Vt +sum.smc-sum(diag(r)))/Vt
omega_h= gsq/Vt
omega <-list(omega_h= omega_h,alpha=alpha,G6 = lambda.6,omega.tot =om.tot ,key = key,title=title,f=f)
class(omega) <- c("psych","bifactor")
return(omega)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.