[R-meta] metafor package in R - Risk ratios using rma.mv()

[Viechtbauer, Wolfgang (SP)]

Hi Olivia,

A couple comments before I get to your actual question:

1) You might want to divide the person-years by some factor. For example:

ti = data$person_years / 1000

The incidence rates then reflect the number of cases per 1000 person years. That way, the rates aren't so small and the same goes for the model coefficients (e.g., the 0.0027 will turn into 2.7, which is a bit easier to communicate).

2) You don't need to write pd_ec$yi and pd_ec$vi when you use data=pd_ec. So, rma.mv(yi, vi, ..., data=pd_ec) is sufficient.

3) You should include random effects for studies (as you have done) and also random effects for estimates within studies. So:

pd_ec$id <- 1:nrow(pd_ec)

random = ~ 1 | study_id / id

In fact, I see that you may have multiple estimates for the same cohort, so one could even use an additional level, that is:

random = ~ 1 | study_id / cohort_id / id

But if you only have 8 rows of data, then I wouldn't do that (so then I would use random = ~ 1 | study_id / id).

4) Now for your actual question: You are analyzing incidence rates, so the coefficients of your model reflect how the (average) rate changes as a function of those moderators (e.g., while the rate is 0.0027 for "1Americas", the rate is 0.0027-0.0025=0.0002 for "Africa"). In other words, the coefficients represent incidence rate differences. If you want ratios, you should analyze log incidence rates (measure="IRLN"). Then the coefficients reflect differences between log incidence rates, so if you exponentiate the coefficients of your model, you get incidence rate ratios.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Olivia Cords
Sent: Friday, 26 July, 2019 21:27
To: r-sig-meta-analysis using r-project.org
Subject: [R-meta] metafor package in R - Risk ratios using rma.mv()

ATTACHMENT(S) REMOVED: metafor_output_0724 (1).png | 702_tb_incidence.xlsx | 702_metafor_tb_incidence (4).r 

Hello!
My name is Olivia and I'm a researcher at Stanford University. Our group is trying to calculate relative risk ratios for the a meta-analysis extracting  incidence rates of tuberculosis disease using the metafor package.

Study Design:
Incidence rates were extracted for each study identified by a systematic review, with some studies reporting multiple rates for different years or locations. In this case, multiple rates were treated as different study cohorts, meaning that the data is clustered ('cohort_id' nested within 'study_id').

Model:
We used the rma.mv() function, inputting the calculated incidence rate as the outcome variable ('pdc$yi), the variance ('pdc$vi'), WHO region ('who_region') and whether the study was conducted through passive or active screening (passive_active') as moderators, and a random effects argument for the study level ('study_id'). We are unclear how to go from the output to Risk Ratios.

Data:
study_id    cohort_id   n_diagnosed person_years   who_region  passive_active
131         34          77          14298          1Americas   1Passive
93          120         5           27750          1Americas   1Passive
93          121         14          277150         1Americas   1Passive
93          122         15          2000           1Americas   1Passive
136         383         2           2000           Africa      2Active
136         383         7           100000         Africa      2Active
187         16          16          100000         Africa      3Not Specified
187         517         2           100000         S.E. Asia   3Not Specified

Code:

library(xlsx)
library(metafor)

data <- read.xlsx('702_tb_incidence.xlsx', 1)
print(data)
#calculating the incidence rate and the variance
pd_ec <- escalc(
   measure = 'IR',
   xi = data$n_diagnosed,
   ti = data$person_years,
   data = data
)
pd_ec

#specifying mixed effects model
#first level cohort incidence rate and variation
#second level study_id
#who_region and passive_active as moderators
m0 <- rma.mv(pd_ec$yi, pd_ec$vi, method='REML', mods = ~ who_region + passive_active,
        random= ~ 1 | study_id,
        tdist=TRUE,
        data=pd_ec)

summary(m0)

Output:
Multivariate Meta-Analysis Model (k = 8; method: REML)

  logLik  Deviance       AIC       BIC      AICc
 15.4961  -30.9922  -20.9922  -24.0607   39.0078  

Variance Components:

            estim    sqrt  nlvls  fixed    factor
sigma^2    0.0000  0.0037      4     no  study_id

Test for Residual Heterogeneity:
QE(df = 4) = 94.4465, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 4) = 3.7931, p-val = 0.1152

Model Results:

                       estimate      se     tval    pval    ci.lb   ci.ub
intrcpt                  0.0027  0.0027   1.0073  0.3708  -0.0047  0.0101   
who_regionAfrica        -0.0025  0.0046  -0.5493  0.6120  -0.0153  0.0102   
who_regionS.E. Asia     -0.0027  0.0046  -0.5797  0.5932  -0.0154  0.0101   
passive_active2Active   -0.0001  0.0053  -0.0167  0.9874  -0.0148  0.0146 

Any advice/insight in much appreciated!

Best wishes,
Olivia