I have applied a piecewise structural equation model to my data using generalised linear models. I have used the glm() R function along with psem in the "piecewiseSEM" R package.

It has worked, and I would Like to present the results, but because I am self taught I am full of self doubt and I would appreciate if someone could look at my code and output and see if I have made any glaring omissions or mistakes.

GLMs have been used because some of the data is not normally distributed and also because I expect/anticipate that some of the responses are non-linear. I have used a piecewise SEM because I don't know whether some of my variables are predictors or responses, and theoretically they might go either way. Further, I sense it is a more mathematically sound approach than lots of linked regressions. I hope this sounds logical.

The models inside the SEM make use of 11 variables, none of which I transformed before inputting them. They are pH, three vegetation variables which are measured in %, one soil variable measured in %, and three further water quality variables with continuous and unbounded scales. The other 3 variables are species abundances, which are bounded proportions (but they are non-integers).

The GLMs with continuous variables as response variables use Gaussian family distributions whilst the GLMs with the species abudances use quasi-poisson families because they are counts but non-integers. Is this correct?

I have used the following code:

vegspec1 <- glm(specabund1 ~ ph + veg1 + soil1 + veg2*veg3, data = mydata,quasipoisson())
vegspec2 <- glm(specabund2 ~ veg1 + soil1, data = mydata,quasipoisson())
water1 <- glm(water1 ~ specabund1 + specabund2 + specabund3, data = mydata)
water2 <- glm(water2 ~ specabund2 + specabund3, data = mydata)
water3 <- glm(water3 ~ specabund2, data = mydata)

For the GLMs. I have then combined this in a piecewise SEM:

sem.8 <- psem(vegspec1, vegspec2, water1, water2, water3)
## summary - non-linearities detected in here

Because I have used GLMs in the SEM, this throws up the following error:

Error: Non-linearities detected in the basis set where P-values are not symmetrical. This can bias the outcome of the tests of directed separation.

Offending independence claims: 
 Saprotroph.Symbiotroph <- Symbiotroph *OR* Saprotroph.Symbiotroph -> Symbiotroph 
Option 1: Specify directionality using argument 'direction = c()' in 'summary'.
Option 2: Remove path from the basis set by specifying as a correlated error using '%~~%' in 'psem'.
Option 3 (recommended): Use argument 'conserve = TRUE' in 'summary' to compute both tests, and return the most conservative P-value.

I chose option 3, because my hypothesis does not support me in specifying a path. Am I right in doing this? I then run:

summary(sem.8, conserve=TRUE)

This gives me a functioning model with results I am very excited about. In the tests of directed seperation, I get a Fisher's C = 71.321 with P-value = 0.057 and on 54 degrees of freedom. This is a bit too close to a significant p value for me, and the model has returned some variables that are not independent. The AIC of this model 115.321.

I consequently run:

summary(update(sem.8, water2 %~~% water3, water1 %~~%water2, water1 %~~%w water3, specabund1 %~~% specabund3),standardize ="scale",conserve=TRUE)

This new model has a Fisher's C = 43.903 with P-value = 0.56 and on 46 degrees of freedom. The AIC is 87.903 so clearly a better model.

The rest I can confidently interpret: there are significant relationships identified, and at the end of the output the total R2 of each variable is listed, and consequently I can load this into a path diagram.

Can anybody answer my specific questions, and is there any obvious errors or problems that I have missed?

  • $\begingroup$ Did you come closer in finding a solution for estimating a ~generalized structural equation model? $\endgroup$
    – Johan
    Commented Jun 14, 2023 at 8:34

1 Answer 1


One thing I'll note - as mentioned by Lefcheck here, "A key point to be made is that the piecewise approach does not absolve the user of all assumptions associated with the statistical tests."

I interpret this as "each of your component GLMs should also be checked for goodness of fit." For each of your models (vegspec1, vegspec2, water1, water2, water3) I would recommend all of the standard checks: R-squared, delta AIC, overdispersion, residuals vs. fitted plot, using DHARMa to simulate residuals, etc. to make sure that they fit well.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.