# A function to create a correlation matrix with a hierarchical structure
# The default values match those of Jensen and Weng
#dropped the call to mvrnorm Nov 28, 2014
#and added back the mu parameter in Sept,2016 (thanks to Alan Robinson for noticing this)
"sim.hierarchical" <-
function (gload=NULL,fload=NULL,n=0,raw=FALSE,mu = NULL) {
cl <- match.call()
# require(MASS)
if(is.null(gload)) gload=matrix(c(.9,.8,.7),nrow=3)
if(is.null(fload)) {fload <-matrix(c(.8,.7,.6,rep(0,9),.7,.6,.5,rep(0,9),.6,.5,.4), ncol=3)}
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
nvar <- dim(fload)[1]
colnames(model) <- rownames(model) <- paste("V",1:nvar,sep="")
if(n>0) {
if(is.null(mu)) mu <- rep(0,nvar)
#observed <- mvrnorm(n = n, mu, Sigma=model, tol = 1e-6, empirical = FALSE)
#the next 3 lines replaces mvrnorm (adapted from mvrnorm, but without the checks)
eX <- eigen(model)
observed <-matrix(rnorm(nvar * n),n)
observed <- t( eX$vectors %*% diag(sqrt(pmax(eX$values, 0)), nvar) %*% t(observed))
observed <- t(t(observed) + mu)
colnames(observed) <- paste("V",1:nvar,sep="")
r <- cor(observed)
if(!raw) {
result <- list(model=model, r= r, N=n,Call=cl) } else {
result <- list(model=model, r= r,observed=observed, N=n,Call=cl)}
class(result) <- c("psych","sim")
return(result)
} else {return( model) }
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.