Calculating average variance extracted (AVE) in R for checking discriminant validity (Fornell-Larcker criterion) We are using lavaan in R to calculate CFAs (confirmatory factor analyses) and SEMs (structural equation models). We now want to test whether two latent constructs can be assumed to be unrelated (discriminant validity).
Refering to 

Fornell, C. & Larcker, D. F. (1981). Evaluating structural equation models with unobservable variables and measurement error. Journal of Marketing Research, 18, 39–50.

we need to calculate the so called average variance extracted (AVE).
I found the following R-help-mailing list post. However, I don't think the output is the AVE I need for the Fornell-Larcker criterion. Values are calculated for each item and not factors.
Does anybody know how to obtain this value within lavaan or any other package?
If it is not possible to automatically obtain the AVE for each construct, can somebody point out how to calculate it based on the lavaan output?
Help would be much appreciated! Thanks.
 A: As it seems the semTools package now provides a function that calculates the AVE based on a model estimated in lavaan:
reliability()
However, the package description notes that the calculation slightly differs from the one proposed in Fornell & Larcker (1981):
"Note that this formula is modified from Fornell & Larcker (1981) 
 in the case that factor variances are not 1."

I am not quite sure what this means for interpretation of the AVE.
A: I developed a function to calculate CR and AVE for the sem package. I am sure you can adapt to lavaan. 
Code follows:
crave <- function(model, digits=2) {
# ig111207
# computes composite reliability and average variance extracted
# usage: crave(model)
# where model is a sem fitted model
# such as sem.wh.1
require(sem)
x <- stdCoef(model, twoheaded=F)
coeff <- x[,2]
paths <- as.character(x[,3])
newl <- strsplit(paths, " <--- ")
newl <- as.data.frame(newl)
newl <- t(newl)
item <- newl[,1]
LV <- newl[,2]

coeff2 <- coeff^2
e2 <- 1-coeff2

y <- data.frame(LV, item, coeff, coeff2, e2)
rownames(y) <- item

cr <- by(y, LV, function(x) sum(x$coeff)^2 / (sum(x$coeff)^2+sum(x$e2)) )
    ave <- by(y, LV, function(x) sum(x$coeff2) / (sum(x$coeff2)+sum(x$e2)) )

z <- data.frame(
  cr=round(as.numeric(cr), digits),
  s1=ifelse(cr<.7,"","*"),
  ave=round(as.numeric(ave), digits),
  s2=ifelse(ave<.5,"","*")
)  

#print(y,digits)
#print(z)
answer <- list(items=y,latent=z)
answer

}

I hope this helps.
Iuri.
