Are mixed models necessary if random effects estimates are close to zero?

Question

I am currently using (general) linear mixed models as a mean to avoid pseudoreplication and control for measures made on various individuals in the same location (transect). I'm now starting to think using mixed models is not necessary and would really appreciate some guidance by more knowledgeable folks.

For clarity: I measured survival and size of several trees nested in transects along multiple shelterbelts (think ~1500 trees within ~250 transects within ~40 shelterbelts). I am trying to assess what factors influence trees survival and size in agroforestry systems. A single transect has several corresponding measures (ie soil texture, width, etc.), which are then shared by every single tree measured in the same transect. The data is then structured as a list of trees, with every ~5-10 neighboring individuals having different response data but otherwise sharing the exact same data structure, which is inherited from their "home transect".

Since I felt this was problematic, I use mixed models where I treat transect ID as a random effect, such as this model: (the following code is in R)

size_mod = lmer(size ~ various_fixed_effects + 
           age + (1|transect_id), data=trees)

After fitting a model, I plotted the random effects using sjPlot::plot_model(size_mod, type="re")

I then saw that random effects estimates were close to zero, with every single std. error bar crossing zero (as seen below).

I then saw creutzml's comment on an unrelated post about singular fits, which made me question the utility to use mixed models at all, hence this post.

So, my questions are;

1. would the random effects estimates be considered different from zero ?
2. are mixed models at all necessary if random effects coefficients (odds-ratio/estimates) are zero ? In my case, they ultimately don't seem to have an impact on the models so I'm really not sure. But then,
3. wouldn't using non-mixed models be considered pseudoreplication ?

Answer 1

To summarize the useful information provided in comments: you don't necessarily have to use a mixed model, but you do have to take the lack of independence among the measurements into account in some way. Robust standard errors and generalized estimating equations (GEE), as discussed by Mc Neish et al. (linked from a comment by @mkt), and generalized least squares are alternatives in many situations. Your situation, with trees nested within transects nested within shelterbelts, would seem best represented by a multi-level model.

The inability to distinguish any particular random-effect estimate from 0 in a mixed model isn't important. More important is the variance among those random-effect values, which are modeled as coming from a Gaussian distribution. For example, consider this adaptation of a simple example from the lmer help page:

library(lme4)
(mod <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy))
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (1 | Subject)
#    Data: sleepstudy
# REML criterion at convergence: 1786.465
# Random effects:
#  Groups   Name        Std.Dev.
#  Subject  (Intercept) 37.12   
#  Residual             30.99   
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
# (Intercept)         Days  
#      251.41        10.47 

Notice that the random effects are reported first, with a standard deviation for the random intercept. In this example, that standard deviation is about the same magnitude as the residual standard deviation. The variance (square of the standard deviation) among those random-effect intercept values is what the model primarily estimates. That's what you should primarily pay attention to, in terms of "random effects" being "different from zero" (Question 1). Yes, you can extract the individual random-intercept estimate from the model, but the fact that you can't distinguish any individual values from 0 doesn't mean that the random effect as a whole is unimportant (Question 2).

It's possible to have such low dispersion among the individual random-effect values that your model can't reliably estimate their variance and you end up with a singular model fit. That does not, however, appear to be the case with your data. As described on this page, linked from a comment by @dipetkov, that's typically seen when you try to fit more combinations of random effects than your data allow. How to proceed in such a case requires evaluation of what in particular is leading to the problem. The lack of independence among correlated observations still needs to be taken into account. Otherwise you do suffer from "pseudoreplication" (Question 3).