Random effects induce heteroskedasticity and correlation in "error terms" in the model
A useful way to look at random effects models is to use some mathematical manipulation to mash them back onto the form of a traditional model without a random effect term. This can be accomplished by absorbing the random effects terms into the error terms of a traditional fixed effect model. We can then examine the properties of the new error term in the traditional fixed effect model form, to see what happens with this variation.
To illustrate this technique, consider a Gaussian linear random effects model with observations $Y_{i,j}$ taken over categories $j=1,...,k$. The model can be written as:
$$Y_{i,j} = \beta_0 + \beta_j + u_j + \varepsilon_{i,j}
\quad \quad \quad
u_j \sim \text{N}(0,\sigma_j^2)
\quad \quad \quad
\varepsilon_{i,j} \sim \text{N}(0,\sigma^2),$$
where the random effects terms and error terms are mutually independent. Since both the error terms and the random effect terms are random variables in the model, we can combine them into a single alternative error term and write the model in a form that does not show a separate random effect. Specifically, we write the model in the traditional fixed effect form:
$$Y_{i,j} = \beta_0 + \beta_j + \eta_{i,j},$$
where we have defined the quantities:
$$\eta_{i,j} = z_j + \varepsilon_{i,j}
\quad \quad \quad \quad \quad
\rho_j = \frac{\sigma_j^2}{\sigma_j^2+\sigma^2}
\quad \quad \quad \quad \quad
\sigma_{j*}^2 = \sigma_j^2+\sigma^2.$$
In this latter form, the error terms $\eta_{i,j}$ are still (jointly) normally distributed with zero mean, but they are no longer homoskedastic and uncorrelated --- they have covariance values:
$$\mathbb{C}(\eta_{i,j},\eta_{i',j'})
= \mathbb{C}(z_j + \varepsilon_{i,j}, z_{j'} + \varepsilon_{i',j'})
= \begin{cases}
\sigma_{j*}^2 & & \text{if } i = i' \text{ and } j = j', \\[6pt]
\rho_j \sigma_{j*}^2 & & \text{if } i \neq i' \text{ and } j = j', \\[6pt]
0 & & \text{otherwise}. \\[6pt]
\end{cases}$$
This form means that there is heteroskedasticity in the model (i.e., with variance $\sigma_{j*}^2$ for observations in category $j$) and the errors within a group are positively correlated (with correlation coefficient $\rho_{j}$).
Should you use a random effects model? As you can see from the above, a Gaussian linear random effects model using categorical predictors is equivalent to a traditional Gaussian linear regression model using categorical predictors, where the latter has heterosedasticity across the categories of observations and positively correlated errors within each category. Consequently, your choice of whether or not to include random effects terms can be framed equivalently as a choice of whether or not to generalise the behaviour of the error terms in the model to allow heteroskedasticity across categories and correlation of error terms within categories.
In the study you are conducting, your forest patches are the categories and your plots within these forest patches are your individual observations. In this case, including a random effect in your model (taken at the level of the forest patches) is equivalent to allowing heteroskedasticity across different forest patches and also having positive correlation of the mammal abundance of different plots within each forest patch. In order to decide whether or not this is appropriate, you merely need to ask yourself if this type of heteroskedasticity/correlation might plausibly occur in this case.
You have also noted that your forest patches were not randomly sampled. This is not necessarily a problem for the random effects model, and it is no more of a problem than to a fixed effects model. In assessing your sampling method you should consider whether your choice of sites was influenced by any of the variables under study, and consider the types of biases this could induce. However, there is nothing inherent in the random effects model (as opposed to the fixed effects model) that presents an analytical difference here.
I would suggest that in the present case this kind of heteroskedasticity/correlation might plausibly occur, owing to the closeness of the plots within a patch and the possible movement behaviour of the mammals under study. Different plots within the same forest patch might plausibly have positively correlated mammal abundance due to the movement of mammals around a forest patch, breeding of mammals from amongst nearby plots within a forest patch, and common conditions/threats to mammals within the same forest patch. You have noted that you tried both the random effects and fixed effects models and conducted a likelihood-ratio test, finding that there was no significant evidence of the presence of random effects. That is perfectly fine, and it is one way to conduct your analysis --- you have an initially plausible model form and then you find that it does not operate better than a simpler model form.
