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

Olivia Cords ocord@ @end|ng |rom @|umn|@@t@n|ord@edu
Thu Aug 8 00:08:15 CEST 2019


Dear Wolfgang,
Thanks so much again for your help. The comments you provided have been
incredibly useful and I have implemented your suggestions. I just wanted to
follow up with you on this topic with a few additional clarifications
regarding incident rates, incident rate ratios, confidence intervals, and
similar issues for prevalence/odds ratios.

*Incidence Rates*
I calculated the incidence rates (measure = "IR", ti =
data$person_years/1000) for tuberculosis stratified by World Health
Organization Region and got the following results from the model.

m0 <- rma.mv(yi, vi, method='REML', mods = formula,
                        random= ~ 1 | study_id/cohort_id,
                        tdist=TRUE,
                        data=pd_ec)

                                          estimate      se     tval    pval
intrcpt                                     3.5210  0.9728   3.6195  0.0003
who_region.1African Region                 14.5115  2.6764   5.4219  <.0001
who_region.1Eastern Mediterranean Region   -0.9006  3.8084  -0.2365  0.8132
who_region.1European Region                 5.2299  2.2226   2.3531  0.0191
who_region.1South-East Asian Region         7.2904  3.7296   1.9547  0.0513
who_region.1Western Pacific Region          2.1339  3.3186   0.6430  0.5206
                                            ci.lb    ci.ub
intrcpt                                    1.6083   5.4337  ***
who_region.1African Region                 9.2492  19.7738  ***
who_region.1Eastern Mediterranean Region  -8.3885   6.5873
who_region.1European Region                0.8600   9.5998    *
who_region.1South-East Asian Region       -0.0426  14.6235    .
who_region.1Western Pacific Region        -4.3910   8.6587

  To find the incidence rates, I added the slopes to the intercept to
produce the following Incident Rates/1000 person years:
Intercpt (Americas) 3.521
African Region 18.03249
Eastern Mediterranean Region 2.620433
European Region 8.750917
South-East Asian Region 10.81145
Western Pacific Region 5.654908

To confirm these results I also ran the model on the dataset subsetted by
region, using the 'estimate' as each region's incidence rate. This resulted
in the following:

Region of the Americas
3.244745079
African Region
19.06.383209
Eastern Mediterranean Region 2.491278082
European Region 7.675655135
South-East Asian Region
11.48738471
Western Pacific Region 5.600276079

*Q1) Why are the Incidence rates different when using the entire dataset
versus first subsetting by region? *

*Rate Ratios:*

Following your suggestion above for calculating rate ratios for each WHO
Region, I used the measure = "IRLN" and exponentiated the coefficients of
the model. I got the following model output:

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

Variance Components:

            estim    sqrt  nlvls  fixed              factor
sigma^2.1  2.3898  1.5459     76     no            study_id
sigma^2.2  0.5833  0.7637    390     no  study_id/cohort_id

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

Test of Moderators (coefficients 2:6):
F(df1 = 5, df2 = 384) = 7.3370, p-val < .0001

Model Results:

                                          estimate      se      tval    pval
intrcpt                                    -7.1326  0.2449  -29.1259  <.0001
who_region.1African Region                  3.1798  0.6332    5.0217  <.0001
who_region.1Eastern Mediterranean Region    0.9541  0.9634    0.9904  0.3226
who_region.1European Region                 1.8665  0.5612    3.3261  0.0010
who_region.1South-East Asian Region         2.6693  0.9360    2.8518  0.0046
who_region.1Western Pacific Region          1.4697  0.8381    1.7536  0.0803
                                            ci.lb    ci.ub
intrcpt                                   -7.6141  -6.6511  ***
who_region.1African Region                 1.9348   4.4248  ***
who_region.1Eastern Mediterranean Region  -0.9401   2.8483
who_region.1European Region                0.7631   2.9698  ***
who_region.1South-East Asian Region        0.8290   4.5096   **
who_region.1Western Pacific Region        -0.1781   3.1176    .


And exponentiating the coefficients, I got these rate ratios:

Intrcpt (Americas) 0.000799
African Region 24.04225
Eastern Mediterranean Region 2.596419
European Region 6.465383
South-East Asian Region 14.42971
Western Pacific Region 4.348128

Based on the incidence rates produced by using the entire dataset in the
model and first subsetting by region, these Rate Ratios don't seem to be
correct. Simply dividing the incidence rates by a comparator (Region of the
Americas) to produce rate ratios would give the following:
Entire Dataset First Subsetted
Region of the Americas
African Region 5.121411531 5.875294245
Eastern Mediterranean Region 0.744106788 0.767788539
European Region 2.485089463 2.365564921
South-East Asian Region 3.070434536 3.540304225
Western Pacific Region 1.604657768 1.725952561


*Q2) Is exponentiating the coefficients (measure = IRLN) the way to
calculate rate ratios? Why are these results so different?*

*Confidence Intervals*

To calculate the 95% confidence intervals for the rate ratios, I first
calculated the standard deviation (SD[ln(IR)] = (1/A1 + 1/A2)^0.5, where A1
and A2 are the number of tuberculosis cases in each region), and then used
the following: 95% CI's = exp[ln(IR) ± 1.96(SD)]).

*Q3) Is this the right approach to calculating the confidence intervals? Is
there a way to calculate the confidence intervals from the model (either
measure = "IRLN" or measure = "IR") output?*

*Prevalence:*
Additionally, we have data for calculating tuberculosis prevalence rate by
WHO region and are trying to calculate odds ratios. For prevalence rate,
all of the code for the model is the same, except the escalc function is as
follows:
            pd_ec <- escalc(measure = 'PR', xi = data_sub$prev_positive,ni
= data_sub$prev_total_n, append = TRUE,
            data = data_sub),

*Q4) For prevalence rates, should we first subset the data to find each
region's prevalence rate, or use the entire dataset (adding intercept to
each region's slope)?*
*Q5)How would we calculate odds ratios using the model?*
*Q6)How would we calculate the confidence intervals for these odds ratios
using the model? *

Thank you so much again for your advice. It's really been indispensable and
is greatly appreciated.

All the best,
Olivia

On Fri, Jul 26, 2019 at 1:44 PM Viechtbauer, Wolfgang (SP) <
wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

> 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
>

	[[alternative HTML version deleted]]



More information about the R-sig-meta-analysis mailing list