[R-meta] confidence interval for I^2 in multilevel model

[Viechtbauer Wolfgang (SP)]

If I understand you correctly, multiple effects within each study are calculated based on different subjects. In that case, this is exactly like the three-level model in Konstantopoulos. In Konstantopoulos, 'studies' are clustered in 'districts', but that is no different than 'nodule classification' clustered in 'studies' (assuming that there is no overlap in the subjects used in computing the outcomes for the different classificatons within a study).

So, an appropriate model would be the one with 'random = ~ 1 | Study/EffectSizeEntry'. It could very well be the case that all of the variance is then at the Study level. That's fine.

But you could actually go a step further, since the different levels of 'Bethesda' are clearly distinct (that part is a bit different than the Konstantopoulos example, where the numbering of studies within districts is entirely arbitrary). So, you could try:

mvregML <- rma.mv(data=PPV_prevdat, yi, vi, random = ~ factor(Bethesda) | Study, struct="UN", mods = ~ Bethesda + prev + Blinded + NRAS12_13 + HRAS12_13 + KRAS61)

This allows for different amounts of heterogeneity in the three levels of Bethesda and different pairwise correlations.

Best,
Wolfgang

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