How to specify curved lines for covariances and custom placement of residual variance when using a custom layout in semPaths

[libbybenson3]

Thank you for the awesome work that went into the creation of the semPlot package and semPaths() - I use it all the time for my own work and for teaching purposes :)

I am assisting a professor with creating teaching materials for a longitudinal structural equation modeling class, and have ran into a few snags.

For example, I am working on cleaning up a cross-lagged model plot. The first two steps I took were to specify a custom layout for the latent and manifest variables, and then to clean up the paths with overlapping values.

Question 1: I then wanted to add curved lines for the autocorrelated errors using the curvePivot argument, but it seems that this argument does not work when a custom layout is specified.

Question 2: I also wanted to try to move some of the residual variance double headed arrows, but couldn't figure out how to do this based on the documentation.

Below, I have included a data generation method, the model I am fitting, and the different plot iterations I went through. Any advice on how to solve question 1 and/or question 2 would be greatly appreciated. Thank you! Libby

# create data
mycormat <- matrix(c(1.00,  0.42, 0.10,  0.12,  0.50,  0.14, -0.07, -0.10,
                     0.42,  1.00, 0.09,  0.08,  0.28,  0.41, -0.07, -0.16,
                     0.10,  0.09, 1.00,  0.38,  0.03,  0.00,  0.43,  0.10,
                     0.12,  0.08, 0.38,  1.00,  0.09, -0.03,  0.15,  0.36,
                     0.50,  0.28, 0.03,  0.09,  1.00,  0.26, -0.04,  0.01,
                     0.14,  0.41, 0.00, -0.03,  0.26,  1.00,  0.01, -0.05,
                     -0.07, -0.07, 0.43,  0.15, -0.04,  0.01,  1.00,  0.20,
                     -0.10, -0.16, 0.10,  0.36,  0.01, -0.05,  0.20,  1.00),
                   nrow = 8, ncol = 8, 
                   dimnames = list (c("y11", "y21", "y31", "y41",
                                      "y12", "y22", "y32", "y42")))

nvars <- 8
nobs <- 500
L <- chol(mycormat)
r = t(L) %*% matrix(rnorm(nvars*nobs), nrow=nvars, ncol=nobs)
r = t(r)
rdata = as.data.frame(r)
names(rdata) = c("y11", "y21", "y31", "y41","y12", "y22", "y32", "y42")

# Specify model 
  model <- 
  '
  # measurement model
    # factor 1, occasion 1
      F11 =~ 1*y11 + y21 
    # factor 2, occasion 1
      F21 =~ 1*y31 + y41
    # factor 1, occasion 2
      F12 =~ 1*y12 + y22
    # factor 2, occasion 2
      F22 =~ 1*y32 + y42
  # allow for auto-correlated errors for same items across occasions
    y11 ~~ y12
    y21 ~~ y22
    y31 ~~ y32
    y41 ~~ y42
  # specify the cross-lagged relations
    F12 ~ F11 + F21
    F22 ~ F11 + F21
  '

# Fit model 
  model_fit <- cfa(model, data=rdata)

# Look at results
  summary(model_fit, fit.measures=TRUE)

# Plot (using semPaths defaults)
   semPaths(model_fit, what = "model", "est",
         sizeLat = 6, sizeMan = 6, sizeInt = 5,
         edge.label.cex = 1)

# New layout I would like to use (to better represent time)
   newlayout <- matrix(c(-.5, .75, # y11
                     -0.25, 0.75, # y21
                     -0.5, -0.75, # y31
                     -0.25, -0.75, # y41
                     0, .75, # y12
                     0.25, 0.75, # y22
                     0, -.75, # y32
                     .25, -.75, # y42
                     -.25, .25, # F11
                     -.25, -.25, # F21
                     0, .25, # F12
                     0, -.25), # F22
                   ncol=2, byrow=TRUE)

# Plot with new layout
   semPaths(model_fit, what = "model", "est",
         sizeLat = 6, sizeMan = 6, sizeInt = 5,
         edge.label.cex = 1, layout = newlayout)

# Now to fix the paths with overlapping numbers
   newedges <- c(.6, .6, .6, .6, # manifest loadings on latent variables
                 .6, .6, .6, .6, # manifest loadings on latent variables
                 .25, .25, .25, .25, # autocorrelations among measurement errors
                 .5, .25, .25, .5, # cross-lag
                 .5, .5, .5, .5, # manifest variance (occasion 1)
                 .5,.5,.5,.5, # manifest variance (occasion 2)
                 .5,.5,.5,.5, # latent factor variance
                 .5,.5,.25,.25, 
                 .25, .25, .5,.5)  

# Plot with new layout and new edges
   semPaths(model_fit, what = "model", "est",
         sizeLat = 6, sizeMan = 6, sizeInt = 5,
         edge.label.cex = 1, layout = newlayout,
         edge.label.position = newedges)

# autocorrelations among measurement errors are hard to see/ misleading - try curved lines
    semPaths(model_fit, what = "model", "est",
         sizeLat = 6, sizeMan = 6, sizeInt = 5,
         edge.label.cex = 1, layout = newlayout,
         edge.label.position = newedges, 
         curvePivot = TRUE)

    # For some reason curvePivot does not work with custom layout set
    # Specifically, I would like to figure out how to have curved lines for:
      # y11 ~~ y12
      # y21 ~~ y22
      # y31 ~~ y32
      # y41 ~~ y42

    # I also see that for y41, the measurement error variance is below the manifest variable, whereas for y21 it is to the right
    # Is it also possible to customize this feature? I like how it looks for y41 much better than y21

[SachaEpskamp]

You can tweak things after plotting manually using qgraph. Do the codes below do what you want?

Graph <- semPaths(model_fit, what = "model", "est",
         sizeLat = 6, sizeMan = 6, sizeInt = 5,
         edge.label.cex = 1, layout = newlayout,
         edge.label.position = newedges, 
         curvePivot = TRUE)

# Extract some info:
Edgelist <- as.data.frame(Graph$Edgelist)
Labels <- Graph$graphAttributes$Nodes$labels

# Make curve object:
Curve <- rep(0, nrow(Edgelist))

# Make some edges curved:
Curve[Edgelist$from == which(Labels == "y12") & Edgelist$to == which(Labels == "y11")] <- Curve[Edgelist$from == which(Labels == "y11") & Edgelist$to == which(Labels == "y12")]  <- 1

Curve[Edgelist$from == which(Labels == "y21") & Edgelist$to == which(Labels == "y22")] <- Curve[Edgelist$from == which(Labels == "y22") & Edgelist$to == which(Labels == "y21")]  <- 1

Curve[Edgelist$from == which(Labels == "y31") & Edgelist$to == which(Labels == "y32")] <- Curve[Edgelist$from == which(Labels == "y32") & Edgelist$to == which(Labels == "y31")]  <- -1

Curve[Edgelist$from == which(Labels == "y41") & Edgelist$to == which(Labels == "y42")] <- Curve[Edgelist$from == which(Labels == "y42") & Edgelist$to == which(Labels == "y41")]  <- -1

# Now the loopRotation object to fix the rotation:
loopRotation <- Graph$graphAttributes$Nodes$loopRotation
loopRotation[Labels == "y31"] <- pi
loopRotation[Labels == "y42"] <- pi
loopRotation[Labels == "y11"] <- 0
loopRotation[Labels == "y21"] <- 0
loopRotation[Labels == "y22"] <- 0

library(qgraph)
qgraph(Graph, curve = Curve, loopRotation = loopRotation)

[libbybenson3]

Hi Sacha, Yes! Your code answers/ solves both of my questions - thank you :)