tips:regression_with_rma

The standard linear regression model is given by $$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \ldots + \beta_p x_{ip} + e_i,$$ where $e_i \sim N(0, \sigma^2)$. Models of this sort can be fitted with the R function `lm()`

. The mixed-effects meta-regression model is given by $$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \ldots + \beta_p x_{ip} + u_i + e_i,$$ where $u_i \sim N(0, \tau^2)$ and $e_i \sim N(0, v_i)$, where $v_i$ are the (approximately) known sampling variances of the observed outcomes or effect size estimates (e.g., standardized mean differences, log odds ratios). Models of this sort can be fitted with the function `rma()`

from the metafor package.

Consequently, if one were to set $v_i = 0$ for all outcomes, then the standard linear regression model and the mixed-effects meta-regression model are actually identical (with $\sigma^2$ denoting the same parameter as $\tau^2$). This equivalence can be demonstrated with an arbitrary dataset:

library(metafor) stackloss

Air.Flow Water.Temp Acid.Conc. stack.loss 1 80 27 89 42 2 80 27 88 37 3 75 25 90 37 4 62 24 87 28 5 62 22 87 18 6 62 23 87 18 7 62 24 93 19 8 62 24 93 20 9 58 23 87 15 10 58 18 80 14 11 58 18 89 14 12 58 17 88 13 13 58 18 82 11 14 58 19 93 12 15 50 18 89 8 16 50 18 86 7 17 50 19 72 8 18 50 19 79 8 19 50 20 80 9 20 56 20 82 15 21 70 20 91 15

See `help(stackloss)`

for more details on this dataset. Most importantly, variable `stack.loss`

is the dependent variable with `Air.Flow`

, `Water.Temp`

, and `Acid.Conc`

serving as potentially relevant predictors.

Now let's fit the standard linear regression model to these data with both the `lm()`

and the `rma()`

functions with:

res.lm <- lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data=stackloss)

and

stackloss$vi <- 0 res.rma <- rma(stack.loss, vi, mods = ~ Air.Flow + Water.Temp + Acid.Conc., data=stackloss)

Note that all sampling variances are set to 0 for `rma()`

(the function will actually issue a warning that the dataset includes outcomes with non-positive sampling variances – which would be rather strange in the meta-analytic context – but this can be safely ignored here).

We can now compare the output from the two models with:

summary(res.lm)

Call: lm(formula = stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data = stackloss) Residuals: Min 1Q Median 3Q Max -7.2377 -1.7117 -0.4551 2.3614 5.6978 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -39.9197 11.8960 -3.356 0.00375 ** Air.Flow 0.7156 0.1349 5.307 5.8e-05 *** Water.Temp 1.2953 0.3680 3.520 0.00263 ** Acid.Conc. -0.1521 0.1563 -0.973 0.34405 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.243 on 17 degrees of freedom Multiple R-squared: 0.9136, Adjusted R-squared: 0.8983 F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09

and

res.rma

Mixed-Effects Model (k = 21; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 10.5194 (SE = 3.6081) tau (square root of estimated tau^2 value): 3.2434 R^2 (amount of heterogeneity accounted for): 89.83% Test of Moderators (coefficient(s) 2:4): QM(df = 3) = 179.7067, p-val < .0001 Model Results: estimate se zval pval ci.lb ci.ub intrcpt -39.9197 11.8960 -3.3557 0.0008 -63.2354 -16.6039 *** Air.Flow 0.7156 0.1349 5.3066 <.0001 0.4513 0.9800 *** Water.Temp 1.2953 0.3680 3.5196 0.0004 0.5740 2.0166 *** Acid.Conc. -0.1521 0.1563 -0.9733 0.3304 -0.4585 0.1542 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

A few things are worth noting:

- The estimated model coefficients, corresponding standard errors, and the test statistics are exactly the same. However,
`lm()`

computes the p-values based on the t-distribution, while`rma()`

uses (by default) the standard normal distribution.

- The omnibus test of the model coefficients conducted by
`lm()`

is an F-test (`F-statistic: 59.9 on 3 and 17 DF, p-value: 3.016e-09`

), while`rma()`

uses by default a chi-square test (`QM(df = 3) = 179.7067, p-val < .0001`

).

- The estimated
`Residual standard error: 3.243`

is the same as the estimate of`tau: 3.2434`

.

- The (pseudo) R-squared value reported by
`rma()`

(`89.83%`

) actually corresponds to the adjusted R-squared value provided by`lm()`

(`0.8983`

).

- The $I^2$ and $H^2$ statistics typically reported by the
`rma()`

function are missing, since these statistics cannot be computed when the dataset includes outcomes with non-positive sampling variances. Similarly, the usually reported results from the $Q$-test for heterogeneity are omitted for the same reason.

To get full correspondence between the two models, we can use the 'Knapp & Hartung' method when fitting the model with the `rma()`

function:

res.rma <- rma(stack.loss, vi, mods = ~ Air.Flow + Water.Temp + Acid.Conc., data=stackloss, test="knha") res.rma

Mixed-Effects Model (k = 21; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 10.5194 (SE = 3.6081) tau (square root of estimated tau^2 value): 3.2434 R^2 (amount of heterogeneity accounted for): 89.83% Test of Moderators (coefficient(s) 2:4): F(df1 = 3, df2 = 17) = 59.9022, p-val < .0001 Model Results: estimate se tval pval ci.lb ci.ub intrcpt -39.9197 11.8960 -3.3557 0.0038 -65.0180 -14.8213 ** Air.Flow 0.7156 0.1349 5.3066 <.0001 0.4311 1.0002 *** Water.Temp 1.2953 0.3680 3.5196 0.0026 0.5188 2.0717 ** Acid.Conc. -0.1521 0.1563 -0.9733 0.3440 -0.4819 0.1776 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Now the p-values and the results from the omnibus test of the moderators also match up completely with the output from the `lm()`

function.

tips/regression_with_rma.txt · Last modified: 2018/11/23 08:58 (external edit)