Solved – Metafor rma.mv function: missing estimates for two levels of a categorical moderator

meta-analysismultilevel-analysisr

I am conducting a meta-analysis with several categorical moderators and one continuous moderator; I am also interested in an interaction between two of the moderators (outcome and scale).

I have used the following code to do this:

Meta1<-rma.mv(yi=G, V=VG, mods = ~ Outcome * Scale + Scale_Type + Profac_def + Gender + Bias, random = list(~ 1 | Author_Num, ~ 1 | Study.in.Author, ~ 1 | Scale.in.Study, ~ 1 | Outcome.in.Study))

Followed by removing the intercept to examine estimates for the level of each factor (and the scale*outcome combinations) with the following:

Meta1.1<-rma.mv(yi=G, V=VG, mods = ~ Outcome * Scale + Scale_Type + Profac_def + Gender + Bias - 1, random = list(~ 1 | Author_Num, ~ 1 | Study.in.Author, ~ 1 | Scale.in.Study, ~ 1 | Outcome.in.Study))

Both outcome and scale have a lot of levels (12 and 19, respectively). However, the output is not showing all levels of the scale factor, or all of the outcome*scale combinations.

In the model with the intercept, one level of the scale moderator is missing, and a large number of combinations. In the model without the intercept, two levels of scale are missing (including the reference category), and again a large number of combinations.

I am very new to R and have no idea what is causing this. If anyone has any solutions and can advise how I can get around this it would be greatly appreciated

Best Answer

If you examine the output carefully, you should have gotten the following warning when fitting your models:

Warning message:
In rma.mv(...) :
  Redundant predictors dropped from the model.

The problem is that the model matrix (which is formed based on the moderators that you include in the model) is not of full rank. Put another way: There is not sufficient data available to estimate each of the coefficients in the model. As a result, the function has dropped moderators/predictors from the model that are redundant (so that the reduced model matrix is of full rank).

You note that Outcome and Scale have 12 and 19 levels, respectively. That alone means that a total of $12\times19 = 228$ coefficients would be need to be estimated, plus additional coefficients for the other moderators. Unless you have thousands of data points (and each combination of Outcome and Scale actually occurs in your dataset), this isn't going to work.

You need to consider simplifying your model, for example by collapsing some of the levels of the Outcome and Scale variables into a smaller number of levels. Even then, you need to check that each combination of levels actually occurs. You can easily check this by examining a contingency table of those factors:

table(Outcome, Scale)

You may also need to reconsider whether it is even possible/realistic to estimate the interaction between those two factors. Maybe you need to just stick to a model with main effects only.

Related Solutions

Solved – Metafor package: Interpreting meta-regression model

Yes, based on what you have shown, I would say that the analysis is sensible. One concern might be the relatively large number of moderator variables (or more specifically, model coefficients) relative to the number of estimates. Right now, you have $105 / 14 = 7.5$ estimates per coefficient (not counting the intercept). Some might want that ratio to be closer to 10 or even 15, but some might also be okay with a ratio of 5. None of these are right or wrong, but the lower the ratio, the more concerned I would be with overfitting.
Indeed, strictly speaking, the PSS version factor fails to be significant at $\alpha = .05$. However, I think you can still discuss this factor -- cautiously. Based on psychometric theory and all else equal, it is to be expected that longer versions would lead to higher reliability, which is indeed what you find here (although the 14-item version does not seem to yield, on average, higher reliability than the 10-item version -- maybe those 4 extra items are not as internally consistent as the rest or maybe there is something else that is different about studies examining the 14-item version that is not captured by all the other moderator variables already included in the model).
It is common practice to examine one moderator at a time. In principle, this is poor practice, since moderator variables are often correlated. So, fitting a model including multiple moderators (as you have done) would be better, as that gets you closer to examining the contribution of a particular moderator variable while controlling for the rest. One reason why this is often not done is that the dataset typically looks like Swiss cheese, with lots of holes (i.e., missing data) in it. After listwise deletion, one then ends up with a (much) smaller dataset (i.e., only the studies with complete information on all moderator variables). Besides the loss of information itself, when this happens, a major concern here is potential bias due to the missingness. Hence, instead, analyses are often conducted one moderator at a time, so that all of the studies providing information on a particular moderator variable can be used. Bias due to missingness may still be an issue here, but maybe less so. But, as mentioned at the beginning, you are then not controlling for other moderator variables, so a "fake" moderator might appear to be relevant simply because it is correlated with a "true" moderator.

There are fancy techniques to deal with missingness (e.g., multiple imputation, full information maximum likelihood estimation), but these methods are poorly developed in the meta-analytic context. Alternatively, you could run the 'full model' analysis and the 'one at a time' analyses and put them side-by-side and hopefully you find some consistency in the conclusions. If so, the discussion section will be easy to write. If not, then good luck ;)

Meta Analysis in R – Interpreting a Multi-Level Mixed Effects Model with Intercept and Without Intercept

Welcome to CV @Am95!

@Wolfgang's linked page does an excellent job explaining what's happening here in detail, but it sounds like perhaps you're struggling to see the proverbial forest for the trees? I'll borrow some of @Wolfgang's notation and try to put his page another way (more in the context of your provided example) in a way that is hopefully helpful.

The Big Picture/Short-Story

The two outputs (once with intercept, and once without it) are spitting out different p-values, and different estimates.

Indeed, but that is because you have specified two different models, that are set up (or "parameterized") to provide answers to two different questions. Therefore, their overarching F-test, included/excluded parameters, estimates, and p-values are (mostly) different as well.

The very quick answers to your questions are that both the significance of (Q1) and the effect size for (Q2) each level are provided by different columns of your Model 2 (let's call it the intercept-removed model). However, most of what you need to report (in terms of what is both useful and normative of meta-analytic moderator testing) is available in Model 1 (let's call it the intercept-present model). Indeed, from my perspective--and it's just my opinion--the intercept-removed model is best thought of as a mere "programming hack" to quickly/painlessly get your estimates broken down by each category level, if you have evidence of differences based on the intercept-present model (or are nonetheless compelled to report separate estimates).

As @Jeremy Miles was attempting to explain in your related previous post, a lot of what's happening here has much more to do with how categorical predictors in linear models work (e.g., simple regression), rather than much of anything special to do with meta-analysis per se. If all you want is an answer to Q1 and Q2, then all you need to know is that the meta-analytic average correlation of each level is in the "estimate" column of the intercept-removed model, and the significance level of the test for each (against an $H_0$ value of 0) is in the adjacent "pval" column of the intercept-removed model.

If you'd like to understand what is mathematically and functionally different between the intercept-present and intercept-removed models, and what utility the intercept-present model provides, read on.

The Detailed Version/Longer-Story

Functionally, the intercept-removed model is allowing you to estimate a unique meta-analytic average correlation for each of LevelA - LevelD. However, it is often of (greater) interest to know whether the meta-analytic averages of these levels differ from one another (e.g., LevelA's correlation of 0.1232 is superficially different [i.e., to the naked eye they are not the same number] from LevelB's correlation of 0.1911, but are they statistically different?). This is what the intercept-present model allows you to determine, and having some statistical evidence of differences between levels is usually a normative precursor to presenting different estimates by level (as in the second table from the intercept-removed model), because if there's not any evidence of differences, then one correlation (e.g., estimating across LevelA-LevelD) will do just fine--no need to complicate things with four correlations.

Specifically, the intercept-present model is estimating the meta-analytic average correlation of LevelA ($\mu_A = 0.1232$ [notice, this estimate and its p-value are identical between models/tables]), while the estimates of LevelB - LevelD are now parameterized to capture the difference in meta-analytic average correlation between LevelA and meta-analytic average correlation of the level on that particular row of the table. For example, if you add the estimated difference between LevelB and LevelA ($\beta_B = 0.0679$) to the intercept value/LevelA's average ($\mu_A = 0.1232$), then you arrive back at LevelB's average (as in the second row of the second table; $\mu_B = 0.1911$).

And so, while the estimate of LevelB in the first table captures the difference between its average and the average of LevelA ($\mu_B - \mu_A$), the estimate of LevelB the second table is the average of LevelB ($\mu_B$). Likewise, the p-value in the first table for LevelB is for the test of the difference between correlations for LevelB and LevelA ($H_0 = \beta_B - \beta_A == 0$), whereas the p-value in the second table for LevelB is for the test of its average correlation against a null value of zero ($H_0 = \mu_B == 0$; LevelA has nothing to do with this test*). And finally, the F-tests for each model also tell you something different: while the F-test for the intercept-removed model tests whether any of the meta-analytic average correlations are different from zero the F-test for the intercept-present model tests whether any of the remaining meta-analytic average correlations are different from LevelA's average correlation. In this way, the F-tests and tabular output should (hopefully) make more sense now: you have four significant meta-analytic correlations for LevelA-LevelD that are all significantly different from 0 (individually and hence the sig. F-test, both for intercept-removed model), but they are all relatively comparable in magnitude so the estimated differences compared to LevelA are all quite small and non-significant (both the intercept-present model).

How Did this Happen?

You might wonder why the intercept-present model places such a heavy emphasis on LevelA (i.e., why are all the comparisons for other levels made against LevelA and not, say, LevelD?). The reason is because either you deliberately specified the model in this fashion, OR (my guess), you let the software make the decision for you, in terms of how to "code" the categorical moderator variable consisting of your four factors/levels. I'll avoid going too into the weeds on how these coding schemes work (Cohen, Cohen, West, and Aiken's (2002) chapter on categorical coding in linear models is a must-read for details), the the gist of the situation is that for any linear model including a categorical predictor, you need $G-1$ code-variables to fully and non-redundantly explain (in this case) the differences in meta-analytic averages between your four levels (LevelA - Levels). There are different coding schemes available, but your analysis has implemented "dummy coding", whereby one level is chosen as the referent group (i.e., as its mean estimated as the intercept em), and each other remaining level has a code-variable created for which the slope captures the difference between that group and the referent group.

In your case, the coding scheme for your effects looks like this... :

Group___dummy_B___dummy_C___dummy_D

LevelA____0__________0__________0

LevelB____1__________0__________0

LevelC____0__________1__________0

LevelD____0__________0__________1

The moderator model looks something like this:

$Y_i = \beta_0 + \beta_1*dummy_B + \beta_2*dummy_C + \beta_3*dummy_D$

What's captured in your estimated terms is therefore:

$\beta_0 == \mu_A$

$\beta_1 == \mu_B - \mu_A$

$\beta_2 == \mu_C - \mu_A$

$\beta_3 == \mu_D - \mu_A$

Which spelt-out means (kinda messilyy, but bear with me...):

Meta-analytic correlation of a given group = LevelA's average + difference-between-A-and-B * (if you want Group B's average) + difference-between-A-and-C * (if you want Group C's average) + difference-between-A-and-D * (if you want Group D's average).

Assuming you did not deliberately choose this approach, R has structured the moderator analysis this way because it has detected four levels of your categorical moderator variable, and it defaults to dummy coding categorical predictors in linear models. It has further defaulted to coding LevelA as your referent group because it needs to choose something, and so does this based on what level of your factor is first in alphabetical order (they all start with "Level", and so A would have chosen before B, C, or D). If you wanted a different way of breaking up the between-Level variance (e.g., comparing to the meta-analytic grand-average correlation instead of LevelA's) or to keep with dummy coding but choosing a different level for the referent category (e.g., comparing all against LevelB), you'd likely need to manually recode (e.g., using effects-coding) or refactor (with a different order of levels) your categorical variable.

Suggested Reading

Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences. Guilford Press.

Best Answer

Related Solutions

Solved – Metafor package: Interpreting meta-regression model

Meta Analysis in R – Interpreting a Multi-Level Mixed Effects Model with Intercept and Without Intercept

Related Question