tips:models_with_or_without_intercept

**This is an old revision of the document!**

A question that comes up frequently is the proper interpretation of meta-regression models with or without an intercept term. Below I illustrate the difference using the dataset for the BCG vaccine meta-analysis (Colditz et al., 1994).

The dataset is called `dat.bcg`

. The variables `tpos and `

tneg` in the dataset indicate the number of tuberculosis cases and non-cases in the treated (vaccinated) groups, respectively. Similarly, variables `

cpos and `cneg`

indicate the number of cases and non-cases in the control (non-vaccinated) groups. Based on these four variables, the (log) risk ratios and corresponding sampling variances for the 13 studies can be obtained with:

library(metafor) dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg) dat

trial author year tpos tneg cpos cneg ablat alloc yi vi 1 1 Aronson 1948 4 119 11 128 44 random -0.8893 0.3256 2 2 Ferguson & Simes 1949 6 300 29 274 55 random -1.5854 0.1946 3 3 Rosenthal et al 1960 3 228 11 209 42 random -1.3481 0.4154 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random -1.4416 0.0200 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate -0.2175 0.0512 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate -0.7861 0.0069 7 7 Vandiviere et al 1973 8 2537 10 619 19 random -1.6209 0.2230 8 8 TPT Madras 1980 505 87886 499 87892 13 random 0.0120 0.0040 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random -0.4694 0.0564 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic -1.3713 0.0730 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic -0.3394 0.0124 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic 0.4459 0.5325 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic -0.0173 0.0714

Values of the log risk ratio (variable `yi`

) below 0 indicate that the infection risk was lower in the vaccinated group compared to the control group.

As we can see above, 7 studies used random assignment to determine which study participants would receive the BCG vaccine, 2 studies used alternating assignment (e.g., first person is vaccinated, second person isn't, third person is vaccinated, next person isn't, and so on), and finally 4 studies used systematic assignment (e.g., where assignment to the vaccination group occurred on the basis of the birth year, as in Comstock et al., 1974).

We will now fit a (mixed-effects) meta-regression model to these data, including the method of treatment allocation (`alloc`

) as a predictor/moderator. Note that `alloc`

is a character variable, which will be treated as a factor in the model. The standard way of handling a factor in regression models is to create dummy variables representing the various levels. One of these "dummies" is then left out of the model, with the corresponding level becoming the "reference" level. Fortunately, we do not have to do this coding ourselves, as R will handle this for us.

res <- rma(yi, vi, mods = ~ factor(alloc), data=dat) res

Mixed-Effects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111) tau (square root of estimated tau^2 value): 0.6013 I^2 (residual heterogeneity / unaccounted variability): 88.77% H^2 (unaccounted variability / sampling variability): 8.91 R^2 (amount of heterogeneity accounted for): 0.00% Test for Residual Heterogeneity: QE(df = 10) = 132.3676, p-val < .0001 Test of Moderators (coefficients 2:3): QM(df = 2) = 1.7675, p-val = 0.4132 Model Results: estimate se zval pval ci.lb ci.ub intrcpt -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468 factor(alloc)random -0.4478 0.5158 -0.8682 0.3853 -1.4588 0.5632 factor(alloc)systematic 0.0890 0.5600 0.1590 0.8737 -1.0086 1.1867 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

By default, `alternate`

was made the reference level (since the first letter "a" comes before "r" and "s" in the alphabet), so the corresponding dummy variable for this level has been left out.^{1)} Therefore, the model intercept (coefficients $b_0$) is the estimated (average) log risk ratio for alternating allocation. The coefficients for `factor(alloc)random`

(coefficient $b_1$) and `factor(alloc)systematic`

(coefficient $b_2$) indicate how much lower/higher the estimated (average) log risk ratio is for random and systematic allocation, respectively, compared to alternating allocation. So these coefficients reflect contrasts between these two levels compared to the reference level. Neither coefficient is significantly different from 0 (as indicated by the corresponding p-values), but this isn't a proper test of the factor as a whole.

The results given under "Test of Moderators" provides a joint or omnibus test of all model coefficients except for the intercept (`coefficients 2:3`

). Careful: The numbering of the coefficients here starts at 1, corresponding to the intercept (`intrcpt`

), with coefficients 2 and 3 corresponding to `factor(alloc)random`

and `factor(alloc)systematic`

. The results indicate that we cannot reject the null hypothesis $H_0: \beta_1 = \beta_2 = 0$. Since we cannot reject the null hypothesis that random and systematic allocation are different from alternating allocation, then this implies that the three-level factor as a whole is not significant.

We can also test one or multiple model coefficients ourselves (and contrasts between them, as we will see further below) using the `anova()`

function. The same omnibus test as above can be obtained with:

anova(res, btt=2:3)

Test of Moderators (coefficients 2:3): QM(df = 2) = 1.7675, p-val = 0.4132

The `btt`

argument stands for "betas to test" and is used to specify which coefficients we want include in the test.

A different way to conduct the same test is to use the `L`

argument, which allows us to specify one or more vectors of numbers, which are multiplied with the model coefficients. In particular, we can use:

anova(res, L=rbind(c(0,1,0),c(0,0,1)))

to test the two hypotheses \begin{align} &H_0: (0 \times \beta_0) + (1 \times \beta_1) + (0 \times \beta_2) = 0 \; \Leftrightarrow \; \beta_1 = 0 \\ &H_0: (0 \times \beta_0) + (0 \times \beta_1) + (1 \times \beta_2) = 0 \; \Leftrightarrow \; \beta_2 = 0 \end{align} and the output then again includes the omnibus test of both of these two hypotheses:

Hypotheses: 1: factor(alloc)random = 0 2: factor(alloc)systematic = 0 Results: estimate se zval pval 1: -0.4478 0.5158 -0.8682 0.3853 2: 0.0890 0.5600 0.1590 0.8737 Omnibus Test of Hypotheses: QM(df = 2) = 1.7675, p-val = 0.4132

As discussed above, the coefficients for `random`

($\beta_1$) and `systematic`

($\beta_2$) reflect contrasts between these two levels and alternating allocation, that is, if we denote $\mu_a$, $\mu_r$, and $\mu_s$ as the (average) log risk ratio for alternating, random, and systematic allocation, respectively, then the model coefficients represent
\begin{align}
&\beta_0 = \mu_a \\
&\beta_1 = \mu_r - \mu_a \\
&\beta_2 = \mu_s - \mu_a.
\end{align}
But what about the contrast between random and systematic allocation? It turns out that we can obtain this from the model as the difference between the $\beta_1$ and $\beta_2$ coefficients. In particular, if we subtract $\beta_1$ from $\beta_2$, then
$$
\beta_2 - \beta_1 = (\mu_r - \mu_a) - (\mu_s - \mu_a) = \mu_r - \mu_s
$$
so this difference reflects how different random allocation is compared to systematic allocation. Using the `anova()`

function, we can obtain this contrast with

anova(res, L=c(0,-1,1))

Hypothesis: 1: -factor(alloc)random + factor(alloc)systematic = 0 Results: estimate se zval pval 1: 0.5369 0.4364 1.2303 0.2186 Test of Hypothesis: QM(df = 1) = 1.5137, p-val = 0.2186

The results of the omnibus test of a factor will be identical regardless of which level is chosen as the reference level. To illustrate this, let's make `random`

allocation the reference level, which we can do using the `relevel()`

function:

res <- rma(yi, vi, mods = ~ relevel(factor(alloc), ref="random"), data=dat) res

Mixed-Effects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111) tau (square root of estimated tau^2 value): 0.6013 I^2 (residual heterogeneity / unaccounted variability): 88.77% H^2 (unaccounted variability / sampling variability): 8.91 R^2 (amount of heterogeneity accounted for): 0.00% Test for Residual Heterogeneity: QE(df = 10) = 132.3676, p-val < .0001 Test of Moderators (coefficients 2:3): QM(df = 2) = 1.7675, p-val = 0.4132 Model Results: estimate se zval pval ci.lb ci.ub intrcpt -0.9658 0.2672 -3.6138 0.0003 -1.4896 -0.4420 *** alternate 0.4478 0.5158 0.8682 0.3853 -0.5632 1.4588 systematic 0.5369 0.4364 1.2303 0.2186 -0.3184 1.3921 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(I shortened the names of the coefficients in the output above to make the table under the `Model Results`

more readable). Now the intercept reflects the estimated (average) log risk ratio for random allocation, while the coefficients for `alternate`

and `systematic`

are again contrasts of these two levels compared to random allocation. Note how the coefficient for systematic allocation is the same as we obtained earlier using `anova(res, L=c(0,-1,1))`

. Moreover, as we can see in the output, the results for the omnibus test of these two coefficients is identical to what we obtained earlier.

Now we will fit the model, but remove the intercept, which can be done with:

res <- rma(yi, vi, mods = ~ factor(alloc) - 1, data=dat) res

Mixed-Effects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111) tau (square root of estimated tau^2 value): 0.6013 I^2 (residual heterogeneity / unaccounted variability): 88.77% H^2 (unaccounted variability / sampling variability): 8.91 Test for Residual Heterogeneity: QE(df = 10) = 132.3676, p-val < .0001 Test of Moderators (coefficients 1:3): QM(df = 3) = 15.9842, p-val = 0.0011 Model Results: estimate se zval pval ci.lb ci.ub factor(alloc)alternate -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468 factor(alloc)random -0.9658 0.2672 -3.6138 0.0003 -1.4896 -0.4420 *** factor(alloc)systematic -0.4289 0.3449 -1.2434 0.2137 -1.1050 0.2472 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

When we remove the intercept from the model, then all dummy variables are included in the model as predictors. Hence, in this example, each coefficient then reflects the estimated (average) log risk ratio for a particular level (i.e., in the notation used earlier, $\hat{\mu}_a$, $\hat{\mu}_r$, and $\hat{\mu}_s$, respectively).^{2)}

The results given under "Test of Moderators" are now for an omnibus test of all three model coefficients. This is a joint test of the three null hypotheses \begin{align} &H_0: \mu_a = 0 \\ &H_0: \mu_r = 0 \\ &H_0: \mu_s = 0 \end{align} which we can also write more compactly as $$ H_0: \mu_a = \mu_r = \mu_s = 0. $$ This test is clearly significant ($p = .0011), which indicates that we can reject the null hypothesis that the (average) log risk ratio is zero for all three methods of allocation.

Again, we could use the `anova()`

function to carry out explicitly the same test with:

anova(res, L=rbind(c(1,0,0),c(0,1,0),c(0,0,1)))

Hypotheses: 1: factor(alloc)alternate = 0 2: factor(alloc)random = 0 3: factor(alloc)systematic = 0 Results: estimate se zval pval 1: -0.5180 0.4412 -1.1740 0.2404 2: -0.9658 0.2672 -3.6138 0.0003 3: -0.4289 0.3449 -1.2434 0.2137 Omnibus Test of Hypotheses: QM(df = 3) = 15.9842, p-val = 0.0011

It is important to realize that this does not test whether there are differences between the different forms of allocation (this is what we tested earlier in the model that included the intercept term). However, we can again use contrasts of the model coefficients to test differences between the levels. For example, let's test the difference between alternating and random allocation and the difference between systematic allocation and random allocation:

anova(res, L=rbind(c(-1,1,0),c(-1,0,1)))

Hypotheses: 1: -factor(alloc)alternate + factor(alloc)random = 0 2: -factor(alloc)alternate + factor(alloc)systematic = 0 Results: estimate se zval pval 1: -0.4478 0.5158 -0.8682 0.3853 2: 0.0890 0.5600 0.1590 0.8737 Omnibus Test of Hypotheses: QM(df = 2) = 1.7675, p-val = 0.4132

These are now the exact same results we obtained earlier for the model that included the intercept term.

What the example above shows is that, whether we remove the intercept or not, we are essentially fitting the same model, but using a different parameterization. The Wikipedia page linked to here discusses the idea of parameterization in the context of geometry, but this is directly relevant, since fundamentally the process of fitting (meta)-regression models can also be conceptualized geometrically (involving projections and vector spaces). An excellent (but quite technical) reference for this perspective is Christensen (2011).

A rather different case arises if the model we are fitting only involves one or multiple continuous moderators. For example, consider the following model:

res <- rma(yi, vi, mods = ~ ablat, data=dat)

Mixed-Effects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.0764 (SE = 0.0591) tau (square root of estimated tau^2 value): 0.2763 I^2 (residual heterogeneity / unaccounted variability): 68.39% H^2 (unaccounted variability / sampling variability): 3.16 R^2 (amount of heterogeneity accounted for): 75.62% Test for Residual Heterogeneity: QE(df = 11) = 30.7331, p-val = 0.0012 Test of Moderators (coefficient 2): QM(df = 1) = 16.3571, p-val < .0001 Model Results: estimate se zval pval ci.lb ci.ub intrcpt 0.2515 0.2491 1.0095 0.3127 -0.2368 0.7397 ablat -0.0291 0.0072 -4.0444 <.0001 -0.0432 -0.0150 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here, the absolute latitude (`ablat`

) of the study site was included as a continuous moderator (ranging from 13 to 55 degrees). The idea was to examine if the higher prevalence of nonpathogenic environmental mycobacteria for sites nearer the equator may provide a natural immunity against tuberculosis, leading to lower treatments effects (i.e., log risk ratios close to zero). The results above support this hypothesis: The intercept, which reflects the estimated (average) log risk ratio at the equator (with a latitude of 0), is not significantly different from 0, while the coefficient/slope for `ablat`

suggests that the estimated (average) log risk ratio becomes increasingly more negative for sites further away from the equator (recall that negative log risk ratios here reflect that the vaccine is reducing the risk of an infection and hence is working).^{3)}

Now what happens if we remove the intercept?

res <- rma(yi, vi, mods = ~ ablat - 1, data=dat)

Mixed-Effects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.0865 (SE = 0.0592) tau (square root of estimated tau^2 value): 0.2940 I^2 (residual heterogeneity / unaccounted variability): 78.35% H^2 (unaccounted variability / sampling variability): 4.62 Test for Residual Heterogeneity: QE(df = 12) = 48.7021, p-val < .0001 Model Results: estimate se zval pval ci.lb ci.ub ablat -0.0226 0.0032 -7.0545 <.0001 -0.0289 -0.0163 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

When the model only includes continuous (i.e., numeric) predictors/moderators, then removing the intercept does just that: it removes the intercept. Hence, the model above forces the intercept to be 0, that is, it assumes that at the equator, the (average) log risk ratio is exactly zero. This seems like a rather strong assumption to make, so I would not recommend doing so. In fact, only in rare cases is it ever appropriate to remove the intercept term from a regression model.

Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., et al. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. *Journal of the American Medical Association, 271*(9), 698–702.

Christensen, R. (2011). *Plane answers to complex questions: The theory of linear models.* New York: Springer.

Since

`alloc`

is a character variable, the coercion to a factor via the `factor()`

function isn't strictly necessary (since this would have happened automatically anyway), but it's better to be explicit. In some cases, a numeric variable may in fact be used to represent levels of a categorical variable, in which case use of the `factor()`

function is crucial.Note how the coefficients for

`alternate`

and `random`

are exactly identical to the intercepts in the two models fitted earlier where each of these levels was chosen as the reference level.We should interpret the intercept only very cautiously, since it is an extrapolation beyond the range of values observed for absolute latitude (i.e., no study was conducted exactly at the equator).

tips/models_with_or_without_intercept.1561895020.txt.gz · Last modified: 2019/06/30 11:43 by 127.0.0.1

Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International