Are there times when a multi-level meta-analytic model is discouraged?

Question

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

Study    d    Variance    Category
1        0    0.1         A
1        5    0.1         B
1        10   0.1         C
2        20   0.1         A
2        25   0.1         B

That is, I'm missing data for Category C for Study 2.

I then fit two meta-analytic models on this data with Category as the moderator in both cases. Analysis1 is a two-level model model, with Study as the random factor,

Analysis1 <- rma.mv(d, Variance, random = ~ 1 | Study, mods = ~ factor(Category) - 1, data=TheData)

which gives me the following output:

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

Variance Components: 

              estim     sqrt  nlvls  fixed  factor
sigma^2    199.9504  14.1404      2     no   Study

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

Test of Moderators (coefficient(s) 1,2,3): 
QM(df = 3) = 626.6563, p-val < .0001

Model Results:

                   estimate       se    zval    pval    ci.lb    ci.ub   
factor(Category)A   10.0000  10.0013  0.9999  0.3174  -9.6021  29.6021   
factor(Category)B   15.0000  10.0013  1.4998  0.1337  -4.6021  34.6021   
factor(Category)C   19.9975  10.0050  1.9987  0.0456   0.3880  39.6070  *

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

Analysis2 is a single-level model (and therefore without any random factors),

Analysis2 <- rma.mv(d, Variance, mods = ~ factor(Category) - 1, data=TheData)

which gives me the following output:

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

Variance Components: none

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

Test of Moderators (coefficient(s) 1,2,3): 
QM(df = 3) = 7500.0000, p-val < .0001

Model Results:

                   estimate      se     zval    pval    ci.lb    ci.ub     
factor(Category)A   10.0000  0.2236  44.7214  <.0001   9.5617  10.4383  ***
factor(Category)B   15.0000  0.2236  67.0820  <.0001  14.5617  15.4383  ***
factor(Category)C   10.0000  0.3162  31.6228  <.0001   9.3802  10.6198  ***

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

As can be seen, when using the single-level model, Analysis2, each category estimate is simply the mean of the values for that category. However, when using the three-level model, Analysis1, the estimates for Category A and B are the same as in the single-level model, while Category C seems to be estimated based on an extrapolation of what its value would/should be in Study 2. That is, what seems to be going on (put in an extremely non-mathematical language, since I haven't fully grasped the inner workings of how these analyses are performed) is that the model goes: "Oh, there's no data for Category C in study 2. Well, Category C was larger than both A and B in Study 1 so I guess it should be larger than Category A and B in Study 2 as well."

This is probably fine for a lot of cases, but are there situations when it's simply not appropritate to let the model extrapolate in this way? For example, I have a dataset where the different categories represent different categories of tests that the participants in psychological studies participated in. That is, the categories are things like verbal, visuo-spatial, and tactile, where each category merely designates the type of the test and not a specific test (for example, the category tactile could be applicable on hundereds of different tests). Further, all categories are never represented within the same study, so the data is just filled with missing data.

Here, it seems quite strange to extrapolate when trying to estimate the actual values of each category, especially since I have so many different categories as well (is it realistic to think that the extrapolation would be any good when I, for example, only have 2 out of 9 possible categories represented for a given study?). It simply seems messy to apply this method on this type of data. Am I on to something here?

Answer 1

This is indeed a peculiar result. Let me try to explain what is going on here. If you look at category A, you see estimates of 0 and 20 and since the sampling variances are the same, we can just take their mean to get the fixed effect estimate, which is 10. So in study 1, the estimate is -10 below that and in study 2 it is +10 above that. Analogously, within category B, the estimates are 5 and 25 for a mean of 15 and again the estimate in study 1 is -10 below that and in study it is +10 above that. So roughly, the random effect is -10 for study 1 and +10 for study 2. This isn't 100% right, but close enough to understand what is happening here.

Now in category C, you only have that one estimate of 10, which comes from study 1. Since the random effect for study 1 is -10, the fixed effect estimate for category C must be roughly 20, so that we end up with the estimate of 10 we actually see. That is why the model estimate is getting pushed up. So, this is not due to some kind of extrapolation of what the estimate would/should have been in study 2, but a simple consequence of having added a random effect at the study level (with such a peculiar dataset).

We can also understand the estimated variance for the study random effects this way. If the random effects are -10 and +10 (and so their mean is 0), then their variance is $((-10 - 0)^2 + (10 - 0)^2)/(2-1) = 200$, which is almost that value of $199.9504$ we get from the model. Note that I am dividing by $n-1$ here, since you are using REML estimation, and I am ignoring the sampling variances, so the model estimate isn't quite the same as that 200 -- but close enough. This will get you even closer:

rma.mv(d, .0001, random = ~ 1 | Study, mods = ~ factor(Category) - 1, data=TheData, control=list(optimizer="optim"))

(the syntax above is a shortcut for constraining all of the sampling variance to .0001; I had to switch optimizers though, because nlminb was having difficulties).

Note that in practice, one should not just have random effects for study, but also at the estimate level. Please take a look at:

http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011

where I discuss the three-level model at length. So one should use:

TheData$id <- 1:5
rma.mv(d, .0001, random = ~ 1 | Study/id, mods = ~ factor(Category) - 1, data=TheData)

This doesn't really get you around the issue observed above (the estimate level variance component is estimated to be 0, so nothing is changed), but it is a common mistake to forget this, so it is worth repeating.

Interestingly, if we just stick to random effects at the estimate level, then things make more intuitive sense:

rma.mv(d, .0001, random = ~ 1 | id, mods = ~ factor(Category) - 1, data=TheData)

So, I don't think the solution is to switch to a model without any random effects, but to understand what is going on and to possibly adjust the random effects structure.