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

interpretationmeta-analysismetaformultilevel-analysisr

Dears,

I am having trouble making sense of two MLMEs I did using rma.mv(). I am using a factor variable with 4 levels as a moderator. The two outputs (once with intercept, and once without it) are spitting out different p-values, and different estimates. The estimates are all Pearson's Rs.

This is the output for the model with intercept (model 1):

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 54) = 1.7411, p-val = 0.1695

level estimate pval
intrcpt 0.1232 0.0150*
LevelB 0.0679 0.2425
LevelC -0.0242 0.6699
LevelD 0.0227 0.0227

This is the output for the model without intercept (model 2):

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 54) = 22.3876, p-val < .0001

level estimate pval
LevelA 0.1232 0.0150*
LevelB 0.1911 <.0001***
LevelC 0.0990 0.0008***
LevelD 0.1459 <.0001***

Q1) I want to know which level is significant in its impact on the outcome? Where is this information present? Model 1 or Model 2?

Q2) I want to know how much is the effect size for each level on the outcome? (i.e. which estimate I should look at?). Again, where is this information, model 1 or 2?

Best Answer

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.