[R] arcsine transformation with metafor

[Antonello Preti]

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