[R-meta] Mismatch between output from sub-group analysis and forest plot
Joao Afonso
jot@|on@o @end|ng |rom gm@||@com
Thu Feb 13 11:47:16 CET 2020
Dear Gerta,
That seems to be again a error from my side as I forgot to provide the
code for the double arcsine transformation of the data. Once again I
apologise for all these mistakes (I admit I am new to this platform
and also to R). Let me one last time try to summarise everything and
also provide the data. The problem I am finding is a mismatch between
the pooled estimate when running the sub-group analysis on a moderator
and the figure displayed in the resulting forest plot.
The data follows below and attached (.cvs).
author ssizeanimal nlameanimal
<chr> <dbl> <dbl>
1 Griffiths et al., 2018 14700 4145
2 Phillips, 1990 162 98
3 Chaplin et al., 2000 178 13
4 Barker et al., 2010 33415 12297
5 Archer et al., 2010 1400 868
6 Amory et al., 2008 1824 636
7 Manning, 2018 224 73
8 Bell et al., 2013 332 40
9 Walker et al., 2008 36 18
10 Weaver, 1997 7700 2310
# ... with 23 more rows
The data was transformed using the double arcsine method as per code below:
ies.da=escalc(xi=nlameanimal, ni=ssizeanimal,
data=prevalence_2020_cow_nomv, measure="PFT", add=0)
pes.da=rma(yi, vi, data=ies.da, method = "DL")
pes=predict(pes.da, transf=transf.ipft.hm,
targ=list(ni=prevalence_2020_cow_nomv$ssizeanimal))
print(pes.da, digits=3)
confint(pes.da, digits=3)
print(pes, digits=3)
As heterogeneity was high sub-group analysis followed. For moderator
"lcmbi" the code used was the following:
ies.da.lcmbi<-ies.da %>%
filter(is.na(lcmbi)==FALSE)
subganal.lcmbi=rma(yi, vi, data=ies.da.lcmbi, mods=~lcmbi, method="DL")
pes.da.lcm=rma(yi, vi, data=ies.da.lcmbi, subset = lcmbi=="Records",
method="DL")
pes.da.records=rma(yi, vi, data=ies.da.lcmbi, subset =
lcmbi=="Locomotion Scoring Method", method="DL")
pes.subg.lcmbi=predict(subganal.lcmbi,
transf=transf.ipft.hm,
targ=list(ni=prevalence_2020_cow_nomv$ssizeanimal))
dat.samevar=data.frame(estimate=c((pes.da.lcm$b)[1], (pes.da.records$b)[1]),
stderror=c((pes.da.lcm$se)[1], (pes.da.records$se)[1]),
tau2=subganal.lcmbi$tau2)
pes.da.lcmbi=rma(estimate, sei=stderror,
method="DL",
data=dat.samevar)
pes.lcmbi=predict(pes.da.lcmbi, transf=transf.ipft.hm,
targ=list(ni=prevalence_2020_cow_nomv$ssizeanimal))
print(pes.da.lcm, digits=3) #display subgroup 1 summary effect size
print(pes.da.records, digits=3) #display subgroup 2 summary effect size
print(subganal.lcmbi, digits=3) #display subgroup analysis results
print(pes.lcmbi, digits=3) #display recomputed summary effect size
Which produces the following output:
Random-Effects Model (k = 11; tau^2 estimator: DL)
tau^2 (estimated amount of total heterogeneity): 0.027 (SE = 0.021)
tau (square root of estimated tau^2 value): 0.164
I^2 (total heterogeneity / total variability): 99.93%
H^2 (total variability / sampling variability): 1428.81
Test for Heterogeneity:
Q(df = 10) = 14288.143, p-val < .001
Model Results:
estimate se zval pval ci.lb ci.ub
0.506 0.051 9.905 <.001 0.406 0.606 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> print(pes.da.records, digits=3) #display subgroup 2 summary effect size
Random-Effects Model (k = 22; tau^2 estimator: DL)
tau^2 (estimated amount of total heterogeneity): 0.011 (SE = 0.007)
tau (square root of estimated tau^2 value): 0.104
I^2 (total heterogeneity / total variability): 99.23%
H^2 (total variability / sampling variability): 130.14
Test for Heterogeneity:
Q(df = 21) = 2732.942, p-val < .001
Model Results:
estimate se zval pval ci.lb ci.ub
0.600 0.024 25.359 <.001 0.553 0.646 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> print(subganal.lcmbi, digits=3) #display subgroup analysis results
Mixed-Effects Model (k = 33; tau^2 estimator: DL)
tau^2 (estimated amount of residual heterogeneity): 0.022 (SE = 0.013)
tau (square root of estimated tau^2 value): 0.148
I^2 (residual heterogeneity / unaccounted variability): 99.82%
H^2 (unaccounted variability / sampling variability): 549.07
R^2 (amount of heterogeneity accounted for): 29.36%
Test for Residual Heterogeneity:
QE(df = 31) = 17021.085, p-val < .001
Test of Moderators (coefficient 2):
QM(df = 1) = 2.899, p-val = 0.089
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.601 0.033 18.449 <.001 0.537 0.665 ***
lcmbiRecords -0.096 0.056 -1.703 0.089 -0.207 0.015 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> print(pes.lcmbi, digits=3) #display recomputed summary effect size
pred ci.lb ci.ub cr.lb cr.ub
0.283 0.206 0.368 0.168 0.415
And to produce the forest plot I use the following code:
subganal.lcmbi=rma(yi, vi, data=ies.da.lcmbi, mods=~lcmbi, method="DL")
pes.summary=metaprop(nlameanimal, ssizeanimal, author, data=ies.da.lcmbi,
sm="PFT",
method.tau="DL",
method.ci="NAsm",
byvar=lcmbi,
tau.common=TRUE,
tau.preset=sqrt(subganal.lcmbi$tau2))
precision=sqrt(ies.da.lcmbi$vi)
jpeg(file="forest_cow_cow_lcmbi.jpeg", width=11, height=12, units="in", res=300)
forest(pes.summary,
xlim=c(0,1),
pscale=1,
rightcols=FALSE,
leftcols=c("studlab", "effect", "ci"),
leftlabs=c("Study", "Prevalence", "95% C.I."),
text.random="Prevalence of Lameness in British Dairy Cattle",
bylab = "Lameness Data Source",
print.byvar = TRUE,
xlab="Prevalence of Lameness", smlab="",
weight.study="random", squaresize=0.5, col.square="navy",
col.diamond="maroon", col.diamond.lines="maroon",
pooled.totals=FALSE,
comb.fixed=FALSE,
fs.hetstat=10,
print.tau2=TRUE,
print.Q=TRUE,
print.pval.Q=TRUE,
print.I2=TRUE,
digits=2,
sortvar = precision)
dev.off()
The forest plot follows attached.
Once again I am truly sorry for all this back and forth. I haven't
made you helping me any easier. But I would be very grateful for any
inputs solving this issue.
On Wed, Feb 12, 2020 at 11:02 PM Dr. Gerta Rücker
<ruecker using imbi.uni-freiburg.de> wrote:
>
> Dear Joao,
>
> Honestly, I don't know. Comparing the code for the forest plot with the
> rest of the code, I find it difficult to identify corresponding
> variables. For example, in the forest plot code the number of events is
> called nlameanimal, but this variable name occurs nowhere in the rest of
> the code where you always use yi. Moreover, you have two subsets,
> "Records" (11 studies) and "Locomotion Scoring Method" (22 studies), but
> it seems they are confused somewhere. I find the code extremely
> difficult to understand without knowing the data.
>
> Best,
>
> Gerta
>
> Am 12.02.2020 um 18:05 schrieb Joao Afonso:
> > Dear Gerta,
> >
> > Thank you so much for all the insights and guidance. I was looking at
> > the figures and there is still one thing I can't make sense of which
> > relates to this part of the code:
> >
> >> print(pes.lcmbi, digits=3) #display recomputed summary effect size
> > This hopefully is the part of the code that provide the
> > back-transformed pooled estimate of prevalence as it outputs the
> > following:
> >
> > pred ci.lb ci.ub cr.lb cr.ub
> > 0.283 0.206 0.368 0.168 0.415
> >
> > I would expect to see this in the forest plot, yet what the plot provides is
> >
> > 0.29 (0.24-0.34)
> >
> > Any thoughts?
> >
> > Many thanks once again. Have a great evening,
> >
> > On Wed, Feb 12, 2020 at 5:00 PM Gerta Ruecker
> > <ruecker using imbi.uni-freiburg.de> wrote:
> >> Dear Joao,
> >>
> >> This is what I suspected in one of my previous e-mails. It is all right. While metafor gives the results transformed, that is, on the arcsin(sqrt())-scale, the forest plot provides results conveniently on the original probability scale (backtransformed). You can (roughly) switch between both measures by using
> >>
> >> arcsin(sqrt(x)) (from the forest plot to the text output), or vice versa
> >> sin^2(x) (from the text output to the forest plot).
> >>
> >> It is only rough here (1) because the package used the Freeman-Tukey transformation which is more complicated and (2) because of rounding error.
> >>
> >> See some calculations inline below.
> >>
> >> Best,
> >>
> >> Gerta
> >>
> >> Am 12.02.2020 um 11:47 schrieb Joao Afonso:
> >>
> >> Dear all,
> >>
> >> Once again I am sorry for posting this message again but I forgot to
> >> send the outputs.
> >>
> >> Again I am conducting a meta-analysis on prevalence and incidence
> >> data. As there is plenty of heterogeneity I am doing a sub-group
> >> analysis. However I
> >> am finding the output of the same and forest plots produced to have
> >> different values for the pooled estimates. I am using the double
> >> arcsine transformation and then have the data back-transformed to
> >> produce the estimates.
> >>
> >> [...]
> >>
> >> Random-Effects Model (k = 11; tau^2 estimator: DL)
> >>
> >> tau^2 (estimated amount of total heterogeneity): 0.027 (SE = 0.021)
> >> tau (square root of estimated tau^2 value): 0.164
> >> I^2 (total heterogeneity / total variability): 99.93%
> >> H^2 (total variability / sampling variability): 1428.81
> >>
> >> Test for Heterogeneity:
> >> Q(df = 10) = 14288.143, p-val < .001
> >>
> >> Model Results:
> >>
> >> estimate se zval pval ci.lb ci.ub
> >> 0.506 0.051 9.905 <.001 0.406 0.606 ***
> >>
> >> This apparently corresponds to the smaller subgroup where the forest plot shows 0.23 [0.0.16; 0.31], but gives the untransformed result on the arcsin(sqrt()) scale. In fact:
> >>
> >> arcsin(sqrt(0.23 [0.16; 0.31]) = 0.5 [0.41; 0.59]
> >>
> >> [...]
> >>
> >>
> >> Random-Effects Model (k = 22; tau^2 estimator: DL)
> >>
> >> tau^2 (estimated amount of total heterogeneity): 0.011 (SE = 0.007)
> >> tau (square root of estimated tau^2 value): 0.104
> >> I^2 (total heterogeneity / total variability): 99.23%
> >> H^2 (total variability / sampling variability): 130.14
> >>
> >> Test for Heterogeneity:
> >> Q(df = 21) = 2732.942, p-val < .001
> >>
> >> Model Results:
> >>
> >> estimate se zval pval ci.lb ci.ub
> >> 0.600 0.024 25.359 <.001 0.553 0.646 ***
> >>
> >> Again: arcsin(sqrt(0.32 [0.26; 0.38]) = 0.60 [0.54; 0.66] (similar to what metafor shows, differences due to rounding and the Freeman-Tukey transformation).
> >>
> >> Analogous for the whole group.
> >>
> >> [...]
> >>
> >>
> >> --
> >>
> >> Dr. rer. nat. Gerta Rücker, Dipl.-Math.
> >>
> >> Institute of Medical Biometry and Statistics,
> >> Faculty of Medicine and Medical Center - University of Freiburg
> >>
> >> Stefan-Meier-Str. 26, D-79104 Freiburg, Germany
> >>
> >> Phone: +49/761/203-6673
> >> Fax: +49/761/203-6680
> >> Mail: ruecker using imbi.uni-freiburg.de
> >> Homepage: https://www.imbi.uni-freiburg.de/persons/ruecker/person_view
> >
> >
--
João Afonso
DVM, MSc Veterinary Epidemiology
PhD Student
Department of Infection and Global Health
University of Liverpool
+351914812305
-------------- next part --------------
A non-text attachment was scrubbed...
Name: prev_cow.cvs
Type: application/octet-stream
Size: 1061 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20200213/7ef204e4/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: forest_cow_cow_lcmbi_brah.jpeg
Type: image/jpeg
Size: 144106 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20200213/7ef204e4/attachment-0001.jpeg>
More information about the R-sig-meta-analysis
mailing list