#adapted from my structure.diagram function to draw the output from esem #August 4, 2016 "esemDiagram" <- "esem.diagram" <- function(esem=NULL,labels=NULL,cut=.3,errors=FALSE,simple=TRUE,regression=FALSE,lr=TRUE, digits=1,e.size=.1,adj=2,main="Exploratory Structural Model", ...){ cancor <- FALSE #a helper function sort.f <- function(x) { nvar <- ncol(x) if(is.null(nvar)) {return(x)} else { nitem <- nrow(x) cluster <- data.frame(item <- seq(1:nitem),clust = rep(0,nitem)) cluster$clust <- apply(abs(x),1,which.max) ord <- sort(cluster$clust,index.return=TRUE) x[1:nitem,] <- x[ord$ix,] rownames(x) <- rownames(x)[ord$ix] return(x)}} #first some default values if(length(class(esem)) > 1) { names <- cs(lm, lavaan,fa,principal, esem,mediate, lmCor, setCor, fa.multi,lm) value <- inherits(esem,names,which=TRUE) # value <- class(x)[2] if(any(value > 0) ) { value <- names[which(value > 0)]} else {value <- "other"} if(value=="lmCor" ) {cancor<- TRUE} else {cancor<- FALSE} } xmodel <- sort.f(esem$loadsX) ymodel <- sort.f(esem$loadsY) Phi <- esem$Phi num.y <- num.x <- 0 #we assume there is nothing there vars<- NULL num.xvar <- dim(xmodel)[1] #how many x variables? num.yvar <- dim(ymodel)[1] num.xfactors <- dim(xmodel)[2] num.yfactors <- dim(ymodel)[2] if(is.null(num.xvar)) num.xvar <- length(xmodel) if(is.null(num.yvar)) num.yvar <- length(ymodel) if(is.null(num.xfactors)) num.xfactors <- 1 if(is.null(num.yfactors)) num.yfactors <- 1 # if(max(num.xvar,num.yvar) < 10) e.size <- e.size* 2 #make the ellipses bigger for small problems e.size <- e.size * 10/ max(num.xvar,num.yvar) if(is.null(labels)) { xvars <- rownames(xmodel)} else { xvars <- vars <- labels} if(is.null(ncol(xmodel))) xvars <- names(xmodel) # if(is.null(vars) ) {xvars <- paste0("x",1:num.xvar) } fact <- colnames(xmodel) if (is.null(fact)) { fact <- paste0("X",1:num.xfactors) } if(is.null(ncol(ymodel))) {yvars <-names(ymodel) } else {yvars <- rownames(ymodel)} if(is.null(yvars)) {yvars <- paste0("y",1:num.y) } if(is.null(labels)) {vars <- c(xvars,yvars)} else {yvars <- labels[(num.xvar+1):(num.xvar+num.y)]} yfact <- colnames(ymodel) if(is.null(yfact)) {yfact <- paste0("Y",1:num.yfactors) } fact <- c(fact,yfact) num.var <- num.xvar + num.y num.factors <- num.xfactors + num.yfactors sem <- matrix(rep(NA),6*(num.var*num.factors + num.factors),ncol=3) #this creates an output model for sem analysis lavaan <- vector("list",num.xfactors + num.yfactors) #create a list for lavaan colnames(sem) <- c("Path","Parameter","Value") var.rect <- list() fact.rect <- list() # length.labels <- 0 # a filler for now #plot.new() is necessary if we have not plotted before #strwd <- try(strwidth(xvars),silent=TRUE) strwd <- try(length.labels <- max(strwidth(xvars),strwidth(yvars),strwidth("abc"))/1.8,silent=TRUE) #although this throws an error if the window is not already open, we don't show it #if (class(strwd) == "try-error" ) {plot.new() } if (inherits(strwd, "try-error" )) {length.labels = max(nchar(xvars),3)/1.8 } #length.labels <- max(strwidth(xvars),strwidth("abc"))/1.8 if(lr) {limx <- c(-(length.labels),max(num.xvar,num.yvar)+3 ) limy <- c(0,max(num.xvar,num.yvar)+1) } else { limy <- c(-(length.labels),max(num.xvar,num.yvar) +3 ) limx <- c(0,max(num.xvar,num.yvar)+1) if( errors) limy <- c(-1,max(num.xvar,num.yvar)+2)} scale.xaxis <- 3 if(lr) {plot(0,type="n",xlim=limx,ylim=limy,frame.plot=FALSE,axes=FALSE,ylab="",xlab="",main=main)} else {plot(0,type="n",xlim=limx,ylim=limy,frame.plot=FALSE,axes=FALSE,ylab="",xlab="",main=main) } #now draw the x part k <- num.factors x.adjust <- num.xvar/ max(num.xvar,num.yvar) for (v in 1:num.xvar) { if(lr) { var.rect[[v]] <- dia.rect(0,(num.xvar-v+1)/x.adjust,xvars[v],xlim=limx,ylim=limy,...) } else { var.rect[[v]] <- dia.rect(v,0,xvars[v],xlim=limy,ylim=limx,...) } } nvar <- num.xvar f.scale <- (num.xvar+ 1)/(num.xfactors+1)/x.adjust factors <- round(xmodel,digits) if(is.null(ncol(factors))) factors <- matrix(factors) if (num.xfactors >0) { for (f in 1:num.xfactors) { if(!regression) {if(lr) {fact.rect[[f]] <- dia.ellipse(limx[2]/scale.xaxis,(num.xfactors+1-f)*f.scale,fact[f],xlim=c(0,nvar),ylim=c(0,nvar),e.size=e.size,...)} else {fact.rect[[f]] <- dia.ellipse(f*f.scale,limy[2]/scale.xaxis,fact[f],ylim=c(0,nvar),xlim=c(0,nvar),e.size=e.size,...) } } else {if(lr) {fact.rect[[f]] <- dia.rect(limx[2]/scale.xaxis,(num.xfactors+1-f)*f.scale,fact[f],xlim=c(0,nvar),ylim=c(0,nvar),...)} else { fact.rect[[f]] <- dia.rect(f*f.scale,limy[2]/scale.xaxis,fact[f],xlim=c(0,nvar),ylim=c(0,nvar),...)} } for (v in 1:num.xvar) { if(is.numeric(factors[v,f])) { if(simple && (abs(factors[v,f]) == max(abs(factors[v,])) ) && (abs(factors[v,f]) > cut) | (!simple && (abs(factors[v,f]) > cut))) { if (!regression) {if(lr){dia.arrow(from=fact.rect[[f]],to=var.rect[[v]]$right,labels =factors[v,f],col=((sign(factors[v,f])<0) +1),lty=((sign(factors[v,f])<0) +1),adj=f %% adj +1) } else {dia.arrow(from=fact.rect[[f]],to=var.rect[[v]]$top,labels =factors[v,f],col=((sign(factors[v,f])<0) +1),lty=((sign(factors[v,f])<0) +1)) } } else {dia.arrow(to=fact.rect[[f]]$left,from=var.rect[[v]]$right,labels =factors[v,f],col=((sign(factors[v,f])<0) +1))} } } else { if (factors[v,f] !="0") { if (!regression) { if(lr) {dia.arrow(from=fact.rect[[f]],to=var.rect[[v]]$right,labels =factors[v,f]) } else {dia.arrow(from=fact.rect[[f]],to=var.rect[[v]]$top,labels =factors[v,f])} } else {if(lr) {dia.arrow(to=fact.rect[[f]],from=var.rect[[v]]$right,labels =factors[v,f])} else {dia.arrow(to=fact.rect[[f]],from=var.rect[[v]]$top,labels =factors[v,f])} } } } } } if (num.xfactors ==1) { lavaan[[1]] <- paste(fact[1],"=~ ") for(i in 1:num.xvar) { sem[i,1] <- paste(fact[1],"->",vars[i],sep="") lavaan[[1]] <- paste0(lavaan[[1]], ' + ', vars[i]) if(is.numeric(factors[i])) {sem[i,2] <- vars[i]} else {sem[i,2] <- factors[i] } }} #end of if num.xfactors ==1 k <- num.xvar+1 k <- 1 for (f in 1:num.xfactors) { #if (!is.numeric(factors[i,f]) || (abs(factors[i,f]) > cut)) lavaan[[f]] <- paste0(fact[f] ," =~ ") for (i in 1:num.xvar) { if((!is.numeric(factors[i,f] ) && (factors[i,f] !="0"))|| ((is.numeric(factors[i,f]) && abs(factors[i,f]) > cut ))) { sem[k,1] <- paste(fact[f],"->",vars[i],sep="") lavaan[[f]] <- paste0(lavaan[[f]], ' + ', vars[i]) if(is.numeric(factors[i,f])) {sem[k,2] <- paste("F",f,vars[i],sep="")} else {sem[k,2] <- factors[i,f]} k <- k+1 } #end of if } } } #end of if num.xfactors >0 if(errors) { for (i in 1:num.xvar) {if(lr) { dia.self(var.rect[[i]],side=2) } else { dia.self(var.rect[[i]],side=1)} sem[k,1] <- paste(vars[i],"<->",vars[i],sep="") sem[k,2] <- paste("x",i,"e",sep="") k <- k+1 } } #do it for y model if(!is.null(ymodel)) { if(lr) { y.adj <- num.yvar/ max(num.xvar,num.yvar) #y.adj <- num.yvar/2 - num.xvar/2 f.yscale <- limy[2]/(num.yfactors+1) y.fadj <- 0} else { # y.adj <- num.xvar/2 - num.yvar/2 y.adj <- num.yvar/2 f.yscale <- limx[2]/(num.yfactors+1) y.fadj <- 0} for (v in 1:num.yvar) { if(lr){ var.rect[[v+num.xvar]] <- dia.rect(limx[2]-.35,limy[2]-v /y.adj,yvars[v],xlim=limx,ylim=limy,...)} else { var.rect[[v+num.xvar]] <- dia.rect(v / y.adj,limx[2],yvars[v],xlim=limy,ylim=limx,...)} } } #we have drawn the y variables, now should we draw the Y factors if(!is.null(ymodel)){ y.factors <- round(ymodel,digits) if(is.null(ncol(y.factors))) { y.factors <- matrix(y.factors) } num.y <- nrow(y.factors) for (f in 1:num.yfactors) {if(lr) { fact.rect[[f+num.xfactors]] <- dia.ellipse(2*limx[2]/scale.xaxis,(num.yfactors+1-f)*f.yscale +y.fadj,yfact[f],xlim=c(0,nvar),ylim=c(0,nvar),e.size=e.size,...)} else { fact.rect[[f+num.xfactors]] <- dia.ellipse(f*f.yscale+ y.fadj,2*limx[2]/scale.xaxis,yfact[f],ylim=c(0,nvar),xlim=c(0,nvar),e.size=e.size,...)} for (v in 1:num.yvar) {if(is.numeric(y.factors[v,f])) { {if(simple && (abs(y.factors[v,f]) == max(abs(y.factors[v,])) ) && (abs(y.factors[v,f]) > cut) | (!simple && (abs(factors[v,f]) > cut))) { if(lr) { dia.arrow(from=fact.rect[[f+num.xfactors]],to=var.rect[[v+num.xvar]]$left,labels =y.factors[v,f],col=((sign(y.factors[v,f])<0) +1),lty=((sign(factors[v,f])<0) +1),adj = f %% adj +1)} else { dia.arrow(from=fact.rect[[f+num.xfactors]],to=var.rect[[v+num.xvar]]$bottom,labels =y.factors[v,f],col=((sign(y.factors[v,f])<0) +1),lty=((sign(factors[v,f])<0) +1))} } } } else {if(factors[v,f] !="0") {if(lr) {dia.arrow(from=fact.rect[[f+num.xfactors]],to=var.rect[[v+num.xvar]]$left,labels =y.factors[v,f] ,lty=((sign(factors[v,f])<0) +1))} else { dia.arrow(from=fact.rect[[f+num.xfactors]],to=var.rect[[v+num.xvar]]$bottom,labels =y.factors[v,f],lty=((sign(factors[v,f])< 0) +1)) } } }} } if (num.yfactors ==1) { lavaan[[num.xfactors +1 ]] <- paste(fact[num.xfactors +1], "=~") for (i in 1:num.y) { sem[k,1] <- paste(fact[1+num.xfactors],"->",yvars[i],sep="") lavaan[[num.xfactors +1]] <- paste0(lavaan[[num.xfactors +1]], ' + ', yvars[i]) if(is.numeric(y.factors[i] ) ) {sem[k,2] <- paste("Fy",yvars[i],sep="")} else {sem[k,2] <- y.factors[i]} k <- k +1 } } else { #end of if num.yfactors ==1 for (i in 1:num.y) { for (f in 1:num.yfactors) { lavaan[[num.xfactors +f ]] <- paste(fact[num.xfactors +f], "=~") if( abs(y.factors[i,f]) > cut ) { sem[k,1] <- paste(fact[f+num.xfactors],"->",vars[i+num.xvar],sep="") lavaan[[num.xfactors +f]] <- paste0(lavaan[[num.xfactors +f]], ' + ', yvars[i]) if(is.numeric(y.factors[i,f])) { sem[k,2] <- paste("Fy",f,vars[i+num.xvar],sep="")} else {sem[k,2] <- y.factors[i,f]} k <- k+1 } #end of if } #end of factor } # end of variable loop } #end of else if # } if(errors) { for (i in 1:num.y) { if(lr) {dia.self(var.rect[[i+num.xvar]],side=3) } else {dia.self(var.rect[[i+num.xvar]],side=3)} sem[k,1] <- paste(vars[i+num.xvar],"<->",vars[i+num.xvar],sep="") sem[k,2] <- paste("y",i,"e",sep="") k <- k+1 }} } #end of if.null(ymodel) # if(!is.null(Rx)) {#draw the correlations between the x variables # for (i in 2:num.xvar) { # for (j in 1:(i-1)) { # if((!is.numeric(Rx[i,j] ) && ((Rx[i,j] !="0")||(Rx[j,i] !="0")))|| ((is.numeric(Rx[i,j]) && abs(Rx[i,j]) > cut ))) { # if (lr) {if(abs(i-j) < 2) { dia.curve(from=var.rect[[j]]$left,to=var.rect[[i]]$left, labels = Rx[i,j],scale=-3*(i-j)/num.xvar)} else { dia.curve(from=var.rect[[j]]$left,to=var.rect[[i]]$left, labels = Rx[i,j],scale=-3*(i-j)/num.xvar)} # } else { # if(abs(i-j) < 2) { dia.curve(from=var.rect[[j]]$bottom,to=var.rect[[i]]$bottom, labels = Rx[i,j],scale=-3*(i-j)/num.xvar)} else {dia.curve(from=var.rect[[j]]$bottom,to=var.rect[[i]]$bottom, labels = Rx[i,j],scale=-3*(i-j)/num.xvar)} # } # }} # } # # } # if(!is.null(Ry)) {#draw the correlations between the y variables # for (i in 2:num.yvar) { # for (j in 1:(i-1)) { # if((!is.numeric(Ry[i,j] ) && ((Ry[i,j] !="0")||(Ry[j,i] !="0")))|| ((is.numeric(Ry[i,j]) && abs(Ry[i,j]) > cut ))) { # if (lr) {if(abs(i-j) < 2) { dia.curve(from=var.rect[[j+num.xvar]]$right,to=var.rect[[i+num.xvar]]$right, labels = Ry[i,j],scale=3*(i-j)/num.xvar)} else {dia.curve(from=var.rect[[j+num.xvar]]$right,to=var.rect[[i+num.xvar]]$right, labels = Ry[i,j],scale=3*(i-j)/num.xvar)} # } else { # if(abs(i-j) < 2) { dia.curve(from=var.rect[[j+num.xvar]]$bottom,to=var.rect[[i+num.xvar]]$bottom, labels = Ry[i,j],scale=3*(i-j)/num.xvar)} else {dia.curve(from=var.rect[[j+num.xvar]]$bottom,to=var.rect[[i+num.xvar]]$bottom, labels = Ry[i,j],scale=3*(i-j)/num.xvar)} # } # }} # } # # } Phi <- round(Phi,digits) if(!regression) { if(!is.null(Phi)) {if (!is.matrix(Phi)) { if(!is.null(FALSE)) {Phi <- matrix(c(1,0,Phi,1),ncol=2)} else {Phi <- matrix(c(1,Phi,Phi,1),ncol=2)}} if(num.xfactors>1) {for (i in 2:num.xfactors) { #first do the correlations within the f set for (j in 1:(i-1)) { {if((!is.numeric(Phi[i,j] ) && ((Phi[i,j] !="0")||(Phi[j,i] !="0")))|| ((is.numeric(Phi[i,j]) && abs(Phi[i,j]) > cut ))) { if(Phi[i,j] == Phi[j,i] ) { if(lr) {dia.curve(from=fact.rect[[i]]$right,to=fact.rect[[j]]$right, labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} else { dia.curve(from=fact.rect[[i]]$top,to=fact.rect[[j]]$top, labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} sem[k,1] <- paste(fact[i],"<->",fact[j],sep="") sem[k,2] <- paste("rF",i,"F",j,sep="") lavaan[[num.xfactors +num.yfactors +1]] <- paste(fact[i], "~~", fact[j])} else {#directed arrows if(Phi[i,j] !="0") { if(lr) { if(abs(i-j) < 2) {dia.arrow(from=fact.rect[[j]],to=fact.rect[[i]], labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} else { dia.curved.arrow(from=fact.rect[[j]]$right,to=fact.rect[[i]]$right, labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} } else { if(abs(i-j) < 2) { dia.arrow(from=fact.rect[[j]],to=fact.rect[[i]], labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} else { dia.curved.arrow(from=fact.rect[[j]]$top,to=fact.rect[[i]]$top, labels = Phi[i,j],scale=2*(i-j)/num.xfactors)} } sem[k,1] <- paste(fact[j]," ->",fact[i],sep="") sem[k,2] <- paste("rF",j,"F",i,sep="")} else { if(lr) { if(abs(i-j) < 2) {dia.arrow(from=fact.rect[[i]],to=fact.rect[[j]], labels = Phi[j,i],scale=2*(i-j)/num.xfactors)} else { dia.curved.arrow(from=fact.rect[[i]]$right,to=fact.rect[[j]]$right, labels = Phi[j,i],scale=2*(i-j)/num.xfactors)} } else { if(abs(i-j) < 2) { dia.arrow(from=fact.rect[[i]],to=fact.rect[[j]], labels = Phi[j,i],scale=2*(i-j)/num.xfactors)} else { dia.curved.arrow(from=fact.rect[[i]]$top,to=fact.rect[[j]]$top, labels = Phi[j,i],scale=2*(i-j)/num.xfactors)} } sem[k,1] <- paste(fact[i],"<-",fact[j],sep="") lavaan[[num.xfactors +num.yfactors +k]] <- paste(fact[j], "~", fact[i]) sem[k,2] <- paste("rF",i,"F",j,sep="")} } } else { sem[k,1] <- paste(fact[i],"<->",fact[j],sep="") if (is.numeric(Phi[i,j])) {sem[k,2] <- paste("rF",i,"F",j,sep="")} else {sem[k,2] <- Phi[i,j] } } k <- k + 1} } } } #end of correlations within the fx set if(!is.null(ymodel)) { if(!cancor) { for (i in 1:num.xfactors) { for (j in 1:num.yfactors) { if((!is.numeric(Phi[j+num.xfactors,i] ) && (Phi[j+num.xfactors,i] !="0"))|| ((is.numeric(Phi[j+num.xfactors,i]) && abs(Phi[j+num.xfactors,i]) > cut ))) { dia.arrow(from=fact.rect[[i]],to=fact.rect[[j+num.xfactors]],Phi[j+num.xfactors,i],adj=i %% adj +1) sem[k,1] <- paste(fact[i],"->",fact[j+num.xfactors],sep="") lavaan[[num.xfactors +num.yfactors +k]] <- paste(fact[j+num.xfactors], "~", fact[i]) } else { sem[k,1] <- paste(fact[i],"<->",fact[j+num.xfactors],sep="")} if (is.numeric(Phi[j+num.xfactors,i])) {sem[k,2] <- paste("rX",i,"Y",j,sep="")} else {sem[k,2] <- Phi[j+num.xfactors,i] } k <- k + 1 } }} else { #we just plot the canoncial correlations which are just between the factors for(i in 1:num.xfactors) { for (j in 1:num.yfactors) { if(Phi[j,i] > cut) dia.arrow(to=fact.rect[[i]],from=fact.rect[[j+num.xfactors]],Phi[j,i],both=TRUE) }} } } } } if(num.factors > 0 ) { for(f in 1:num.factors) { sem[k,1] <- paste(fact[f],"<->",fact[f],sep="") sem[k,3] <- "1" k <- k+1 } model=sem[1:(k-1),] class(model) <- "mod" #suggested by John Fox to make the output cleaner lavaan <- unlist(lavaan) lavaan <- noquote(lavaan) result <- list(sem=model,lavaan=lavaan) invisible(result) } } #use the esemDiagram function to draw a canonical correlation diagram cancorDiagram <- function(fit, cut=.1,digits=2, nfX = NULL, nfY =NULL,simple=FALSE,e.size=.1,main="Canonical Correlation",...) { esem <- list() if(is.null(nfX)) nfX <- NCOL(fit$Xmat ) if(is.null(nfY)) nfY <- NCOL(fit$Ymat ) esem$loadsX <- fit$Xmat[,1:nfX,drop=FALSE] esem$loadsY <- fit$Ymat[,1:nfY,drop=FALSE] if(is.complex(fit$cancor)) fit$cancor <- Re(fit$cancor) #don't work with imaginary roots esem$Phi <- diag(fit$cancor) # if(is.complex(esem$Phi)) esem$Phi <- Re(esem$Phi) class(esem) <- class(fit) esemDiagram(esem, cut=cut,digits=digits,simple=simple,e.size=e.size, main=main,...) }