#=== CURRENT ANALYSIS - no variance stabilizing transformation ===========================================
library(metafor)

dat <- read.table("ir.cancer.rec.csv", sep=";", header=TRUE)
dat <- escalc(measure="IR", xi=xi, ti=ti, data=dat)

# Because the incidence rate meta-analysis shows the incidence rates as proportions, 
# i multiplied the effect x100 size to make it show the numbers as % per year
dat$yi.p <-dat$yi*100
dat$vi.p <-dat$vi*10000

res <- rma(yi=yi.p, vi=vi.p, slab=dat$study, data=dat)
par(mar=c(4,2,1,1))

forest(res, xlab="Incidence rate of cancer recurrence", 
       psize=(c((dat[,2]^0.5)/10)),
       xlim=c(-10,18), alim=c(-5,10), ylim=c(-1,43), refline=0, digits=2, cex=0.8)
text(-10, 42, "Study", pos=4, cex=0.8)
text( 18, 42, "% per year [95% CI]", pos=2, cex=0.8)

#=== NEW ANALYSIS - with FT (double arcsine)==============================================================
# To perform a Freeman-Tukey (double arcsine) transformation, I followed the steps from
# http://www.metafor-project.org/doku.php/analyses:miller1978

dat <- read.table("ir.cancer.rec.csv", sep=";", header=TRUE)
dat <- escalc(measure="IRFT", xi=xi, ti=ti, data=dat)
dat[1:10,]

# Making them into percentages
dat$yi.p <-dat$yi*100
dat$vi.p <-dat$vi*10000
dat[1:10,]

# Back-transforming
transf.irft(dat$yi, dat$ti)

# The following transformation was explained for proportion meta-analysis (IPFT). 
# Is it also applicable for incidence-rates (IRFT)?

res2 <- rma(yi=dat$yi, vi=dat$vi, method="FE", data=dat)
res2
pred <- predict(res2, transf=transf.ipft.hm, targs=list(ni=dat$ni))
pred

# Making them into percentages
pred.p <- pred$pred*1000
ci.lb.p <- pred$ci.lb*1000
ci.ub.p <- pred$ci.ub*1000
pred.p
ci.lb.p
ci.ub.p

dat.back <- summary(dat, transf=transf.ipft, ni=dat$ni)
dat.back
# Again, transforming the proportion values into more understandable values (%)
dat.back$yi.p <-dat.back$yi*100
dat.back$vi.p <-dat.back$vi*10000
dat.back$ci.lb.p <-dat.back$ci.lb*100
dat.back$ci.ub.p <-dat.back$ci.ub*100
dat.back

forest(dat.back$yi.p, ci.lb=dat.back$ci.lb.p, ci.ub=dat.back$ci.ub.p,
       xlab="Incidence rate of cancer recurrence", psize=(c((dat.back[,2]^0.5)/10)),
       xlim=c(-4,15), alim=c(0,10), ylim=c(-1,43), refline=NA, digits=2, cex=0.8)
abline(h=0.5)
text(-4, 42, "Study", pos=4, cex=0.8)
text( 15, 42, "% per year [95% CI]", pos=2, cex=0.8)

# Using percentage values for the effect size estimate as well
addpoly(pred.p, ci.lb=ci.lb.p, ci.ub=ci.ub.p, row=-0.5, digits=2, mlab="RE Model", cex=0.8)