Repeated measures on hundreds of subjects: too many row effects?

[gerverska]

Hello again @JenniNiku and @BertvanderVeen!

I had a question related to the (over?)specification of row effects in models. I have a study in which the same 491 trees had their needle microbiomes surveyed at 2 sites (each with 4 spatially distinct blocks, totaling 8 blocks). Each tree was sourced from one of six seed regions and it's a categorical variable of interest. Certain trees occur within certain blocks at each garden. Sampling year is coded as a numeric variable, and I'm explicitly interested in how fungal communities change over the sampling period, accounting for interactions with site and region. For 491 trees and 3 sampling years, there are a total of 1473 samples. I have 14 fungal taxa within my species matrix, with the total matrix containing around 50 million read counts across these samples and taxa.

This is my current model:

model.gxext <- gllvm(tab[, 2:(ncol(tab) - 2)],
                     meta |> select(year, site, region),
                     formula = ~ year * site + year * region + site * region,
                     family = 'negative.binomial',
                     num.lv = 0,
                     studyDesign = meta |> select(block, tree),
                     row.eff = ~ (1 | block / tree),
                     offset = offset,
                     control.start = list(n.init = 10))

When I had initially attempted formula = ~ year * site * region, I got a warning:

  Determinant of the variance-covariance matix is zero. Please double check your model for e.g. overfitting or lack of convergence. 

I then tried formula = ~ year * site + year * region + site * region, (which I thought was essentially the same as first formula I tried), but this time I didn't get any error and residuals plots looked fine.

However, whenever I re-run the same model without seed on different dates, even the formula = ~ year * site + year * region + site * region model can produce this warning, so it seems like there are "lucky" seeds.

Still troubling, running coefplot() with models that do converge without warning generates confidence intervals that are practically non-existent. Examining summary() output reveals that many coefficients have the same, very low P values. Standard errors are also very low (1x10^-7 or so), but seem to vary "naturally".

Is this some indication that the model is still over-specified but just isn't throwing any errors/warnings, or that the P values just aren't being calculated reliably? I know the 491 tree row effects are very high, but it's not apparent how else I might deal with the repeated measures issue and still examine how the abundances of different fungal taxa change over time. Ultimately, if a model has the capacity to (easily) produce this warning, is warning-free output likely still suspect in certain ways?

Thanks for your help! I can share follow up info, too.

[BertvanderVeen]

What is the scale estimate for the row-effects? model$params$sigma. If that is sufficiently far away from zero I am not worried about the row effects.

It does sound like something fishy is going on with your model, but diagnosing it can be difficult. Start by making sure to set n.init = 3 or so, and set gradient.check = TRUE for some more help with convergence checking.

Have you checked that you have sufficient (non-zero) observations in each of the combination of categories in the covariates, and that a negative-binomial response distribution is necessary?

[gerverska]

Ah, an update: looking at model$params$sigma for when reltol = 1e-12, values seem to increase:

> model_e12$params$sigma

tree:block      block 
0.32202550 0.07414664

Although this output still gives the large gradients (>0.05) warning, is this model still better than the one I got with default reltol, which also coverged but has lower sigma values (see previous reply)?

[gerverska]

The estimates I get from the lower reltol model also seem pretty different (see OTU0002, OTU0005, OTU0008)

# Lower reltol ####
> reltol_e12 <- summary(model_e12)$Coef.tableX

> reltol_e12 |> head()

                Estimate   Std. Error     z value     Pr(>|z|)
year:OTU0001 -0.31085026 4.975308e-05 -6247.85905 0.000000e+00
year:OTU0002 -0.39750498 1.469578e-04 -2704.89112 0.000000e+00
year:OTU0003  0.27268358 7.989876e-05  3412.86388 0.000000e+00
year:OTU0004  0.23563912 7.816223e-05  3014.74393 0.000000e+00
year:OTU0005 -0.00486743 4.847735e-04   -10.04063 1.010312e-23
year:OTU0008 -1.37172573 4.303985e-04 -3187.10597 0.000000e+00

# Higher (default) reltol ####
> reltol_e8 <- summary(model_e8)$Coef.tableX

> reltol_e8 |> head()

                Estimate   Std. Error     z value Pr(>|z|)
year:OTU0001 -0.27848896 4.894943e-05 -5689.32011        0
year:OTU0002 -0.75944263 1.427724e-04 -5319.25269        0
year:OTU0003  0.30817478 8.096538e-05  3806.25366        0
year:OTU0004  0.15084529 7.653858e-05  1970.83989        0
year:OTU0005 -0.07726671 5.182191e-04  -149.10047        0
year:OTU0008 -0.02480134 4.589724e-04   -54.03666        0

But in both cases you can see the low P values...

[BertvanderVeen]

ok, great! With the number of observations, is that the number of non-zero observations in each category?

If you increase reltol and you get a warning that the maximum iterations has reached, please do increase the number of iterations. There is a fair chance the model has not actually converged, though falls warning are also common. Increasing the maximum iterations might not change anything, you at least need to check.

Another thing you can try is changing the optimizer and see if that gives you a better/different model fit. Generally, small -p-values, a symptom of small standard errors, indicate inaccuracy when calculating the Hessian. This usually indicates that one of the blocks in the asymptotic joint precision matrix is singular, which in turn does tend to indicate a problem in the model fit.

As a last resort, after trying changing the number of iterations and optimizer, you can set method = "EVA" and/or method = "LA". This will fit the model with a different parameterisation, and might help with some of the issues you are encountering.

[gerverska]

Thanks Bert! I'll try your recommendations and report back.

And sorry if I'm not being clear about the study design. I've attached a figure that hopefully communicates the sampling design better. The numbers I gave earlier were the number of trees for every year, site, and region combination, and I neglected to account for block.

So just to refresh on the current model properties:

formula = ~ year * site + year * region + site * region,
row.eff = ~ (1 | block / tree)

Do you think the unbalanced nature of the blocks is creating some of the convergence/large gradient issues?

[BertvanderVeen]

No, that would surprise me.