%matplotlib inline

import numpy as np
import pandas as pd
from scipy import stats, optimize

from statsmodels.regression.linear_model import WLS
from statsmodels.genmod.generalized_linear_model import GLM

from statsmodels.stats.meta_analysis import (
    effectsize_smd, effectsize_2proportions, combine_effects,
    _fit_tau_iterative, _fit_tau_mm, _fit_tau_iter_mm)

# increase line length for pandas
pd.set_option('display.width', 100)

"""
Explanation: Meta-Analysis in statsmodels
Statsmodels include basic methods for meta-analysis. This notebook illustrates the current usage.
Status: The results have been verified against R meta and metafor packages. However, the API is still experimental and will still change. Some options for additional methods that are available in R meta and metafor are missing.
The support for meta-analysis has 3 parts:

effect size functions: this currently includes
  effectsize_smd computes effect size and their standard errors for standardized mean difference,
effectsize_2proportions computes effect sizes for comparing two independent proportions using risk difference, (log) risk ratio, (log) odds-ratio or arcsine square root transformation
The combine_effects computes fixed and random effects estimate for the overall mean or effect. The returned results instance includes a forest plot function.
helper functions to estimate the random effect variance, tau-squared

The estimate of the overall effect size in combine_effects can also be performed using WLS or GLM with var_weights.
Finally, the meta-analysis functions currently do not include the Mantel-Hanszel method. However, the fixed effects results can be computed directly using StratifiedTable as illustrated below.
End of explanation
"""


data = [
["Carroll", 94, 22,60,92, 20,60],
["Grant", 98, 21,65, 92,22, 65],
["Peck", 98, 28, 40,88 ,26, 40],
["Donat", 94,19, 200, 82,17, 200],
["Stewart", 98, 21,50, 88,22 , 45],
["Young", 96,21,85, 92 ,22, 85]]
colnames = ["study","mean_t","sd_t","n_t","mean_c","sd_c","n_c"]
rownames = [i[0] for i in data]
dframe1 = pd.DataFrame(data, columns=colnames)
rownames

mean2, sd2, nobs2, mean1, sd1, nobs1 = np.asarray(dframe1[["mean_t","sd_t","n_t","mean_c","sd_c","n_c"]]).T
rownames = dframe1["study"]
rownames.tolist()

np.array(nobs1 + nobs2)

"""
Explanation: Example
End of explanation
"""


eff, var_eff = effectsize_smd(mean2, sd2, nobs2, mean1, sd1, nobs1)

"""
Explanation: estimate effect size standardized mean difference
End of explanation
"""


res3 = combine_effects(eff, var_eff, method_re="chi2", use_t=True, row_names=rownames)
# TODO: we still need better information about conf_int of individual samples
# We don't have enough information in the model for individual confidence intervals
# if those are not based on normal distribution.
res3.conf_int_samples(nobs=np.array(nobs1 + nobs2))
print(res3.summary_frame())

res3.cache_ci

res3.method_re

fig = res3.plot_forest()
fig.set_figheight(6)
fig.set_figwidth(6)

res3 = combine_effects(eff, var_eff, method_re="chi2", use_t=False, row_names=rownames)
# TODO: we still need better information about conf_int of individual samples
# We don't have enough information in the model for individual confidence intervals
# if those are not based on normal distribution.
res3.conf_int_samples(nobs=np.array(nobs1 + nobs2))
print(res3.summary_frame())

"""
Explanation: Using one-step chi2, DerSimonian-Laird estimate for random effects variance tau
Method option for random effect method_re="chi2" or method_re="dl", both names are accepted.
This is commonly referred to as the DerSimonian-Laird method, it is based on a moment estimator based on pearson chi2 from the fixed effects estimate.
End of explanation
"""


res4 = combine_effects(eff, var_eff, method_re="iterated", use_t=False, row_names=rownames)
res4_df = res4.summary_frame()
print("method RE:", res4.method_re)
print(res4.summary_frame())
fig = res4.plot_forest()


"""
Explanation: Using iterated, Paule-Mandel estimate for random effects variance tau
The method commonly referred to as Paule-Mandel estimate is a method of moment estimate for the random effects variance that iterates between mean and variance estimate until convergence.
End of explanation
"""


eff = np.array([61.00, 61.40, 62.21, 62.30, 62.34, 62.60, 62.70,
                62.84, 65.90])
var_eff = np.array([0.2025, 1.2100, 0.0900, 0.2025, 0.3844, 0.5625,
                    0.0676, 0.0225, 1.8225])
rownames = ['PTB', 'NMi', 'NIMC', 'KRISS', 'LGC', 'NRC', 'IRMM', 'NIST', 'LNE']

res2_DL = combine_effects(eff, var_eff, method_re="dl", use_t=True, row_names=rownames)
print("method RE:", res2_DL.method_re)
print(res2_DL.summary_frame())
fig = res2_DL.plot_forest()
fig.set_figheight(6)
fig.set_figwidth(6)

res2_PM = combine_effects(eff, var_eff, method_re="pm", use_t=True, row_names=rownames)
print("method RE:", res2_PM.method_re)
print(res2_PM.summary_frame())
fig = res2_PM.plot_forest()
fig.set_figheight(6)
fig.set_figwidth(6)

"""
Explanation: Example Kacker interlaboratory mean
In this example the effect size is the mean of measurements in a lab. We combine the estimates from several labs to estimate and overall average.
End of explanation
"""


import io


ss = """\
    study,nei,nci,e1i,c1i,e2i,c2i,e3i,c3i,e4i,c4i
    1,19,22,16.0,20.0,11,12,4.0,8.0,4,3
    2,34,35,22.0,22.0,18,12,15.0,8.0,15,6
    3,72,68,44.0,40.0,21,15,10.0,3.0,3,0
    4,22,20,19.0,12.0,14,5,5.0,4.0,2,3
    5,70,32,62.0,27.0,42,13,26.0,6.0,15,5
    6,183,94,130.0,65.0,80,33,47.0,14.0,30,11
    7,26,50,24.0,30.0,13,18,5.0,10.0,3,9
    8,61,55,51.0,44.0,37,30,19.0,19.0,11,15
    9,36,25,30.0,17.0,23,12,13.0,4.0,10,4
    10,45,35,43.0,35.0,19,14,8.0,4.0,6,0
    11,246,208,169.0,139.0,106,76,67.0,42.0,51,35
    12,386,141,279.0,97.0,170,46,97.0,21.0,73,8
    13,59,32,56.0,30.0,34,17,21.0,9.0,20,7
    14,45,15,42.0,10.0,18,3,9.0,1.0,9,1
    15,14,18,14.0,18.0,13,14,12.0,13.0,9,12
    16,26,19,21.0,15.0,12,10,6.0,4.0,5,1
    17,74,75,,,42,40,,,23,30"""
df3 = pd.read_csv(io.StringIO(ss))
df_12y = df3[["e2i", "nei", "c2i", "nci"]]
# TODO: currently 1 is reference, switch labels
count1, nobs1, count2, nobs2 = df_12y.values.T
dta = df_12y.values.T

eff, var_eff = effectsize_2proportions(*dta, statistic="rd")

eff, var_eff

res5 = combine_effects(eff, var_eff, method_re="iterated", use_t=False)#, row_names=rownames)
res5_df = res5.summary_frame()
print("method RE:", res5.method_re)
print("RE variance tau2:", res5.tau2)
print(res5.summary_frame())
fig = res5.plot_forest()
fig.set_figheight(8)
fig.set_figwidth(6)

"""
Explanation: Meta-analysis of proportions
In the following example the random effect variance tau is estimated to be zero. 
I then change two counts in the data, so the second example has random effects variance greater than zero.
End of explanation
"""


dta_c = dta.copy()
dta_c.T[0, 0] = 18
dta_c.T[1, 0] = 22
dta_c.T

eff, var_eff = effectsize_2proportions(*dta_c, statistic="rd")
res5 = combine_effects(eff, var_eff, method_re="iterated", use_t=False)#, row_names=rownames)
res5_df = res5.summary_frame()
print("method RE:", res5.method_re)
print(res5.summary_frame())
fig = res5.plot_forest()
fig.set_figheight(8)
fig.set_figwidth(6)

res5 = combine_effects(eff, var_eff, method_re="chi2", use_t=False)
res5_df = res5.summary_frame()
print("method RE:", res5.method_re)
print(res5.summary_frame())
fig = res5.plot_forest()
fig.set_figheight(8)
fig.set_figwidth(6)

"""
Explanation: changing data to have positive random effects variance
End of explanation
"""


from statsmodels.genmod.generalized_linear_model import GLM

eff, var_eff = effectsize_2proportions(*dta_c, statistic="or")
res = combine_effects(eff, var_eff, method_re="chi2", use_t=False)
res_frame = res.summary_frame()
print(res_frame.iloc[-4:])

"""
Explanation: Replicate fixed effect analysis using GLM with var_weights
combine_effects computes weighted average estimates which can be replicated using GLM with var_weights or with WLS.
The scale option in GLM.fit can be used to replicate fixed meta-analysis with fixed and with HKSJ/WLS scale
End of explanation
"""


weights = 1 / var_eff
mod_glm = GLM(eff, np.ones(len(eff)),
              var_weights=weights)
res_glm = mod_glm.fit(scale=1.)
print(res_glm.summary().tables[1])

# check results
res_glm.scale, res_glm.conf_int() - res_frame.loc["fixed effect", ["ci_low", "ci_upp"]].values

"""
Explanation: We need to fix scale=1 in order to replicate standard errors for the usual meta-analysis.
End of explanation
"""


res_glm = mod_glm.fit(scale="x2")
print(res_glm.summary().tables[1])

# check results
res_glm.scale, res_glm.conf_int() - res_frame.loc["fixed effect", ["ci_low", "ci_upp"]].values

"""
Explanation: Using HKSJ variance adjustment in meta-analysis is equivalent to estimating the scale using pearson chi2, which is also the default for the gaussian family.
End of explanation
"""


t, nt, c, nc = dta_c
counts = np.column_stack([t, nt - t, c, nc - c])
ctables = counts.T.reshape(2, 2, -1)
ctables[:, :, 0]

counts[0]

dta_c.T[0]

import statsmodels.stats.api as smstats

st = smstats.StratifiedTable(ctables.astype(np.float64))

"""
Explanation: Mantel-Hanszel odds-ratio using contingency tables
The fixed effect for the log-odds-ratio using the Mantel-Hanszel can be directly computed using StratifiedTable.
We need to create a 2 x 2 x k contingency table to be used with StratifiedTable.
End of explanation
"""


st.logodds_pooled, st.logodds_pooled - 0.4428186730553189  # R meta

st.logodds_pooled_se, st.logodds_pooled_se - 0.08928560091027186  # R meta

st.logodds_pooled_confint()

print(st.test_equal_odds())

print(st.test_null_odds())

"""
Explanation: compare pooled log-odds-ratio and standard error to R meta package
End of explanation
"""


ctables.sum(1)

nt, nc

"""
Explanation: check conversion to stratified contingency table
Row sums of each table are the sample sizes for treatment and control experiments
End of explanation
"""


print(st.summary())

"""
Explanation: Results from R meta package
```

res_mb_hk = metabin(e2i, nei, c2i, nci, data=dat2, sm="OR", Q.Cochrane=FALSE, method="MH", method.tau="DL", hakn=FALSE, backtransf=FALSE)
res_mb_hk
     logOR            95%-CI %W(fixed) %W(random)
1   2.7081 [ 0.5265; 4.8896]       0.3        0.7
2   1.2567 [ 0.2658; 2.2476]       2.1        3.2
3   0.3749 [-0.3911; 1.1410]       5.4        5.4
4   1.6582 [ 0.3245; 2.9920]       0.9        1.8
5   0.7850 [-0.0673; 1.6372]       3.5        4.4
6   0.3617 [-0.1528; 0.8762]      12.1       11.8
7   0.5754 [-0.3861; 1.5368]       3.0        3.4
8   0.2505 [-0.4881; 0.9892]       6.1        5.8
9   0.6506 [-0.3877; 1.6889]       2.5        3.0
10  0.0918 [-0.8067; 0.9903]       4.5        3.9
11  0.2739 [-0.1047; 0.6525]      23.1       21.4
12  0.4858 [ 0.0804; 0.8911]      18.6       18.8
13  0.1823 [-0.6830; 1.0476]       4.6        4.2
14  0.9808 [-0.4178; 2.3795]       1.3        1.6
15  1.3122 [-1.0055; 3.6299]       0.4        0.6
16 -0.2595 [-1.4450; 0.9260]       3.1        2.3
17  0.1384 [-0.5076; 0.7844]       8.5        7.6

Number of studies combined: k = 17
                  logOR           95%-CI    z  p-value

Fixed effect model   0.4428 [0.2678; 0.6178] 4.96 < 0.0001
Random effects model 0.4295 [0.2504; 0.6086] 4.70 < 0.0001
Quantifying heterogeneity:
 tau^2 = 0.0017 [0.0000; 0.4589]; tau = 0.0410 [0.0000; 0.6774];
 I^2 = 1.1% [0.0%; 51.6%]; H = 1.01 [1.00; 1.44]
Test of heterogeneity:
     Q d.f. p-value
 16.18   16  0.4404
Details on meta-analytical method:
- Mantel-Haenszel method
- DerSimonian-Laird estimator for tau^2
- Jackson method for confidence interval of tau^2 and tau

res_mb_hk$TE.fixed
[1] 0.4428186730553189
res_mb_hk$seTE.fixed
[1] 0.08928560091027186
c(res_mb_hk$lower.fixed, res_mb_hk$upper.fixed)
[1] 0.2678221109331694 0.6178152351774684

```
End of explanation
"""


