8
$\begingroup$

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).

Random effects estimates

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 ?
$\endgroup$
10
  • 2
    $\begingroup$ (+1) Welcome, this is a well-formulated question! About your model: have you accounted for the shelterbelts in there somehow? I would think that would merit another random effect (possibly nested, depending on your naming convention). $\endgroup$
    – mkt
    Jul 6, 2022 at 18:45
  • 3
    $\begingroup$ Personally, I would keep the random effects in there based on our a priori expectations of spatial autocorrelation, and not remove them unless there were problems with the fitting. I don't see much downside of having some weak random effects in the model, even if we expect them to be much stronger. $\endgroup$
    – mkt
    Jul 6, 2022 at 18:48
  • 3
    $\begingroup$ The random effects in your model are random intercepts since E{size_transect} = re_transect + betas * various_fixed_effects + beta * age. Your plot is a (forest) plot of the re_transect estimates. These random intercepts might be close to zero but are not all equal, so their estimated variance is not zero. Ignoring the structure of the data (trees are nested within transects) means underestimating the residual variance. $\endgroup$
    – dipetkov
    Jul 6, 2022 at 23:20
  • 2
  • 2
    $\begingroup$ To oversimplify a bit, the random effects are one way of accounting for non-independence between points. So if there is spatial autocorrelation at multiple scales (farm/shelterbelt/transect), I would begin by including all of those as random intercepts at least. Random slopes may be appropriate, but they also make the models more complex and might lead to singular fits. The naming convention you mention helps you decide between whether to specify the random effects as crossed or nested (the transect would be crossed, in your case). $\endgroup$
    – mkt
    Jul 7, 2022 at 6:25

1 Answer 1

6
$\begingroup$

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).

$\endgroup$
1
  • 2
    $\begingroup$ It might also be helpful to mention that the standard errors of the random effects estimates depend on the group sizes: the more trees from a transect, the smaller standard error its re_transect has. So we could guess that transect 1211 has more trees sampled than transect 914. $\endgroup$
    – dipetkov
    Jul 10, 2022 at 8:12

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.