Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
-----Original Message-----
From: Porter, Kyle [mailto:Kyle.Porter at osumc.edu]
Sent: Friday, August 25, 2017 18:17
To: Porter, Kyle; Viechtbauer Wolfgang (SP); 'r-sig-meta-analysis at r-project.org'
Subject: RE: confidence interval for I^2 in multilevel model
I think I may have ended up with an over-specified model. The result was an estimate of zero for sigma^2.2 for the Study/EffectSizeEntry factor. (EffectSizeEntry is the row number in my data; each effect size was given a unique identifier. It equates to each Study/Bethesda combination). My thought was this specification was analogous to the three-level in Konstantopoulos with my EffectSizeEntry equating to his Study level and my Study equating to his District level.
mvregML <- rma.mv(data=PPV_prevdat, yi, vi, random = ~ 1 | Study/EffectSizeEntry, mods = ~ Bethesda + prev + Blinded + NRAS12_13 + HRAS12_13 + KRAS61)
Even if I take out the moderators I get the zero estimate for sigma^2.2 for the Study/EffectSizeEntry factor.
mvregML <- rma.mv(data=PPV_prevdat, yi, vi, random = ~ 1 | Study/EffectSizeEntry)
-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Porter, Kyle
Sent: Friday, August 25, 2017 11:40 AM
To: 'Viechtbauer Wolfgang (SP)'; 'r-sig-meta-analysis at r-project.org'
Subject: Re: [R-meta] confidence interval for I^2 in multilevel model
Wolfgang,
Thank you for the response. I should have provided more information about the specifics of the meta-analysis we are working on. The situation is different than the three-level scenario in Konstantopoulos; we have multiple levels of results within study rather than results from multiple studies that are clustered together such as "district" in Konstantopoulos. I was attempting to account for the correlation in the results within study (in my code "ID" represents the study). I changed "ID" to "Study" in the message below for clarity.
Our outcome is a single proportion (positive predictive value from the number of true positives (TP) and positive tests (TP + FP, designated PosTest)) and I am using the arcsin transformation.
For each study, we have separated the effect sizes into up to 3 categories based on nodule classification (Bethesda) as either 3, 4, or 5. So the full dataset has up to 3 rows of data per study with an effect size (PPV) and an indicator (Bethesda) for each specific nodule classification level.
Sample data (first 7 studies; moderator/covariates not shown) Study,Bethesda,TP,PosTest
1,3,4,5
1,4,4,5
1,5,2,2
2,4,3,5
3,3,2,3
3,5,10,10
4,5,1,1
5,4,22,31
6,4,2,15
7,4,12,21
7,5,1,1
I have run separate analyses for each Bethesda level to evaluate the proportional of residual heterogeneity. I was able to get I^2 and the corresponding confidence intervals for these models:
mregB3 <- rma.uni(measure="PAS", data=B3, xi=TP, ni=posTest, method="REML", mods = ~ factor(Blinded) + prev + factor(HRAS12_13) + factor(KRAS61) + factor(NRAS12_13))
mregB3
confint(mregB3)
where B3 is a data frame subset containing only Bethesda level 3 nodules. I did the same for 4 and 5 separately.
We would also like to estimate I^2 with confidence intervals using all classifications combined (3, 4, 5). Now we have multiple data points per study, so my thinking was I need to account for potential correlation among these within-study effect sizes. Having thought further about this since my initial post, I wonder if the existing study-level random effect may be sufficient along with the indicator variables for the Bethesda factor to account for the within-study correlation. But if I run:
mvregAll2 <- rma(measure="PAS",data=PPVdat, xi=TP, ni=posTest, method="REML", mods = ~ factor(Bethesda) + factor(Blinded) + factor(NRAS12_13) + factor(HRAS12_13) + factor(KRAS61))
then the model is treating each effect size independently, as if it were from a separate study.
So back to my original attempt:
PPVdat <- escalc(measure="PAS", data=RAS_PPV, xi=TP, ni=posTest) mvregAll <- rma.mv(data=PPVdat, yi, vi, method="REML", random = ~ 1 | Study, mods = ~ factor(Bethesda) + factor(Blinded) + factor(NRAS12_13) + factor(HRAS12_13) + factor(KRAS61))
The results of mvregAll2 and mvregAll are similar, but I think mvregAll does have the clustered relationship I am looking to include, but perhaps the random effect is conflating study level variation with "sub-study" level variation (variation across each row of data). Perhaps I need to add a variance component at the effect size level (as if they were separate studies) and then also the study level variance component and specify a model in the three-level form.
I will attempt that and report back. Any other feedback is appreciated in the meantime. (Also thanks for pointing out that struct="DIAG" was not necessary and not functional in this model; I am a SAS user with much less R experience so my "translation" of model specification was off).
Thanks,
Kyle