McDonald has proposed coefficient omega as an estimate of the general factor saturation of a test. Zinbarg, Revelle, Yovel and Li (2005) compare McDonald's Omega to Chronbach's alpha and Revelle's beta. They conclude that omega is the best estimate.
To find omega it is necessary to do a factor analysis of the original data set, rotate the factors obliquely, do a Schmid Leiman transformation, and then find omega. Here we present code to do that. This code is included in the psych package of routines for personality research that may be loaded from the repository http://personality-project.org/r. Note that this package is my first attempt at building packages and leaves a lot to be desired.
#R Functions to find McDonald's coefficient omega as well as Cronbach's alpha
#Omega is an estimate of the general factor saturation of a test
#Alpha is an estimate of the mean split half reliability of a test and in the case
#of no group factors, is also an estimate of the general factor saturation of a test
#see Zinbarg, Revelle, Yovel, and Li, Psychometrika, 2005
#
#In addition, we provide a function to do a Schmid Leiman transformation
#of a data matrix to find general and orthogonal group factor loadings.
#(See Schmid, J. & Leiman, J.M. (1957), Psychometrika.)
#
#An additional routine generates hierarchical structured matrices to follow
#an example by Jensen and Weng, Intelligence, 1994
#William Revelle
#revelle@northwestern.edu
#May 29, 2005
#updated November 5, 2005
#
#Requires the GPArotation package (downloadable from CRAN)
#
#
######
#Omega
######
"omega" <-
function(m,nfactors=3,pc="pa",...) {
#m is a correlation matrix, or if not, the correlation matrix is found
#nfactors is the number of factors to extract
require(GPArotation)
nvar <-dim(m)[2]
if(dim(m)[1] != dim(m)[2]) m <- cor(m,use="pairwise")
gf<-schmid(m,nfactors,pc,...)
Vt <- sum(m) #find the total variance in the scale
Vitem <-sum(diag(m)) #
gload <- gf$sl[,1]
gsq <- (sum(gload))^2
alpha <- ((Vt-Vitem)/Vt)*(nvar/(nvar-1))
omega <- list(omega= gsq/Vt,alpha=alpha,schmid=gf)
}
#assuming all variables in the model are positively loaded on g, then omega is the ratio of the squared sum
#of the gloadings / the sum of the correlations
#######
#Schmid Leiman
#######
#corrected estimate of communality, May 21, 2007
"schmid" <-
function (model, nfactors = 3, pc = "pa",...)
{
#model is a correlation matrix, or if not, the correlation matrix is found
#nfactors is the number of factors to extract
require(GPArotation)
nvar <-dim(model)[2]
if(dim(model)[1] != dim(model)[2]) model <- cor(model,use="pairwise")
if (pc=="pc") {
fact <- principal(model, nfactors,...)
} else {if (pc=="pa") {fact <- factor.pa(model, nfactors,...) } else {
fact <- factanal(x, covmat = model, factors = nfactors,...)
}}
orth.load <- loadings(fact)
obminfact <- oblimin(orth.load)
rownames(obminfact$loadings) <- attr(model,"dimnames")[[1]]
fload <- obminfact$loadings
factr <- t(obminfact$Th) %*% (obminfact$Th)
gfactor <- factanal(x, covmat = factr, factors = 1)
gload <- loadings(gfactor)
gprimaryload <- fload %*% gload
colnames(gprimaryload) <- "g factor"
u2 <- 1 - diag(orth.load %*% t(orth.load))
h2 <- 1 - u2
uniq <- 1 - fload^2
guniq <- 1 - gprimaryload^2
Ig <- matrix(0, ncol = nfactors, nrow = nfactors)
diag(Ig) <- gload
primeload <- fload %*% Ig
uniq2 <- 1 - uniq - primeload^2
sm <- sqrt(uniq2)
schmid <- list(sl = cbind(gprimaryload, sm, h2, u2), orthog = fact$loadings, oblique = fload,
fcor = factr, gloading = gload)
}
######
#generate a schmid leiman hierarchical matrix following Jensen and Weng, 1994
make.hierarchical <- function (x) {
gload<-matrix(c(.9,.8,.7),nrow=3) # a higher order factor matrix
fload <-matrix(c( #a lower order (oblique) factor matrix
.8,0,0,
.7,0,.0,
.6,0,.0,
0,.7,.0,
0,.6,.0,
0,.5,0,
0,0,.6,
0,0,.5,
0,0,.4), ncol=3,byrow=TRUE)
fcor <- gload %*% t(gload) #the factor correlation matrix
diag(fcor) <-1 #put ones on the diagonal
model <- fload%*% fcor %*% t(fload) #the model correlation matrix for oblique factors
diag(model)<- 1 # put ones along the diagonal
make.hierarchical <- model }
#sample output for the Jensen and Weng problem
round(sm,2)
g factor Factor1 Factor2 Factor3 h2 u2
[1,] 0.72 0.35 0.00 0.00 0.64 0.36
[2,] 0.63 0.31 0.00 0.00 0.49 0.51
[3,] 0.54 0.26 0.00 0.00 0.36 0.64
[4,] 0.56 0.00 0.42 0.00 0.49 0.51
[5,] 0.48 0.00 0.36 0.00 0.36 0.64
[6,] 0.40 0.00 0.30 0.00 0.25 0.75
[7,] 0.42 0.00 0.00 0.43 0.36 0.64
[8,] 0.35 0.00 0.00 0.36 0.25 0.75
[9,] 0.28 0.00 0.00 0.29 0.16 0.84
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1.00 0.56 0.48 0.40 0.35 0.29 0.30 0.25 0.20
[2,] 0.56 1.00 0.42 0.35 0.30 0.25 0.26 0.22 0.18
[3,] 0.48 0.42 1.00 0.30 0.26 0.22 0.23 0.19 0.15
[4,] 0.40 0.35 0.30 1.00 0.42 0.35 0.24 0.20 0.16
[5,] 0.35 0.30 0.26 0.42 1.00 0.30 0.20 0.17 0.13
[6,] 0.29 0.25 0.22 0.35 0.30 1.00 0.17 0.14 0.11
[7,] 0.30 0.26 0.23 0.24 0.20 0.17 1.00 0.30 0.24
[8,] 0.25 0.22 0.19 0.20 0.17 0.14 0.30 1.00 0.20
[9,] 0.20 0.18 0.15 0.16 0.13 0.11 0.24 0.20 1.00
> sum(model)
[1] 27.9762
> 19.18435/27.9762
[1] 0.6857382
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.64 0.56 0.48 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.56 0.49 0.42 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.48 0.42 0.36 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.49 0.42 0.35 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.42 0.36 0.30 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.35 0.30 0.25 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.36 0.30 0.24
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.30 0.25 0.20
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.24 0.20 0.16
gload<-matrix(c(.9,.8,.7),nrow=3) # a higher order factor matrix
fload <-matrix(c( #a lower order (oblique) factor matrix
.8,0,0,
.7,0,.0,
.6,0,.0,
0,.7,.0,
0,.6,.0,
0,.5,0,
0,0,.6,
0,0,.5,
0,0,.4), ncol=3,byrow=TRUE)
fcor <- gload %*% t(gload) #the factor correlation matrix
diag(fcor) <-1 #put ones on the diagonal
model <- fload%*% fcor %*% t(fload) #the model correlation matrix for oblique factors
diag(model)<- 1 # put ones along the diagonal