Composite reliability (CFA)

3,080 views
Skip to first unread message

Jorge Sinval

unread,
May 22, 2017, 1:03:04 PM5/22/17
to lavaan
Hello!

I would like to know an easy way to calculate the composite reliability (CR) in order to easily produce my reports with Rmarkdown. I know that I can calculate with:

l<-inspect(fit,"coef")$lambda
v<-diag(inspect(fit,"coef")$theta)

#Composite Reliability
cr<-sum(l)^2/(sum(l)^2+sum(v))

But how can I calculate it for all factor with a single function? 

Thanks!

kma...@aol.com

unread,
May 23, 2017, 10:24:15 AM5/23/17
to lavaan
Jorge,

require(lavaan)

# The Holzinger and Swineford (1939) example
HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '


fit <- lavaan(HS.model, data=HolzingerSwineford1939,
              auto.var=TRUE, auto.fix.first=TRUE,
              auto.cov.lv.x=TRUE)

             

my.CR <- function(fit, ...){             
  l<-inspect(fit,"coef")$lambda
  v<-diag(inspect(fit,"coef")$theta)

  ll <- colSums(l)^2
  vm <- v %*% t(rep(1,length(ll)))
  sl <- sign(l^2)
  vv <- colSums(sl * vm)
  rr <- ll / (ll + vv)
  return(rr)
}

my.CR(fit)

   visual   textual     speed
0.6734319 0.8872648 0.8516294

If I am not mistaken, however, the equation that you gave only applies when the model has simple structure (no cross-loadings).  You would need to modify the function to handle more complex cases.

Keith
------------------------
Keith A. Markus
John Jay College of Criminal Justice, CUNY
http://jjcweb.jjay.cuny.edu/kmarkus
Frontiers of Test Validity Theory: Measurement, Causation and Meaning.
http://www.routledge.com/books/details/9781841692203/

Stephen Martin

unread,
May 23, 2017, 2:31:01 PM5/23/17
to lavaan
library(semTools)
reliability
(fit)

Jorge Sinval

unread,
May 23, 2017, 3:02:37 PM5/23/17
to lav...@googlegroups.com
Hello Keith!

I used this approach:

l<-inspect(fit.PBQ.2001.WLSMV,"coef")$lambda
l1<-l[1:12]

v<-diag(inspect(fit.PBQ.2001.WLSMV,"coef")$theta)
v1<-v[1:12]

#Composite Reliability
cr<-sum(l1)^2/(sum(l1)^2+sum(v1))
cr

But as you can see, I would have to change for each factor the interval of items. With your function, it works like a charm!

Thank you! Btw your book is excellent!


*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/TpWINJo_CRI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+unsubscribe@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at https://groups.google.com/group/lavaan.
For more options, visit https://groups.google.com/d/optout.

Jorge Sinval

unread,
May 23, 2017, 3:03:41 PM5/23/17
to lav...@googlegroups.com
Sorry Martin, the realiability function of semTools, only gives alfa, Omega1:3; and AVE, I want CR.

*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

--

Stephen Martin

unread,
May 23, 2017, 3:35:50 PM5/23/17
to lavaan
Yes, but, iirc, omega is a composite reliability. If you look at the formula for omega 1 (IIRC), it is the formula for composite reliability, but it includes the possibility of residual covariances or cross loadings. If you don't have residual covariances or cross loadings, the formula collapses to the CR formula.



On Tuesday, May 23, 2017 at 2:03:41 PM UTC-5, Jorge Sinval wrote:
Sorry Martin, the realiability function of semTools, only gives alfa, Omega1:3; and AVE, I want CR.

*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

On Tue, May 23, 2017 at 7:31 PM, Stephen Martin <hwki...@gmail.com> wrote:
library(semTools)
reliability
(fit)



On Monday, May 22, 2017 at 12:03:04 PM UTC-5, Jorge Sinval wrote:
Hello!

I would like to know an easy way to calculate the composite reliability (CR) in order to easily produce my reports with Rmarkdown. I know that I can calculate with:

l<-inspect(fit,"coef")$lambda
v<-diag(inspect(fit,"coef")$theta)

#Composite Reliability
cr<-sum(l)^2/(sum(l)^2+sum(v))

But how can I calculate it for all factor with a single function? 

Thanks!

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/TpWINJo_CRI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.

kma...@aol.com

unread,
May 23, 2017, 5:10:23 PM5/23/17
to lavaan
Jorge,
The equation that you gave assumes that the factor variances are fixed to 1 and thus drop out.  Stephen's post helped me to notice that this was not the case in the example from the lavaan() help page.  If you change the model specification accordingly, the CR values match Omega in the reliability() output.

+ my.CR(fit)
   visual   textual     speed
0.6253175 0.8851753 0.6877599
+ reliability(fit)
          visual   textual     speed     total
alpha  0.6261171 0.8827069 0.6884550 0.7604886
omega  0.6253175 0.8851753 0.6877599 0.8453349
omega2 0.6253175 0.8851753 0.6877599 0.8453349
omega3 0.6120040 0.8850604 0.6858414 0.8596195
avevar 0.3705584 0.7210162 0.4244882 0.5145871

Jorge Sinval

unread,
May 24, 2017, 11:01:57 AM5/24/17
to lavaan
Yes, you and Stephen are right, in the package says: 

"Note that this formula is modified from Fornell & Larcker (1981) in the case that factor variances are not 1."

However, to use the criteria for the Fornell and Larcker (1981) proposal, wich is more correct? The simplified formula? Or the Omega 1 values?

CR by Fornell and Larcker (1981):


kma...@aol.com

unread,
May 25, 2017, 12:57:29 AM5/25/17
to lavaan
Jorge,
I am not sure that I understand the question.  However, I notice that the equation you gave for CR is identical to equation 6.20b for omega given by McDonald 1999.  If both reliability estimates are conceptualized as factor-expained variance divided by total score variance then it is just a matter of two equations for the same quantity: one equation being a special case of a more general equation, based on simplifying assumptions.  When the assumptions are met, they give the same result and when the assumptions are violated the simplified equation is not appropriate.

Keith

McDonald, R. P. (1999).  Test theory: A unified treatment.  Mahwah, NJ:  Erlbaum.

Jorge Sinval

unread,
May 25, 2017, 7:57:04 AM5/25/17
to lav...@googlegroups.com
So the Fornell and Larcker (1981) CR is exactly the same thing that McDonald's ω... 
Inline image 1


By the way, it's interesting to notice that that McDonald (2011) didn't cite once Fornell and Larcker (1981).
When we are assessing the convergent validity, which is more correct reporting this values as CR or as McDonald's ω?

Which is the utility of the other two 
ωs that the semTools' reliability function give us?

Sorry for my ignorance

McDonald, R. P. (2011). Test theory: A unified treatment. New York, NY, USA: Routledge. https://doi.org/10.4324/9781410601087


*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/TpWINJo_CRI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+unsubscribe@googlegroups.com.

kma...@aol.com

unread,
May 25, 2017, 11:13:04 AM5/25/17
to lavaan
Jorge,
Your first question might be a good question for SEMNET.

The help file for the semTools reliability function provides a detailed description of the three omega coefficients that it returns.  It would be hard to top that.

https://www.rdocumentation.org/packages/semTools/versions/0.4-11/topics/reliability

The help file also provides additional references.

Keith
----------------------------------------------------------

Jorge Sinval

unread,
May 25, 2017, 11:24:21 AM5/25/17
to lav...@googlegroups.com
So, after looking at the documentation (and being rigorous), the omega that is equal to the CR formula is the Raykov (2001) omega (McDonald as you pointed, also suggests it). It attributes the omega 3 (hierarchical omega) to McDonald, and not the first one. Well, I think I will report the CR=ω as Raykov (2001). Would you consider it correct?

Thanks!

*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

Stephen Martin

unread,
May 25, 2017, 1:33:32 PM5/25/17
to lavaan
I'm honestly a bit confused about McDonald's omega. I understand hierarchical omega, essentially the reliability of a general factor after accounting for various subfactors. I'm not certain how the formula in reliability() achieves that, but I can say that hierarchical omega is not conceptually the same as omega1 and omega2.
Omega1 and omega2 are more similar to each other than they are to hierarchical omega.

All reliability is "true variance" over "some metric of total variance".
Omega1 and omega2 differ in how they specify total variance.
Omega1: Total variance is the explained variance from the factor of interest plus residual error plus residual covariance.
Omega2: Total variance is obtained from the full implied covariance matrix.

These can have important differences. As mentioned in the docs, Omega1's reliability uses variance after controlling for other factors. Omega2's reliability uses total implied variance from all factors. Omega1, then, is conceptually similar to "reliability of a measure after accounting for another factor", and omega2 is "reliability of a measure without controlling for another factor". I think, anyway.

Of these, I think omega2 makes much more sense, and I tend to prefer it. Composite reliability should be a metric of "if I were to compute sumscores/meanscores, what percent of variance is due to the underlying factor", and to me omega2 is more useful for this usecase, because if you *were* to compute sum scores, I don't think you would compute them after controlling for other factors, which omega1 is implicitly assuming.

Hierarchical omega, as I understand it, is quite different from these. By using the observed covariance matrix, it's likely that hierarchical omega will be a conservative estimate anyway. And hierarchical omega, I think anyway, is meant to be used for cases where you have subscales, but you want to know the reliability of measuring a general factor (bifactor model) if you were to sum across all items and treat the entire scale as one scale, not multiple subscales. If hierarchical omega is high, it means your scale appears to reliably measure one construct even after accounting for subscales. I'm not 100% sure of this, to be frank, but that question is not particularly interesting to me, so I tend to stick with omega 2.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.

kma...@aol.com

unread,
May 26, 2017, 11:23:29 AM5/26/17
to lavaan
The topic of this thread has drifted substantially but here are a few additional thoughts.

If the question is which estimates to report, there is a growing literature of simulation studies offering guidance.

If the question is how to report a given estimate, that was the question for which I suggested SEMNET.

There is probably some room for different perspectives.  I find omega3 attractive precisely because it will be more conservative.  If substantial difference emerge between the three omegas, that may indicate violations of the assumptions of omega1 and/or omega2.

The advantages gained from choosing one estimate over another may be relatively small in relation to the benefits of reporting an interval estimate instead of just a point estimate.  I am aware of three basic options:
1. Compute the reliability estimate as a defined parameter in the model syntax (you can check the point estimate against reliability() as a test of accuracy).
2. Use the boot package to bootstrap a standard error around the point estimate from reliability()
3. Use the ci.reliability() function from the MBESS package.  (A CI for alpha is also available from the rela package.)

Keith

Jorge Sinval

unread,
Jun 29, 2017, 5:34:54 AM6/29/17
to lav...@googlegroups.com
Thank you all, this was a very informative and clarifying discussion. As Keith pointed it drifted substantially, and just to finish, I have one additional contribution. McDonald (1999) presents "his" Omega has a measure to have evidence of discriminant and convergent validity (pp.211-213).

*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
 Nanos gigantum humeris insidentes.

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/TpWINJo_CRI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+unsubscribe@googlegroups.com.

João Marôco

unread,
May 29, 2018, 11:00:28 AM5/29/18
to lavaan
Dear Keith et al.

Using

fit <- lavaan(HS.model, data=HolzingerSwineford1939,
              auto.var=TRUE, auto.fix.first=TRUE,
              auto.cov.lv.x=F)

I believe that if you use non standardized coeficientes ([that is what you get with inspect(fit, "coef")$lambda],
your function should be

my.CR <- function(fit, ...){              
  l<-inspect(fit,"coef")$lambda
  v<-diag(inspect(fit,"coef")$theta)
  
  ll <- colSums(l)^2*diag(inspect(fit,"coef")$psi)
  vm <- v %*% t(rep(1,length(ll)))
  sl <- sign(l^2)
  vv <- colSums(sl * vm)
  rr <- ll / (ll + vv)
  return(rr)
}

As you can see, I am changing
ll <- colSums(l)^2
to
ll <- colSums(l)^2*diag(inspect(fit,"coef")$psi)


Can you please confirm this?
Best,
João

kma...@aol.com

unread,
May 30, 2018, 10:01:31 AM5/30/18
to lavaan
Joao & Jorge,
  Thanks for reviving this old thread and catching the inconsistency between the model specification and the omega calculation.  Jorge and I exchanged email about the thread off list yesterday.  For clarity, I decided to post three analyses of data provided as a worked example by Rod McDonald in his 1999 Test Theory book.

# McDonald1999, Table 6.3, Page 82
McDonald1999 <- c(2.5666, 1.560, 1.487, 1.195, 1.425, 1.560, 2.493, 1.283, 0.845, 1.313, 1.487, 1.283, 2.462, 1.127, 1.313, 1.195, 0.845, 1.127, 2.769, 1.323, 1.425, 1.313, 1.313, 1.323, 3.356)
McDonald1999Cov <- matrix(data=McDonald1999, ncol=5, nrow=5)
McDonald1999Names <- paste0('item', 1:5)
dimnames
(McDonald1999Cov) <- list(McDonald1999Names, McDonald1999Names)
McDonald1999Cov

# Unit loading
McDonaldModelLoading <- '
  LV =~ 1*item1 + item2 + item3 + item4 + item5
  LV ~~ LV
  item1 ~~ item1
  item2 ~~ item2
  item3 ~~ item3
  item4 ~~ item4
  item5 ~~ item5
'
# end model

# Unit variance
McDonaldModelVariance <- '
  LV =~ item1 + item2 + item3 + item4 + item5
  LV ~~ 1*LV
  item1 ~~ item1
  item2 ~~ item2
  item3 ~~ item3
  item4 ~~ item4
  item5 ~~ item5
'
# end model

# Fit models
fitUnitLoading
<- lavaan(model=McDonaldModelLoading, sample.cov=McDonald1999Cov, sample.nobs=215, estimator='ULS')
fitUnitVariance
<- lavaan(model=McDonaldModelVariance, sample.cov=McDonald1999Cov, sample.nobs=215, estimator='ULS')

# Extract estimates
c1
<-lavInspect(fitUnitLoading,"coef")
l1
<-c1$lambda
v1
<-diag(c1$theta)
p1
<-c1$psi[1,1]

c2
<-lavInspect(fitUnitVariance,"coef")
l2
<-c2$lambda
v2
<-diag(c2$theta)
p2
<-c2$psi[1,1]

c3
<-lavInspect(fitUnitLoading,"std")
l3
<-c3$lambda
v3
<-diag(c3$theta)
p3
<-c3$psi[1,1]

# Visually inspect estimates
l1
;l2;l3
v1
;v2;v3
p1
;p2;p3


# Compute omega
omega1
<- ((sum(l1)^2)*p1) / (((sum(l1)^2)*p1) + sum(v1))
omega2
<- (sum(l2)^2) / ((sum(l2)^2) + sum(v2))
omega3
<- (sum(l3)^2) / ((sum(l3)^2) + sum(v3)) # not appropriate

# McDonald 1999 page 90 calculates .8189 or .8188
omega1
; omega2; omega3

omega1b
<- ((sum(l2)^2)*p2) / (((sum(l2)^2)*p2) + sum(v2))
omega1b


The estimates are as follows.  Note that I used ULS estimation to match the worked example.

> omega1; omega2; omega3
[1] 0.8189322
[1] 0.8189322
[1] 0.8222746

The take away:
1. Method 1 (multiplying by the factor variance) is a general method that will work for either a unit loading or a unit factor variance.

2. Method 2 (no such multiplication) assumes that unit factor variance but will return the same result when that assumption is correct.

3. Method 3 uses standardized loading and is not appropriate as an estimate of raw omega because the raw total score effectively weights the items by their variances whereas the standardized solution gives them equal weights.  This method instead estimates the total score derived from converting the items to z scores before computing the total score.  This is analogous to standardized alpha.

João Marôco

unread,
May 30, 2018, 10:34:38 AM5/30/18
to lav...@googlegroups.com
Dear Keith, 
Thank you for getting back on this. I agree with you. Method 1 is the most appropriate as a general formula. That is why I made the suggestion for the correction. No offense meant. 
Warm regards, 

João Marôco
[Sent from my not that smart smartphone with more errors than usual]

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

Jorge Sinval

unread,
May 30, 2018, 12:23:12 PM5/30/18
to lav...@googlegroups.com
Thank you!
*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
"The someone's life is limited; the honor and respect last forever." - Miyamoto Musashi


You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/TpWINJo_CRI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.

kma...@aol.com

unread,
May 30, 2018, 8:45:38 PM5/30/18
to lavaan
Joao,
  No offense taken.  :-)

Chao Xu

unread,
May 31, 2018, 3:51:29 PM5/31/18
to lavaan
Hi Jorge,

Did you standardize the estimated parameters before you compute the composite reliability? The CR is just another way to say how much sample variances the factors account for.

Chao


On Monday, May 22, 2017 at 12:03:04 PM UTC-5, Jorge Sinval wrote:

Jorge Sinval

unread,
May 31, 2018, 5:51:12 PM5/31/18
to lavaan
Yes, used "std". 

Chao Xu

unread,
May 31, 2018, 6:13:23 PM5/31/18
to lavaan
Thanks for clarification.

Jorge Sinval

unread,
May 31, 2018, 6:19:12 PM5/31/18
to lav...@googlegroups.com
Use the formula that Professor Marôco sugested.
*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
"The someone's life is limited; the honor and respect last forever." - Miyamoto Musashi

--

Chao Xu

unread,
May 31, 2018, 6:27:03 PM5/31/18
to lavaan
I see.

It is also possible to specify what="std" in the inspect function. If your factor loading matrix is congeneric, then the squared value of each loading + its corresponding residual estimate should be 1. Of course this is the case of fitting model to polychoric correlation matrix if your data are ordinal.

Jorge Sinval

unread,
May 31, 2018, 6:31:58 PM5/31/18
to lav...@googlegroups.com
Yes, the formula of Fornell, and Larcker (1981).
*************************************************
Best Regards

Jorge Sinval (BS Psych, MS Psych)
Double PhD Student in Psychology (FFCLRP-USP & FPCEUP)
Mendeley Advisor
*************************************************
"The someone's life is limited; the honor and respect last forever." - Miyamoto Musashi

kma...@aol.com

unread,
May 31, 2018, 11:43:17 PM5/31/18
to lavaan
Chao,
I am not sure what the source of confusion is.  Do not use the standardized loadings unless you are interested in the reliability of a composite formed by summing the z scores of the item scores.  Under most circumstances, you will be interested in the reliability of a composite formed by summing the raw item scores.  For that, you need the raw loadings.

The simplified formula with raw loadings is perfectly okay so long as you fix the factor variance to one.  Otherwise, use the more general formula.
Reply all
Reply to author
Forward
0 new messages