# I'm trying to adapt to my own data the double arcsine transformation
according to Miller (1978) as described in the metafor site
# I've understood all passages (I think)
# I'm able to build a forest plot with the correct data
# I would like to build a funnel plot and a radial plot
# However, the rule that is helpful to build a forest plot does not work
for the radial or the funnel plot
# When I use the results of the fixed (or random) effects model, as
expected, the estimates are wrong (about two times the correct estimates)
# How can the radial and funnel plot be built for the double arcsine
transformation?
# Thank you in advance, Antonello
# Here the code I've used
library(metafor)
# The data used by Miller (1978) to illustrate the transformation and its
inversion can be re-created with:
dat <- data.frame(xi=c(3, 6, 10, 1), ni=c(11, 17, 21, 6))
dat$pi <- with(dat, xi/ni)
dat <- escalc(measure="PFT", xi=xi, ni=ni, data=dat)
dat
# The yi values are the Freeman-Tukey (double arcsine) transformed
proportions,
# while the vi values are the corresponding sampling variances.
# We can check whether the back-transformation works for individual values
with:
transf.ipft(dat$yi, dat$ni)
# As described by Miller (1978), we can aggregate the transformed values,
# either by computing an unweighted or a weighted mean (with
inverse-variance weights).
# The unweighted mean can be obtained with:
res <- rma(yi, vi, method="FE", data=dat, weighted=FALSE)
res
# The value reported by Miller for the unweighted mean is again twice as
large as the value given above,
# which we can confirm with:
predict(res, transf=function(x) x*2)
# To back-transform the average, a value for the sample size is needed.
# Miller suggests to use the harmonic mean of the individual sample sizes
in the inversion formula.
# This can be accomplished with:
predict(res, transf=transf.ipft.hm, targs=list(ni=dat$ni))
# Since the true proportions appear to be homogeneous (e.g., Q(3)=2.18,
p=.54),
# a more efficient estimate of the true proportion can be obtained by using
inverse-variance weights.
# For this, we first synthesize the transformed values with:
res <- rma(yi, vi, method="FE", data=dat)
res
# Again, the value reported by Miller is twice as large.
# We can back-transform the estimated transformed average with:
predict(res, transf=transf.ipft.hm, targs=list(ni=dat$ni))
# When using the Freeman-Tukey transformation,
# an additional complication arises when using the back-transformation.
# To back-transform the individual proportions, we need the individual
sample sizes:
transf.ipft(dat$yi, dat$ni)
# To back-transform the estimated average transformed proportion,
# we need to use the harmonic mean of the sample sizes:
res <- rma(yi, vi, method="FE", data=dat)
pred <- predict(res, transf=transf.ipft.hm, targs=list(ni=dat$ni))
pred
# To build a forest plot we need to first obtain the CI bounds of the
individual studies with:
dat.back <- summary(dat, transf=transf.ipft, ni=dat$ni)
# Now the back-transformation is applied to each transformed proportion
with the study-specific sample sizes.
# The yi values are now the back-transformed values (i.e., the raw
proportions)
# and the ci.lb and ci.ub values are the back-transformed 95% CI bounds.1)
# Finally, we can create the forest plot by directly passing the observed
outcomes (i.e., proportions)
# and the CI bounds to the function. Then the back-transformed average with
the corresponding CI bounds obtained earlier
# can be added to the plot with the addpoly() function. We add a couple
tweaks to make the final forest plot look nice:
forest(dat.back$yi, ci.lb=dat.back$ci.lb, ci.ub=dat.back$ci.ub,
xlim=c(-.5,1.8), alim=c(0,1), ylim=c(-1,8), refline=NA, digits=3,
xlab="Proportion")
addpoly(pred$pred, ci.lb=pred$ci.lb, ci.ub=pred$ci.ub, row=-0.5, digits=3,
mlab="FE Model", efac=1.3)
abline(h=0.5)
text(-0.5, 7, "Study", pos=4)
text( 1.8, 7, "Proportion [95% CI]", pos=2)
### However, this is wrong
radial(res)
funnel(res)
## sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=Italian_Italy.1252 LC_CTYPE=Italian_Italy.1252
[3] LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C
[5] LC_TIME=Italian_Italy.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] metafor_1.9-3 Matrix_1.1-0 Formula_1.1-1
loaded via a namespace (and not attached):
[1] grid_3.0.2 lattice_0.20-24
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
