> #================== Deek's Funnel Plot ============================#
> library(metafor)
> 
> ## read data
> df <- read.csv("chang_cbcl.csv")
> 
> ## calculate effect size as lnDOR
> data1 <- escalc(measure="OR", ai=TP, bi=FP, ci=FN, di=TN, data=df)
> data1
                   Study Year  TP FP  FN  TN  n1  n2      X1.4n1      X1.4n2       added       ESS     yi     vi 
1           Derks_F 2006 2006  20 36   5 155  25 191 0.010000000 0.001308901 0.011308901  88.42593 2.8462 0.2842 
2           Derks_M 2006 2006  20 14   7 151  27 165 0.009259259 0.001515152 0.010774411  92.81250 3.4280 0.2709 
3           Edwards 2015 2015  38 17  11  29  49  46 0.005102041 0.005434783 0.010536823  94.90526 1.7738 0.2105 
4           Eiraldi 2000 2000  88  6  85  30 173  36 0.001445087 0.006944444 0.008389531 119.19617 1.6441 0.2231 
5            Elkins 2014 2014  17  2   6  21  23  23 0.010869565 0.010869565 0.021739130  46.00000 3.3928 0.7731 
6            Gjevik 2015 2015  16 19   1  19  17  38 0.014705882 0.006578947 0.021284830  46.98182 2.7726 1.1678 
7             Gould 1993 1993  24 33  17 234  41 267 0.006097561 0.000936330 0.007033891 142.16883 2.3037 0.1351 
8           Hudziak 2004 2004  55 14  11 103  66 117 0.003787879 0.002136752 0.005924631 168.78689 3.6051 0.1902 
9               Kim 2005 2005  24  8   9   5  33  13 0.007575758 0.019230769 0.026806527  37.30435 0.5108 0.4778 
10          Lampert 2004 2004 341 92 114 216 455 308 0.000549451 0.000811688 0.001361139 734.67890 1.9492 0.0272 
11           Rishel 2005 2005  55 30  22 129  77 159 0.003246753 0.001572327 0.004819080 207.50847 2.3749 0.1047 
12 Roessner_germany 2007 2007 240 36   8  35 248  71 0.001008065 0.003521127 0.004529191 220.78997 3.3730 0.1855 
13  Roessner_brazil 2007 2007 131 63  23  72 154 135 0.001623377 0.001851852 0.003475228 287.75087 1.8732 0.0809 
14         Schwarte 2005 2005  39 11  22  13  61  24 0.004098361 0.010416667 0.014515027  68.89412 0.7396 0.2389 
15            Tripp 2006 2006  69 39  39  37 108  76 0.002314815 0.003289474 0.005604288 178.43478 0.5179 0.0928 
16       Wassenberg 2004 2004  16  8   3  43  19  51 0.013157895 0.004901961 0.018059856  55.37143 3.3557 0.5441 
> 
> ## make funnel plot, lnDOR against 1/root(EES)
> fplot <- funnel(data1$yi, data1$vi, ni = data1$added, yaxis="sqrtni", atransf=exp,
+                 at=log(c(1, 10, 100, 1000)))
> ## save as a new data
> newf <- fplot
> 
> ## transform lnDOR back to DOR by exponentiation
> newf$exp_x <- exp(newf$x)
> 
> ## test asymmetry 
> ## regression of lnDOR with 1/root(ESS) weighted by ESS
> newf_lm <- lm(y ~ x, weights = data1$ESS, data = newf)
> summary(newf_lm)

Call:
lm(formula = y ~ x, data = newf, weights = data1$ESS)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-0.91976 -0.03411  0.24473  0.43088  0.61162 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 0.061018   0.023038   2.649   0.0191 *
x           0.005032   0.009747   0.516   0.6137  
---
Signif. codes:  0 ．***・ 0.001 ．**・ 0.01 ．*・ 0.05 ．.・ 0.1 ． ・ 1

Residual standard error: 0.421 on 14 degrees of freedom
Multiple R-squared:  0.01868,	Adjusted R-squared:  -0.05141 
F-statistic: 0.2665 on 1 and 14 DF,  p-value: 0.6137

> 
> ## plot a new funnel plot using ggplot
> ggplot(data = newf,
+        mapping = aes(x = newf$exp_x, y = newf$y)) +
+     geom_point(color = "blue",
+                size = 2) +
+     geom_text(aes(label = c(1:nrow(df))),
+               hjust = -0.6) +
+     labs(title = "Deek's Funnel Plot Asymmetry Test",
+          x = "Diagnostic Odds Ratio",
+          y = "1/root(ESS)") +
+     scale_y_reverse(breaks = (seq(0, 20, 5) / 100), minor_breaks = NULL) +
+     scale_x_log10(breaks = c(1, 10, 100, 1000), minor_breaks = NULL) +
+     expand_limits(x = c(1, 1000), y = c(0.2, 0)) +
+     theme(panel.grid.major.x = element_blank(),
+           plot.title = element_text(hjust = 0.5))