Is my full model with all significant estimates really the most parsimonious model?

I am performing a model selection based on BIC (Bayesian information criterion). I start with the following model fitted with rma.mv() in the 'metafor' package:

res<-rma.mv(yi~Nadd*(MAT_e+MAP_e+Ndep_e+Fert+Duration.yr),vi,data=pfg,random=list(~ 1 | Experiment,~1|Species),method="ML")


I get the following output:

Multivariate Meta-Analysis Model (k = 773; method: ML)

Variance Components:

estim    sqrt  nlvls  fixed      factor
sigma^2.1  0.6877  0.8293     37     no  Experiment
sigma^2.2  4.3396  2.0832    364     no     Species

Test for Residual Heterogeneity:
QE(df = 763) = 12090.4603, p-val < .0001

Test of Moderators (coefficient(s) 2:10):
QM(df = 9) = 74.6533, p-val < .0001

Model Results:

estimate      se     zval    pval    ci.lb    ci.ub
intrcpt            -0.8997  0.5141  -1.7501  0.0801  -1.9074   0.1079    .
Nadd                0.0019  0.0006   3.0630  0.0022   0.0007   0.0032   **
MAT_e              -0.1019  0.0326  -3.1251  0.0018  -0.1659  -0.0380   **
Ndep_e              0.0270  0.0205   1.3173  0.1877  -0.0131   0.0671
FertbNH4            0.6857  0.3939   1.7407  0.0817  -0.0864   1.4579    .
Duration.yr         0.0810  0.0441   1.8363  0.0663  -0.0055   0.1674    *
Nadd:MAT_e          0.0005  0.0002   2.5786  0.0099   0.0001   0.0010   **
Nadd:Ndep_e        -0.0003  0.0001  -2.3002  0.0074  -0.0005  -0.0000   **
Nadd:FertbNH4      -0.0061  0.0015  -4.1442  <.0001  -0.0090  -0.0032  ***
Nadd:Duration.yr   -0.0004  0.0001  -3.9385  <.0001  -0.0006  -0.0002  ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


I can't go further in the selection as you always get higher BIC by dropping any others of the terms. It is very wierd that all interaction terms are significant to me. Could this be interpeted as overfitting due to the large number of levels at the random part of the model (1|Species)?

• Although they are statistically significant the actual estimates look rather small (not sure what scale they are on) so i suspect your diagnosis may be correct: you have lots of data. – mdewey Nov 22 '17 at 12:04

You are adding random effects/intercepts for each level of Experiment and Species, but you are not adding random effects/intercepts corresponding to each individual estimate. That is a common mistake in fitting models of this type (in essence, you are assuming that the underlying true effects within experiments and the same species are homogeneous, which may not be the case). See also:
pfg\$Id <- 1:nrow(pfg)