More moderators in a meta-analysis equals greater chance of reaching statistical significance?

Question

I have the following example data, stored in the variable TheData,

Study    ID    Category    Cohens_d    Variance   
1        1     A           0           0.1    
1        2     B           5           0.1
1        3     C           10          0.1
2        4     A           20          0.1
2        5     B           25          0.1

which I want to run a meta-analytic model with moderators on.

If I'm using the rma.mvfunction from the Metafor package, and run

rma.mv(Cohens_d, Variance, random = ~ 1 | ID, mods = ~ factor(Category) - 1, data=TheData)

I get the following result

Multivariate Meta-Analysis Model (k = 5; method: REML)

Variance Components: 

              estim     sqrt  nlvls  fixed  factor
sigma^2    199.9003  14.1386      5     no      ID

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

Test of Moderators (coefficient(s) 1,2,3): 
QM(df = 3) = 3.7500, p-val = 0.2898

Model Results:

                   estimate       se    zval    pval     ci.lb    ci.ub   
factor(Category)A   10.0000  10.0000  1.0000  0.3173   -9.5997  29.5997   
factor(Category)B   15.0000  10.0000  1.5000  0.1336   -4.5997  34.5997   
factor(Category)C   10.0000  14.1421  0.7071  0.4795  -17.7181  37.7181   

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

which is what I want, namely that each category is estimated to the average of the observed values (see this question for a discussion of using different random effects structures). So far, so good. However, if I add data for a fourth category, D, like so

Study    ID    Category    Cohens_d    Variance   
1        1     A           0           0.1    
1        2     B           5           0.1
1        3     C           10          0.1
1        4     D           0           0.1
2        5     A           20          0.1
2        6     B           25          0.1
2        7     D           10          0.1

this gives me the following result

Multivariate Meta-Analysis Model (k = 7; method: REML)

Variance Components: 

              estim     sqrt  nlvls  fixed  factor
sigma^2    137.4002  11.7218      7     no      ID

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

Test of Moderators (coefficient(s) 1,2,3,4): 
QM(df = 4) = 5.5454, p-val = 0.2358

Model Results:

                   estimate       se    zval    pval     ci.lb    ci.ub   
factor(Category)A   10.0000   8.2916  1.2060  0.2278   -6.2512  26.2512   
factor(Category)B   15.0000   8.2916  1.8091  0.0704   -1.2512  31.2512  .
factor(Category)C   10.0000  11.7260  0.8528  0.3938  -12.9826  32.9826   
factor(Category)D    2.5000   8.2916  0.3015  0.7630  -13.7512  18.7512   

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

Note that here, the standard error (and subsequently the p-values) of my original categories have changed, or more specifically shrinked. Also note the p-value for the overall test of moderators which shrinks as well when adding another category (in this case from 0.29 to 0.24).

Is this how the model is supposed to work, or have I somehow distorted it with my strange dataset? From this result, it seems like I could, if I'm not happy with my p-values (let's say I would like cateogory B to be significantly different from zero), just chuck in a couple of extra categories and call it a day (or make sure that I include a couple extra moderators from the start if I don't want to modify my analyses after the fact).

I'm not sure if I'm bumping up against some grander statistical problem/principle here, something that's not specifically related to only meta-analyses, but nevertheless, I'm still confused.

Answer 1

The SEs decrease because the estimate of (residual) heterogeneity has decreased, so that makes perfect sense. Try this dataset:

TheData <- structure(list(Study = c(1L, 1L, 1L, 1L, 2L, 2L, 2L), ID = 1:7, 
Category = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 4L), .Label = c("A", 
"B", "C", "D"), class = "factor"), Cohens_d = c(0, 5, 10, 
-5, 20, 25, 15), Variance = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 
0.1)), .Names = c("Study", "ID", "Category", "Cohens_d",  "Variance"), row.names = c(NA, -7L), class = "data.frame")

Then you get:

Multivariate Meta-Analysis Model (k = 7; method: REML)

Variance Components: 

              estim     sqrt  nlvls  fixed  factor
sigma^2    199.9003  14.1386      7     no      ID

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

Test of Moderators (coefficient(s) 1,2,3,4): 
QM(df = 4) = 4.0000, p-val = 0.4060

Model Results:

                   estimate       se    zval    pval     ci.lb    ci.ub   
factor(Category)A   10.0000  10.0000  1.0000  0.3173   -9.5997  29.5997   
factor(Category)B   15.0000  10.0000  1.5000  0.1336   -4.5997  34.5997   
factor(Category)C   10.0000  14.1421  0.7071  0.4795  -17.7181  37.7181   
factor(Category)D    5.0000  10.0000  0.5000  0.6171  -14.5997  24.5997   

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

Note that this includes category D, but now the estimate of residual heterogeneity is unchanged. And so are the SEs for categories A, B, and C.

Naturally, the omnibus test is changed, since it now includes an additional coefficient. But if you fit the model with btt=1:3 (so the omnibus test only includes the first three coefficients), then you get the same result.