[R-meta] metafor package in R - Risk ratios using rma.mv()
Viechtbauer, Wolfgang (SP)
wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Fri Jul 26 22:44:02 CEST 2019
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.
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
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.
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').
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.
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
data <- read.xlsx('702_tb_incidence.xlsx', 1)
#calculating the incidence rate and the variance
pd_ec <- escalc(
measure = 'IR',
xi = data$n_diagnosed,
ti = data$person_years,
data = data
#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,
Multivariate Meta-Analysis Model (k = 8; method: REML)
logLik Deviance AIC BIC AICc
15.4961 -30.9922 -20.9922 -24.0607 39.0078
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
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!
More information about the R-sig-meta-analysis