[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,
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>
```