Fourth corner analysis with fuzzy coded data

[joopcoolen]

Hello. I love the package, thanks for developing and maintaining it.
I am working with functional trait data, which is fuzzy coded, meaning each trait has multiple colums for the different modalities of the trait. Within a trait the sum of the modalities is 1. A trait can have a value between 0 and 1 in multiple modalities.
I would like to include these data in a fourth corner analysis as outlined here: https://cran.r-project.org/web/packages/gllvm/vignettes/vignette1.html#:~:text=The%20fourth%20corner%20model%20can,random%20slopes%20in%20the%20model.
If I, however, include the fuzzy coded data as traits, I get the following error: Standard errors for parameters could not be calculated, due to singular fit.
I am assuming this is caused by the 0 and 1 in the trait data as without the traits in the model, there is no error. I could not find any information online on how to do this correctly.

Can you please suggest an approach to include fuzzy coded trait data?
Thanks

[BertvanderVeen]

Thanks for posting! I do not see why fuzzy coding would cause any issues here. I'm happy to be proven wrong, but I see no reason why it should be related.

The error you receive can occur in any type of model fitted in the gllvm package, and is often related to the exact model implemented, its convergence, and perhaps the amount of information in your data.

Could you provide some more information on the module you have fitted, and the amount of data that you have available?

[joopcoolen]

Hi Bert, thanks for replying.
This is the model I have fitted:

M1.zinb <- gllvm(y = specdat,
                 X = covdat,
                 TR = traitdat,
                 studyDesign = ExtraVariables,
                 row.eff = ~ (1 | fDate) + (1 | fcast) + (1 | fDepth) + (1 | fsample),
                 formula = ~ (farea + distco.s) + farea : hab.rock ,
                 num.lv = 1,
                 family = "ZINB",
                 control.start = list(n.init = 5,
                                      jitter.var = 0.1),
                 starting.val = "random")

Specdat is data on 40 species (all occuring >10 times in a samples), with 632 samples.
covdat contains farea which is factor data with 4 levels, and distco.s which is standardised continuous data.
traitdat is a single trait (for now) with 6 modalities. I also tried standardising this but this gave the same error
extravariables are all factor data

head(specdat[1:10,1:5]), only first 5 species

Ammodytidae Callionymus lyra Cetorhinus maximus Chirolophis ascanii Ciliata mustela
s1 16110 0 0 0 0
s2 6379 0 0 0 0
s3 32722 0 0 0 0
s4 0 0 0 0 0
s5 0 0 0 0 0
s6 2500 0 0 0 0

head(covdat):

farea distco.s
s1 4_high 0.05872815
s2 2_mid -0.50434692
s3 2_mid -0.36708039
s4 2_mid -0.36708039
s5 1_low 0.51033520
s6 1_low 0.51033520

head(traitdat)

hab.sand hab.gravel hab.rock hab.migrate hab.seagrass hab.shell
Ammodytidae 1 0.0 0.0 0 0 0.00
Callionymus lyra 1 0.0 0.0 0 0 0.00
Cetorhinus maximus 0 0.0 0.0 1 0 0.00
Chirolophis ascanii 0 0.0 1.0 0 0 0.00
Ciliata mustela 0 0.3 0.3 0 0 0.34
Ciliata septentrionalis 1 0.0 0.0 0 0 0.00

head(ExtraVariables)

fDate fcast fDepth fsample
s1 8102023 Cast109 Surface cruise1_41_8102023_Cast109_high2_Surface
s2 9102023 Cast141 Surface cruise1_90_9102023_Cast141_mid1_Surface
s3 4102023 Cast23 Bottom cruise1_23_4102023_Cast23_mid5_Bottom
s4 4102023 Cast23 Surface cruise1_20_4102023_Cast23_mid5_Surface
s5 4102023 Cast27 Bottom cruise1_38_4102023_Cast27_low6_Bottom
s6 4102023 Cast27 Surface cruise1_44_4102023_Cast27_low6_Surface

specdat is strongly zero inflated, so I applied a ZINB

[BertvanderVeen]

Nothing obvious jumps out, except perhaps that you seem to have both fixed effects species intercepts included, as well as random. The former is the default, and can be collected to a single fixed intercept by using the beta0com argument in the gllvm function. The random intercepts are included due to the lack of a 0+ notation in the random effects part of formula.

Having both included can cause the variance of the random effect to be estimated near zero, which could explain the error you're witnessing. For now, I suggest you try adding beta0com = TRUE and refit the model.

If that's not the cause, we can reiterate!

[joopcoolen]

Thanks for the help so far, I am already happy to read that the fuzzy coding in itself is probably not the cause of my issue.
I added beta0com = TRUE to the model as below, but still the singular fit response is given.
Could you refer me to documentation on single fixed intercept vs fixed effects species intercepts in gllvm context as I cant find any explanation on this.
As I considered that maybe I was asking too much from the data with the model, I also tried the same model but with only row.eff = ~ (1 | fsample) and when that gave the singular fit response again, I also tried a model with row.eff = ~ (1 | fsample) as well as formula = ~ farea : hab.rock only, but again this gave the singular fit response.

M1.zinb <- gllvm(y = specdat,
                 X = covdat,
                 TR = traitdat,
                 studyDesign = ExtraVariables,
                 row.eff = ~ (1 | fDate) + (1 | fcast) + (1 | fDepth) + (1 | fsample),
                 formula = ~ (farea + distco.s) + farea : hab.rock ,
                 num.lv = 1,
                 family = "ZINB",
                 beta0com = TRUE,
                 control.start = list(n.init = 5,
                                      jitter.var = 0.1),
                 starting.val = "random")

[AlainZuur]

Joop...I'm not into fourth corner models....but what happens if you use an ordinary negative binomial distribution instead of a ZINB? I find that a (numerically) very challenging distribution in other packages. Plenty of zeros doesn't mean that you have to do a ZINB.

Alain

[BertvanderVeen]

There is no (extensive) documentation regarding the intercepts yet, beyond that in ?gllvm, but I can try to answer any questions you might have.

"singular fit" is a generic error that flags an issue in the calculation of standard errors. Occasionally, this can be caused by a bug in the software, but usually it is caused by a parameter estimate that is on the boundary of the valid space with very small variance (such as for the variances of the random effects). You might be able to spot those via the summary method of the model object.

I can try to look further into a solution for you tomorrow. Unfortunately I am traveling today and do not have access to a computer right now, which makes that a tad difficult.

Have you tried fitting your model with starting.val = "res"?

[joopcoolen]

Hi Alain, changing the distribution indeed works, great! Now the model runs as it should, also with when I include all the modalities.
Bert and Alain thank you both for the help, much appreciated!

[AlainZuur]

Now simulate data from your model and assess whether the NB gllvm can cope with the % of zeros that you have in the observed data.

The model that you originally used may have had various components to model the zeros. For example, if one level of your random effects has only zeros, then a relative large negative random effect can model those zeros. But so can the pi of the ZINB. And so can the theta of the variance term. Hence...plenty of parameters in your model that can all go for the same zeros.....and therefore you get numerical estimation issues. That is not a gllvm-specific problem; you will also get that with glmmTMB or INLA.

Alain

[joopcoolen]

I did that using testZeroInflation() from DHARMa, and all species in the data are now zero inflated, while with ZINB (without the traits in the model) this was not the case.
So ZINB is nog working and negative binomial causes zero inflation. Any suggestions on what to do now?

[AlainZuur]

try and figure out which species is causing it? Alain

…

On 18 Dec 2024 at 12:49 +0000, joopcoolen ***@***.***>, wrote: I did that using testZeroInflation() from DHARMa, and all species in the data are now zero inflated, while with ZINB (without the traits in the model) this was not the case. So ZINB is nog working and negative binomial causes zero inflation. Any suggestions on what to do now? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you commented.Message ID: ***@***.***>

[BertvanderVeen]

Joop, I tracked your issue down to a bug in the calculation of the standard errors. I have now fixed this, so that if you reinstall the package from the development branch, it should fit without error. Apologies for the inconvenience, and thanks for getting in touch about this.