The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools


tips:rma_vs_lm_lme_lmer

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
tips:rma_vs_lm_lme_lmer [2022/08/03 11:37] Wolfgang Viechtbauertips:rma_vs_lm_lme_lmer [2023/11/14 08:00] (current) Wolfgang Viechtbauer
Line 1: Line 1:
 ===== A Comparison of the rma() and the lm(), lme(), and lmer() Functions ===== ===== A Comparison of the rma() and the lm(), lme(), and lmer() Functions =====
  
-In essence, the commonly used meta-analytic models are just special cases of the linear (mixed-effects) model with the only peculiar aspect that the variances of the error terms (i.e., the sampling variances) are known. Not surprisingly, the question therefore comes up occasionally why one cannot use the ''lm()'', ''lme()'', and ''lmer()'' functions for conducting meta-analyses. After all, the functions can be used to fit linear (mixed-effects) models and both functions allow the user to specify the sampling variances via the ''weights'' argument. So it seems that one should also be able to fit meta-analytic models with these functions. Below, I describe and illustrate how the models fitted via the ''lm()'', ''lme()'', and ''lmer()'' functions differ from the models fitted by the ''rma()'' function and why the those functions are therefore not suitable for fitting meta-analytic models.+Commonly used meta-analytic models are just special cases of the general linear (mixed-effects) model with the only peculiar aspect that the variances of the error terms (i.e., the sampling variances) are known. Not surprisingly, the question therefore comes up occasionally why the ''lm()'', ''lme()'', and ''lmer()'' functions cannot be used to conduct a meta-analysis. After all, the functions can be used to fit various types of linear (mixed-effects) models and both functions allow the user to specify the sampling variances via the ''weights'' argument. So it seems that one should also be able to fit meta-analytic models with these functions. Below, I describe and illustrate how the models fitted via the ''lm()'', ''lme()'', and ''lmer()'' functions differ from the models fitted by the ''rma()'' function and why the those functions are therefore not suitable for fitting meta-analytic models.
  
 ==== Equal-Effects Model ==== ==== Equal-Effects Model ====
  
 Let's start with an equal-effects model. As an example, consider the data from the meta-analysis by Molloy et al. (2014) on the relationship between conscientiousness and medication adherence. For each study, we can compute the r-to-z transformed correlation coefficient and corresponding sampling variance with: Let's start with an equal-effects model. As an example, consider the data from the meta-analysis by Molloy et al. (2014) on the relationship between conscientiousness and medication adherence. For each study, we can compute the r-to-z transformed correlation coefficient and corresponding sampling variance with:
 +
 <code rsplus> <code rsplus>
 library(metafor) library(metafor)
Line 30: Line 31:
 16 Wiebe & Christensen 1997  65  0.040  0.0400 0.0161 16 Wiebe & Christensen 1997  65  0.040  0.0400 0.0161
 </code> </code>
 +
 The ''yi'' values are the r-to-z transformed correlations and the ''vi'' values the corresponding sampling variances. The ''yi'' values are the r-to-z transformed correlations and the ''vi'' values the corresponding sampling variances.
  
 We can now fit an equal-effects model to these data with: We can now fit an equal-effects model to these data with:
 +
 <code rsplus> <code rsplus>
 res.ee <- rma(yi, vi, data=dat, method="EE") res.ee <- rma(yi, vi, data=dat, method="EE")
Line 56: Line 59:
  
 Now let's try to fit the same model with the ''lm()'' function by specifying the inverse of the sampling variances as weights: Now let's try to fit the same model with the ''lm()'' function by specifying the inverse of the sampling variances as weights:
 +
 <code rsplus> <code rsplus>
 res.lm <- lm(yi ~ 1, weights = 1/vi, data=dat) res.lm <- lm(yi ~ 1, weights = 1/vi, data=dat)
Line 76: Line 80:
 Residual standard error: 1.595 on 15 degrees of freedom Residual standard error: 1.595 on 15 degrees of freedom
 </code> </code>
 +
 Two things are of note here: Two things are of note here:
   - The estimated intercept (''0.12518'') is exactly the same as the model coefficient obtained earlier with the ''rma()'' function (the value reported by the ''rma()'' function is rounded to 4 decimals, but that can be changed with ''print(res.ee, digits=5)'' yielding the same value).   - The estimated intercept (''0.12518'') is exactly the same as the model coefficient obtained earlier with the ''rma()'' function (the value reported by the ''rma()'' function is rounded to 4 decimals, but that can be changed with ''print(res.ee, digits=5)'' yielding the same value).
-  - The standard error (of the estimated intercept) is different (and hence, the test statistic and p-value also differs).+  - The standard error (of the estimated intercept) is different (and hence, the test statistic and p-value also differ).
  
 The reason for this discrepancy is that the model fitted by the ''lm()'' function assumes that the weights are not exactly known, but only up to a proportionality constant (namely $\sigma^2_e$, that is, the error or residual variance), which is estimated and then factored into the weights. The reason for this discrepancy is that the model fitted by the ''lm()'' function assumes that the weights are not exactly known, but only up to a proportionality constant (namely $\sigma^2_e$, that is, the error or residual variance), which is estimated and then factored into the weights.
  
 This can be demonstrated by extracting the estimated error variance from the ''lm'' object, multiplying the sampling variances by that value, and re-fitting the model with the ''rma()'' function: This can be demonstrated by extracting the estimated error variance from the ''lm'' object, multiplying the sampling variances by that value, and re-fitting the model with the ''rma()'' function:
 +
 <code rsplus> <code rsplus>
 rma(yi, vi*sigma(res.lm)^2, data=dat, method="EE") rma(yi, vi*sigma(res.lm)^2, data=dat, method="EE")
Line 103: Line 109:
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 </code> </code>
 +
 Now these results are exactly the same as those obtained by the ''lm()'' function. Note that multiplying the sampling variances by some constant value does not affect the model estimate, but only the corresponding standard error. Now these results are exactly the same as those obtained by the ''lm()'' function. Note that multiplying the sampling variances by some constant value does not affect the model estimate, but only the corresponding standard error.
  
 Therefore, if we want to obtain the same standard error from ''lm()'' as from ''rma()'', we must adjust the standard error by $\hat{\sigma}_e$ (i.e., the square-root of the estimated error variance or what is called the "residual standard error" in the output from the ''lm()'' function): Therefore, if we want to obtain the same standard error from ''lm()'' as from ''rma()'', we must adjust the standard error by $\hat{\sigma}_e$ (i.e., the square-root of the estimated error variance or what is called the "residual standard error" in the output from the ''lm()'' function):
 +
 <code rsplus> <code rsplus>
 coef(summary(res.ee)) coef(summary(res.ee))
Line 119: Line 127:
 [1] 0.01699805 [1] 0.01699805
 </code> </code>
 +
 +This is now the exact same standard error of the pooled estimate from an equal-effects model.
  
 ==== Random-Effects Model ==== ==== Random-Effects Model ====
  
 A random-effects model can be fitted to the same data using the ''rma()'' function with: A random-effects model can be fitted to the same data using the ''rma()'' function with:
 +
 <code rsplus> <code rsplus>
 res.re <- rma(yi, vi, data=dat) res.re <- rma(yi, vi, data=dat)
Line 148: Line 159:
  
 Now let's try to fit the same model with the ''lme()'' function from the ''nlme'' package. Here, the ''weights'' argument is used to specify a variance function with fixed variances (note: the name of the ''weights'' argument is a bit of a misnomer for the ''lme()'' function, as one does not use it to specify the weights -- as in the ''lm()'' function -- but a variance function structure). In addition, we need to add the random effects for each study via the ''random'' argument. Now let's try to fit the same model with the ''lme()'' function from the ''nlme'' package. Here, the ''weights'' argument is used to specify a variance function with fixed variances (note: the name of the ''weights'' argument is a bit of a misnomer for the ''lme()'' function, as one does not use it to specify the weights -- as in the ''lm()'' function -- but a variance function structure). In addition, we need to add the random effects for each study via the ''random'' argument.
 +
 <code rsplus> <code rsplus>
 library(nlme) library(nlme)
Line 181: Line 193:
  
 Alternatively, we could use the ''lmer()'' function from the ''lme4'' package. The syntax would be: Alternatively, we could use the ''lmer()'' function from the ''lme4'' package. The syntax would be:
 +
 <code rsplus> <code rsplus>
 library(lme4) library(lme4)
Line 211: Line 224:
 (Intercept)  0.14156    0.03072   4.608 (Intercept)  0.14156    0.03072   4.608
 </code> </code>
 +
 Note that some of the default checks need to be switched off before ''lmer()'' will go through with fitting this model. Note that some of the default checks need to be switched off before ''lmer()'' will go through with fitting this model.
  
 A couple things are of note here: A couple things are of note here:
   - The estimated intercept differs from the model coefficient obtained earlier with the ''rma()'' function.   - The estimated intercept differs from the model coefficient obtained earlier with the ''rma()'' function.
-  - The standard error is also different (and hence, the test statistic and p-value also differs).+  - The standard error is also different (and hence, the test statistic and p-value also differ).
   - The estimated standard deviation of the random effects also differs (''0.0682'' for ''lme()'' and ''lmer()'' compared to ''0.0901'' for ''rma()'').   - The estimated standard deviation of the random effects also differs (''0.0682'' for ''lme()'' and ''lmer()'' compared to ''0.0901'' for ''rma()'').
  
-These discrepancies are due to the exact same reason described earlier. The ''lme()'' and ''lmer()'' functions assume that the sampling variances are not exactly known, but again just up to a proportionality constant, namely the residual variance.+These discrepancies arise for the same reason described earlier. In other words, the ''lme()'' and ''lmer()'' functions assume that the sampling variances are not exactly known, but again just up to a proportionality constant.
  
 To illustrate this, we can again factor in that constant into the sampling variances and refit the model with ''rma()'': To illustrate this, we can again factor in that constant into the sampling variances and refit the model with ''rma()'':
 +
 <code rsplus> <code rsplus>
 rma(yi, vi*sigma(res.lme)^2, data=dat) rma(yi, vi*sigma(res.lme)^2, data=dat)
Line 243: Line 258:
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 </code> </code>
 +
 Now these results match exactly what was obtained with the ''lme()'' function. Now these results match exactly what was obtained with the ''lme()'' function.
  
Line 249: Line 265:
 ==== Meta-Regression ==== ==== Meta-Regression ====
  
-The examples above show models without any covariates/predictors/moderators (i.e., the models only include an intercept term). However, the exact same discrepancy will be found when including covariates into the models. In this context, Thompson and Sharp (1999) also describe the difference discussed above. They denote the model fitted by ''lm()'' as a model that accounts for potential heterogeneity by a multiplicative factor (i.e., the residual variance), while the random-effects model fitted via the ''rma()'' function is a model that accounts for potential heterogeneity by an additive factor. They note that the "rationale for using a multiplicative factor for variance inflation is weak" (p. 2705) and clearly recommend to use a model with an additive heterogeneity component.+The examples above show models without any covariates/predictors/moderators (i.e., the models only include an intercept term). However, the exact same discrepancy will be found when including covariates in the models. In this context, Thompson and Sharp (1999) also describe the difference discussed above. They denote the model fitted by ''lm()'' as a model that accounts for potential heterogeneity by a multiplicative factor (i.e., the residual variance), while the random-effects model fitted via the ''rma()'' function is a model that accounts for potential heterogeneity by an additive factor. They note that the "rationale for using a multiplicative factor for variance inflation is weak" (p. 2705) and clearly recommend to use a model with an additive heterogeneity component.
  
 The possibility of using both an additive and multiplicative component simultaneously was not discussed by Thompson and Sharp (1999). This is in fact what the model fitted with the ''lme()'' function does. However, as far as I am aware of, no motivation supporting such a model has ever been given. In fact, it is my impression that those who use the ''lme()'' and ''lmer()'' functions for meta-analytic purposes are typically not aware of the issues discussed above. The possibility of using both an additive and multiplicative component simultaneously was not discussed by Thompson and Sharp (1999). This is in fact what the model fitted with the ''lme()'' function does. However, as far as I am aware of, no motivation supporting such a model has ever been given. In fact, it is my impression that those who use the ''lme()'' and ''lmer()'' functions for meta-analytic purposes are typically not aware of the issues discussed above.
Line 255: Line 271:
 ==== More Complex Models ==== ==== More Complex Models ====
  
-Note that the same issues apply when fitting more complex meta-analytic models (e.g., multivariate/multilevel models). The ''rma.mv()'' function in the metafor package can be used to fit such models (see the [[:analyses#multivariate_multilevel_meta-analysis_models|analysis examples]] section for some illustrations). Using the ''lme()'' and ''lmer()'' functions to fit such models will again lead to a mix of additive and multiplicative variance components, which is not how multivariate/multilevel meta-analytic models are typically defined.+Note that the same issue arises when fitting more complex meta-analytic models (e.g., multivariate/multilevel models). The ''rma.mv()'' function in the metafor package can be used to fit such models (see the [[:analyses#multivariate_multilevel_meta-analysis_models|analysis examples]] section for some illustrations). Using the ''lme()'' and ''lmer()'' functions to fit such models will again lead to a mix of additive and multiplicative variance components, which is not how multivariate/multilevel meta-analytic models are typically defined.
  
 ==== S-Plus Version of lme() ==== ==== S-Plus Version of lme() ====
  
-The ''nlme'' package was developed by José Pinheiro and Douglas Bates for both R and S-Plus. Interestingly, the S-Plus version has a special ''control'' argument that allows the user to set the error variance to a fixed constant. By setting this to 1, one can fit the exact same model as the ''rma()'' does:+The ''nlme'' package was developed by José Pinheiro and Douglas Bates for both R and S-Plus. Interestingly, the S-Plus version has a special ''control'' argument that allows the user to set the error variance to a fixed constant. By setting this to 1, one can fit the exact same model as the ''rma()'' function does: 
 <code rsplus> <code rsplus>
 res.lme <- lme(yi ~ 1, random = ~ 1 | study, weights = varFixed(~ vi), control=lmeControl(sigma = 1), data=dat) res.lme <- lme(yi ~ 1, random = ~ 1 | study, weights = varFixed(~ vi), control=lmeControl(sigma = 1), data=dat)
Line 289: Line 306:
 Number of Groups: 16 Number of Groups: 16
 </code> </code>
 +
 These are the exact same results as obtained earlier with the ''rma()'' function. Unfortunately, the R version of the ''nlme'' package does not provide this functionality. These are the exact same results as obtained earlier with the ''rma()'' function. Unfortunately, the R version of the ''nlme'' package does not provide this functionality.
  
-**Update:** The R version of the ''nlme'' package //does// allow the use of the ''lmeControl(sigma = 1)'' control argument (this was added in version 3.1-123, which was released 2016-01-17). However, using this does not yield the same results as obtained above (the results are close but not the same). The results given above have also been compared to those obtained with Stata (metareg command) and SAS (proc mixed) and are definitely correct. This issue has also been raised on the R-sig-ME mailing list ([[https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q2/024862.html|see here]]). See also this [[https://www.jepusto.com/bug-in-nlme-with-fixed-sigma/|blog post]] by James Pustejovsky for some further comparisons between ''lme()'' and the ''rma()'' and ''rma.mv()'' functions.+**Update:** The R version of the ''nlme'' package //does// allow the use of the ''lmeControl(sigma = 1)'' control argument (this was added in version 3.1-123, which was released 2016-01-17). However, using this does not yield the same results as obtained above (the results are close but not the same). The results given above (for S-Plus and ''rma()''have also been compared to those obtained with Stata (metareg command) and SAS (proc mixed) and are definitely correct. This issue has also been raised on the R-sig-ME mailing list ([[https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q2/024862.html|see here]]). See also this [[https://www.jepusto.com/bug-in-nlme-with-fixed-sigma/|blog post]] by James Pustejovsky for some further comparisons between the ''lme()''''rma()''and ''rma.mv()'' functions.
  
 ==== Summary ==== ==== Summary ====
tips/rma_vs_lm_lme_lmer.txt · Last modified: 2023/11/14 08:00 by Wolfgang Viechtbauer