Accounting for uncertainty in a mixed-effects regression I have calculated an effect size along a dataset of experiments distributed worldwide through a mixed-effects meta-regression. The effect size in the dataset depends on climate (y ~ temperature + precipitation):
meta <- metafor::rma(es, var, data=df, mods= ~ 1 + precip + temp))
As I have data for temperature and precipitation for virtually all points on Earth, I have upscaled this effect and standard error (SE) globally, that is, applying the model to predict the effect for all points on Earth, thus creating a gridded map:
pred <- predict(meta, newmods = cbind(s.df$precipitation, s.df$temperature))

This upscaled map of SE thus represent the meta-analysis derived uncertainty. What this SE does not reflect is the uncertainty derived from the fact that the number of experiments in the dataset used to run the meta-regression is obviously limited, yet I predict the effect globally, even in areas not covered by experiments. In particular, it is very obvious in the dataset that tropical ecosystems are not well represented, with just a couple of experiments in the dataset. Therefore, the upscaled effect should be more uncertain in tropical regions, but in the current version, SE is not particularly higher in the tropics than in other areas, because the model just depends on climate, not ecosystem representativeness. I would like to add a further level of uncertainty to SE based on the number of studies per type of ecosystem in the dataset. The idea is that ecosystems that are poorly sampled should have a higher SE than ecosystems with plenty of measurements in the dataset. 
Ecosystem type does not appear to explain variation in the effect, but this could just be the consequence of the low sample size of some ecosystem types. Regardless, ecosystem type should be forced to be included in the model somehow so that poorly sampled ecosystems show larger uncertainty. How can I do this? Should I include ecosystem type as a random effect?
 A: If you add ecosystem as a main/interaction random effect then the model will indeed incorporate additional parameters for the variability of the dependent variable as function of ecosystems. This allows you to incorporate additional modes of uncertainty for predictions. This will be based on: 


*

*different predicted $\hat{y}$ Since you assume that ecosystems can have an effect, the model will also, likely, make the effects different. (But as you say this does not appear to be so much the case and your fit will give little significant effects. This is no problem) If you make some manual custom made model then you can have the effect size the same and just the error of the estimates vary as function of the ecosystem (heteroskedasticity, which is easy if you apply fixed weights for different ecosystems and a bit more difficult if the weights are unknown)

*and also different standard error for $\hat{y}$, because the means, or $E(\hat{y})$, are being estimated with variations of sample size.

A: To address this comment explicitly as an extension of your original question: 

Would you add 'Biome' as a moderator or as a random effect? i.e.: mod1
  <- rma.mv(es, var, data=df, random= ~1|Biome, mods= ~ 1 + MAP + MAT)
  or mod2 <- rma.mv(es, var, data=df, mods= ~ 1 + MAP + MAT + Biome)?

My instinct would be to incorporate biome as a moderator (i.e., as a predictor of your random effects). Using the predict function and feeding your model the relevant data (which would now include three variables), and specifying the levels= argument, you should get a series of estimates and intervals (specified to your level of uncertainty). For biomes that contain more variability in their effect sizes (either as a function of small samples or as a result of greater variation overall) you should see wider intervals. Those intervals can be incorporated into simple tables or even included in graphics, depending on your final goals.
