Message-ID: <dd2497e1d608440ca57073ed8fe37225@UM-MAIL3214.unimaas.nl>
Date: 2019-10-05T21:49:17Z
From: Wolfgang Viechtbauer
Subject: [R-meta] metafor package in R - Risk ratios using rma.mv()
In-Reply-To: <CA+N-+=cpjGa7t1_YKRRze_V4xkcrhurkoVxQJhXrWJsNzENMcw@mail.gmail.com>
Dear Olivia and Leo,
rma.glmm() does Poisson and logistic models for incidences and proportions, respectively, but it cannot fit models with multilevel strutures (which is what 'random = ~ 1 | study_id/cohort_id' implies you are dealing with). You can however fit such models directly using glmer() from the 'lme4' package. Alternatively, you can stick to rma.mv() (which does allow you to account for the multilevel structure as you are doing right now) and switch to analyzing log incidence rates (IRLN) and logit transformed proportions (PLO) (note: for proportions, the logit transformation is usually better than using log-transformed proportions).
Once you fit such models and then apply the respective back-transformations (exp for IRLN, transf.ilogit for PLO), you are guaranteed to get non-negative estimates and confidence interval bounds. That's just the nature of these transformations. For example, exp() will give you something that is >=0 no matter what value you plug in.
Best,
Wolfgang
-----Original Message-----
From: Leo Martinez [mailto:leomarti at stanford.edu]
Sent: Saturday, 05 October, 2019 0:28
To: Viechtbauer, Wolfgang (SP)
Cc: Olivia Cords; r-sig-meta-analysis at r-project.org
Subject: Re: [R-meta] metafor package in R - Risk ratios using rma.mv()
ATTACHMENT(S) REMOVED: 1003_PR_IR.xlsx
Dear Wolfgang and All,?
Thanks so much for your prior help.?
We have calculated incidence and prevalence rates from a mixed-effects model using the?rma.mv?command. We are attaching the results below.
Variable??????? Cohorts (n)???? Incidence Rate/100k???? Lower 95% CI??? Upper 95% CI
World Health Organization Region??????? ??????? ??????? ???????
??? Americas??? 225???? 324.5?? 166.7?? 482.2
??? African???? 30????? 1906.4? 1134.1? 2678.7
??? Eastern Mediterranean?????? 24????? 249.1?? 32.3??? 466.0
??? European??? 33????? 767.6?? 407.8?? 1127.3
??? South-East Asia???? 48????? 1148.7? 628.6?? 1668.9
??? Western Pacific???? 30????? 560.0?? -131.7? 1251.7
??????? Cohorts (n)???? Prevalence Rate Lower 95% CI??? Upper 95% CI
World Health Organization Region??????? ??????? ??????? ???????
??? Americas??? 33????? 1.7???? 0.9???? 2.5
??? African???? 31????? 2.8???? 1.4???? 4.1
??? Eastern Mediterranean?????? 7?????? 1.9???? 0.5???? 3.4
??? European??? 21????? 1.9???? 0.8???? 3.0
??? South-East Asian??? 10????? 2.4???? -1.1??? 6.0
??? Western Pacific???? 54????? 1.2???? -0.1??? 2.4
Unfortunately we have some negative confidence intervals for some of our incidence and prevalence estimates. We would like to not have any negative confidence intervals?and therefore would like to switch the models that we are using.
Is there a way to keep our code (which we have put below for both incidence and prevalence) and run a poisson model for incidence and a binary or beta model for prevalence so that we no longer have a negative confidence interval for some of our variables? We noticed that when we run the model using log transformed incident and prevalence rates, the confidence intervals are positive. We are also wondering what the difference is between using PR/IR versus PLN/IRLN for fitting the model, and why the latter would result in all positive confidence intervals.
Thank you again for all your help!?
Best
Olivia and Leo
Code for incidence rates:?
#data subsetted by WHO region
pd_ec <- escalc(measure = 'IR', xi = data_sub$inc_positive,ti = data_sub$inc_person_years, append = TRUE,?
? ? ? ? ? ? data = data_sub)?
m0 <-?rma.mv(yi, vi, method='REML', mods = formula,
? ? ? ? ? ? ? ? ? ? ? ? random= ~ 1 | study_id/cohort_id,
? ? ? ? ? ? ? ? ? ? ? ? tdist=TRUE,
? ? ? ? ? ? ? ? ? ? ? ? data=pd_ec)
Code for prevalence rates:?
#data subsetted by WHO region
pd_ec <- escalc(
? ? ? ? ? ? measure = 'PR', xi = data_sub$prev_positive,ni = data_sub$prev_total_n, append = TRUE,
? ? ? ? ? ? data = data_sub)
?m0 <-?rma.mv(yi, vi, method='REML', mods = formula,
? ? ? ? ? ? ? ? ? ? ? ? random= ~ 1 | study_id/cohort_id,
? ? ? ? ? ? ? ? ? ? ? ? tdist=TRUE,
? ? ? ? ? ? ? ? ? ? ? ? data=pd_ec)
Thank you for any help you can provide!
Best
Leo and Olivia
Leonardo Martinez, PhD, MPH
Stanford University School of Medicine
Division of Infectious Diseases and Geographic Medicine
300 Pasteur Drive, Lane Building, Stanford, CA 94305
Phone: +1.202.769.8090 ?
Email:?leomarti at stanford.edu;?chopotin at gmail.com
Website