piecewiseSEM icon indicating copy to clipboard operation
piecewiseSEM copied to clipboard

Model convergence and scale()

Open luroy opened this issue 10 months ago • 0 comments

Dear all, I have just started to use psem package to perform Piecewise Structural Equation Modeling. I therefore apologize in advance if my following comment is due to my unfamiliarity with this package. I'm trying to run the following psem:

SEM<-psem(bm_lc_mod1.1<-lm(Tmi ~ 1, data = data),
                 bm_lc_mod1.3<-lm(Tma ~ p1, data = data),
                 bm_lcpwc_mod1.1b<-glm(Rich ~ p1 + Pra1 + Tmi + Tma, family = "poisson", data = data),
                 bm_lcpwc_mod1.2b<-glm.nb(AB ~ Autoc_AB + p2 + Pra2 + p3, data = data),
                 bm_lcpwc_mod1.3<-lm(Co1 ~ p1 + Pra1, data = data),
                 bm_lcpwc_mod1.4b<-lm(Co2 ~ Autoc_Co2 + Pra2, data = data),
                 bm_c_mod1.1q<-glmer(Tx_V~AB+Co1+(1|obs_V), weights=N_V, data = data, family="binomial"),
                 bm_c_mod1.2q<-glmer(Tx_P~Rich+Autoc_Tx_P+(1|obs_P),  weights=N_P, data = data, family="binomial")
)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
2: Some predictor variables are on very different scales: consider rescaling 
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

When I look at the summary, I get the following error

summary(SEM)
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: 
 Co1 <- Tmi *OR* Co1 -> Tmi 
 
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.

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: 
 Co1 <- Tmi *OR* Co1 -> Tmi 
 
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.

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

I've chosen the first option, and when I call up the summary, I get the following warning message:

summary(SEM,  direction = c("Co1 <- Tmi"))

Structural Equation Model of SEM 

Call:
  Tmi ~ 1
  Tma ~ p1
  Rich ~ p1 + Pra1 + Tmi + Tma
  AB ~ Autoc_AB + p2 + Pra2 + p3
  Co1 ~ p1 + Pra1
  Co2 ~ Autoc_Co2 + Pra2
  Tx_V ~ AB + Co1
  Tx_P ~ Rich + Autoc_Tx_P

    AIC
[..., not shown]

---
Tests of directed separation:
[..., not shown]

---
Global goodness-of-fit:
Chi-Squared = NA with P-value = NA and on 77 degrees of freedom
Fisher's C = 400.808 with P-value = 0 and on 146 degrees of freedom

---
Coefficients:

[..., not shown]
---
Individual R-squared:

[..., not shown]

Warning message:
Check model convergence: log-likelihood estimates lead to negative Chi-squared! 

This convergence error (that leads to Chi-Squared = NA and associated p-value = NA) is linked to bm_c_mod1.1q and bm_c_mod1.2q, which are binomial regressions with a observation-level random effect (ORLE) to correct overdispersion.

These two models, when run independently of the psem, give the following errors:

bm_c_mod1.1q<-glmer(Tx_V~AB+Co1+(1|obs_V), weights=N_V, data = data, family="binomial")
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
bm_c_mod1.2q<-glmer(Tx_P~Rich+Autoc_Tx_P+(1|obs_P),  weights=N_P, data = data, family="binomial")
Warning messages:
1: Some predictor variables are on very different scales: consider rescaling 
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
 Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?

These errors disappear when I scale the fixed effects (using center=FALSE to obtain the same p-values).

bm_c_mod1.1q<-glmer(Tx_V~as.numeric(scale(AB, center=FALSE))+as.numeric(scale(Co1, center=FALSE))+(1|obs_V), weights=N_V, data = data, family="binomial")
bm_c_mod1.2q<-glmer(Tx_P~as.numeric(scale(Rich, center=FALSE))+as.numeric(scale(Autoc_Tx_P, center=FALSE))+(1|obs_P),  weights=N_P, data = data, family="binomial")

Oddly enough, if I choose the third option (recommended), I get the similar results, but the Global goodness-of-fit is different and the warning does not appear.

Chi-Squared = 0 with P-value = 1 and on 78 degrees of freedom
Fisher's C = 211.94 with P-value = 0 and on 134 degrees of freedom

So I tried - first - to scale(... center=FALSE) directly in the psem :

SEM_b<-psem(bm_lc_mod1.1<-lm(Tmi ~ 1, data = data),
                 bm_lc_mod1.3<-lm(Tma ~ p1, data = data),
                 bm_lcpwc_mod1.1b<-glm(Rich ~ p1 + Pra1 + Tmi + Tma, family = "poisson", data = data),
                 bm_lcpwc_mod1.2b<-glm.nb(AB ~ Autoc_AB + p2 + Pra2 + p3, data = data),
                 bm_lcpwc_mod1.3<-lm(Co1 ~ p1 + Pra1, data = data),
                 bm_lcpwc_mod1.4b<-lm(Co2 ~ Autoc_Co2 + Pra2, data = data),
                 bm_c_mod1.1q<-glmer(Tx_V~as.numeric(scale(AB, center=FALSE))+as.numeric(scale(Co1, center=FALSE))+(1|obs_V), weights=N_V, data = data, family="binomial"),
                 bm_c_mod1.2q<-glmer(Tx_P~as.numeric(scale(Rich, center=FALSE))+as.numeric(scale(Autoc_Tx_P, center=FALSE))+(1|obs_P),  weights=N_P, data = data, family="binomial")
)

But when I call the summary (specifying the direction of the error, which is the same as above), I get:

summary(SEM_b,  direction = c("Co1 <- Tmi"))
  |======================================================================                                                                                  |  46%Error in data.frame(Independ.Claim = paste(b[[i]][2], "~", rhs), ct[,  : 
  l  arguments imply differing number of rows: 1, 0

So I tried, secondly, to continue with the first SEM (SEM_R4_ENS, despite the convergence (and NA) problem). The summary shows 8 significant separation tests, which I have integrated in the form of correlated errors.

SEM_2<-update(SEM_R4_ENS, ... %~~% ..., 
                                   ... %~~% ..., 
                                   ... %~~% ...,
                                   ... %~~% ..., 
                                   ... %~~% ..., 
                                   ... %~~% ..., 
                                   ... %~~%  ..., 
                                   ... %~~% ...)

Here, when running summary(SEM_2, direction = c("Co1 <- Tmi")), there is no significant tests of directed separation, and Global goodness-of-fit is Chi-Squared = NA with P-value = NA and on 69 degrees of freedom Fisher's C = 94.43 with P-value = 0.992 and on 130 degrees of freedom

And when running summary(SEM_2, conserve = TRUE), there is also no ignificant tests of directed separation, and Global goodness-of-fit is Chi-Squared = 0 with P-value = 1 and on 70 degrees of freedom Fisher's C = 90.678 with P-value = 0.979 and on 120 degrees of freedom.

Since both summaries indicate that Fisher's C is non significant, can I consider my psem() to be complete, despite the convergence problems identified above?

Thank you for your attention to this matter,

Léa

luroy avatar Apr 23 '24 14:04 luroy