Skip to content

Commit ab1a974

Browse files
author
paolinomanghi
committed
first commit
0 parents  commit ab1a974

1 file changed

Lines changed: 326 additions & 0 deletions

File tree

‎meta_analyses.py‎

Lines changed: 326 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,326 @@
1+
#!/usr/bin/env python
2+
3+
import sys, os
4+
import pandas as pd
5+
import numpy as np
6+
from scipy import stats as sts
7+
import statsmodels.api as sm
8+
import statsmodels.formula.api as smf
9+
from statsmodels.stats.multitest import fdrcorrection
10+
11+
12+
class generalized_meta_analysis(object):
13+
def __init__( self, \
14+
effects, \
15+
variances, \
16+
study_pvalues, \
17+
study_names, \
18+
n_controls, \
19+
n_cases, \
20+
response_var, \
21+
HET="PM", \
22+
overlap_mat=None, \
23+
cov="standardized_mean_diff"):
24+
25+
print("Caller of meta-analysis on %s" %response_var)
26+
print("Heterogeneity: %s" %HET)
27+
print("Number studies: %i" %len(effects))
28+
print("Correlation matrix foreseen: ", not(overlap_mat is None))
29+
30+
self.effects = np.array(effects, dtype=np.float64)
31+
self.n_cases = n_cases
32+
self.n_controls = n_controls
33+
34+
if not any([(x is None) for x in self.n_cases]):
35+
self.tot_n_cases = np.sum(self.n_cases)
36+
if not any([(x is None) for x in self.n_controls]):
37+
self.tot_n_ctrs = np.sum(self.n_controls)
38+
39+
self.study_names = study_names
40+
print("Studies: " + " ".join(self.study_names))
41+
42+
self.n = len(self.study_names)
43+
self.n_studies = self.n
44+
self.variances = np.array(variances, dtype=np.float64)
45+
print("Variances: is the sum zero? -> ", self.variances)
46+
self.HET = HET
47+
self.response_var = response_var
48+
self.devs = np.sqrt(self.variances)
49+
self.var_covar = None
50+
51+
if not any([(x is None) for x in study_pvalues]):
52+
self.study_pvalues = np.array(study_pvalues, dtype=np.float64)
53+
54+
if overlap_mat is None:
55+
self.w = np.array( [(1./v) for v in self.variances], dtype=np.float64 )
56+
self.effects_are_iid = True
57+
58+
else:
59+
self.overlap_mat = np.array(overlap_mat, dtype=np.float64)
60+
self.var_covar = np.eye( len( self.variances ) ) * self.variances
61+
print("Overlap of samples across studies in a two-by-two matrix: ", self.overlap_mat)
62+
63+
if cov == "standardized_mean_difference":
64+
for ith in range(len(self.effects)):
65+
for jth in range(len(self.effects)):
66+
if ith != jth:
67+
d_ith, d_jth = self.effects[ith], self.effects[jth]
68+
n0 = self.overlap_mat[ith, jth]
69+
if n0 > 0:
70+
N = self.n_cases[ith] + self.n_cases[jth] + ( 2 * n0 )
71+
self.var_covar[ith, jth] += ((d_ith*d_jth) / (2*( N-3 ))) + (1./n0)
72+
73+
elif cov == "linsullivan":
74+
for ith in range(len(self.effects)):
75+
for jth in range(len(self.effects)):
76+
if ith != jth:
77+
se_ij = self.devs[ith] * self.devs[jth]
78+
n0 = self.overlap_mat[ith, jth]
79+
if n0 > 0:
80+
term_a = self.n_cases[ith] * self.n_cases[jth]
81+
term_b = self.n_controls[ith] * self.n_controls[jth]
82+
N = (self.n_controls[ith] + self.n_cases[ith]) * (self.n_controls[jth] + self.n_cases[jth])
83+
self.var_covar[ith, jth] += ((n0 * np.sqrt(term_a/term_b)) / np.sqrt(N)) * (se_ij)
84+
85+
elif cov == "precomputed":
86+
for ith in range(len(self.effects)):
87+
for jth in range(len(self.effects)):
88+
if ith != jth:
89+
self.var_covar[ith, jth] += self.overlap_mat[ith, jth]
90+
91+
else:
92+
raise NotImplementedError("Cov = %s is not implemented." %cov)
93+
94+
print("Parameters following the correlated structure: ")
95+
print("Reconstructed Var-Covar matrix: ", self.var_covar)
96+
#### THESE TWO ARE FOR THE COVARIANCE MATRIX OF THE WEIGTHS
97+
self.e = np.ones(len(self.variances), dtype=np.float64)
98+
99+
#print(self.e)
100+
#print(self.var_covar)
101+
#print(np.dot(self.e, self.var_covar))
102+
#print("/", np.dot(np.dot(self.e, self.var_covar), self.e))
103+
104+
inv_cov_mat = np.linalg.inv(self.var_covar)
105+
106+
self.w = np.dot(self.e, inv_cov_mat) / np.dot(np.dot(self.e, inv_cov_mat), self.e)
107+
self.effects_are_iid = False
108+
print("Weights: ", " ".join(list(map(str, self.w))))
109+
110+
mu_bar = np.sum(a*b for a,b in zip(self.w, self.effects))/np.sum(self.w)
111+
self.Q = np.sum(a*b for a,b in zip(self.w, [(x - mu_bar)**2 for x in self.effects]))
112+
self.Qtest = 2.*(1 - sts.chi2.cdf(np.abs(self.Q), len(self.effects)-1))
113+
#### H = np.sqrt(self.Q/(self.n - 1))
114+
115+
#print("variances: ", self.variances)
116+
#print(self.Q, " Q")
117+
#print(len(variances) - 1, "len var minus one")
118+
119+
self.I2 = np.max([0., (self.Q-(len(self.variances)-1))/float(self.Q)])
120+
self.t2_PM, self.t2PM_conv = paule_mandel_tau(self.effects, self.variances)
121+
self.t2_DL = ((self.Q - self.n + 1) / self.scaling( self.w )) if (self.Q > (self.n-1)) else 0.
122+
123+
if self.effects_are_iid:
124+
if self.HET == "PM": self.W = [(1./float(v+self.t2_PM)) for v in self.variances]
125+
elif self.HET.startswith("FIX"): self.W = [(1./float(v)) for v in self.variances]
126+
else: self.W = [(1./float(v+self.t2_DL)) for v in self.variances]
127+
print("Weights: ", " ".join(list(map(str, self.W))))
128+
self.RE = np.sum(self.W*self.effects)/float(np.sum(self.W))
129+
self.RE_Var = 1./float(np.sum(self.W))
130+
131+
else:
132+
self.W = self.w
133+
self.RE = np.sum(self.W*self.effects)
134+
self.RE_Var = 0.
135+
for ith in range(len(self.effects)):
136+
for jth in range(len(self.effects)):
137+
#if ith != jth: ## EITHER WE JUST SUM THE UPPER TRIANGLE
138+
139+
if ith != jth:
140+
self.RE_Var += (self.W[ith] * self.W[jth] * self.var_covar[ith, jth]) ## WE DO NOT MULTIPLY BY TWO BECAUSE WE ARE ADDING DOUBLE THE TABLE
141+
else:
142+
self.RE_Var += (self.W[ith] * self.W[jth] * self.variances[ith])
143+
144+
print("Random/Fixed Effect model main coefficient: ", self.RE)
145+
print("First round of computation of variance led to: ", self.RE_Var, end="\n" if self.effects_are_iid else " ")
146+
147+
print("The effect variance: ", self.RE_Var)
148+
self.stdErr = np.sqrt(self.RE_Var)
149+
self.Zscore = self.RE/self.stdErr
150+
print("Meta-analysis Zeta score: ", self.Zscore)
151+
self.Pval = 2.*(1 - sts.norm.cdf(np.abs(self.Zscore)))
152+
print("Meta-analysis p value: ", self.Pval)
153+
self.conf_int = [self.RE - 1.96*self.stdErr, self.RE + 1.96*self.stdErr]
154+
print("\n*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*.\n")
155+
156+
157+
def tot_var(self, Effects, Weights):
158+
Q = np.sum(Weights * [x**2 for x in Effects]) - ((np.sum(Weights*Effects)**2)/np.sum(Weights))
159+
return Q
160+
161+
def scaling(self, W):
162+
C = np.sum(W) - (np.sum([w**2 for w in W])/float(np.sum(W)))
163+
return C
164+
165+
def tau_squared_DL(self, Q, df, C):
166+
return (Q-df)/float(C) if (Q>df) else 0.
167+
168+
def CombinedEffect(self):
169+
return np.sum(self.W*self.effects)/float(np.sum(self.W))
170+
171+
def pretty_one_feat_print(self):
172+
_,fdr = fdrcorrection(self.study_pvalues)
173+
pp = { \
174+
"effect": list(self.effects) + [self.RE],
175+
"se": list(self.devs) + [self.stdErr],
176+
"p-val": list(self.study_pvalues) + [self.Pval],
177+
"q-val": list(fdr) + [self.Pval],
178+
"n_ctrs": list(self.n_controls) + [self.tot_n_ctrs],
179+
"n_cases": list(self.n_cases) + [self.tot_n_cases],
180+
"response": self.response_var,
181+
}
182+
return pd.DataFrame(pp, index=list(self.study_names) + ["summary"])
183+
184+
def pretty_print(self):
185+
NS = {}
186+
for eff,std,P,study in zip(self.effects, self.devs, self.study_pvalues, self.study_names):
187+
NS[str(study) + "_Effect"] = eff
188+
NS[str(study) + "_Pvalue"] = P
189+
NS[str(study) + "_SE"] = std
190+
191+
NS["RE_Effect"] = self.RE
192+
NS["RE_Pvalue"] = self.Pval
193+
NS["RE_stdErr"] = self.stdErr
194+
NS["RE_conf_int"] = ";".join(list(map(str,self.conf_int)))
195+
NS["RE_Var"] = self.RE_Var
196+
NS["Zscore"] = self.Zscore
197+
198+
if self.effects_are_iid:
199+
NS["Tau2_DL"] = self.t2_DL
200+
NS["Tau2_PM"] = self.t2_PM
201+
NS["I2"] = self.I2
202+
NS["Q"] = self.Qtest
203+
NS = pd.DataFrame(NS, index=[self.response_var])
204+
return NS
205+
206+
207+
class correlation_meta_analysis(object):
208+
def __init__(self, rhos, ers, n_studies, pvals, studies, response_name, het="PM"):
209+
self.HET = het
210+
self.responseName = response_name
211+
self.studies = studies
212+
self.study_pvalues = pvals
213+
214+
self.effects = np.arctanh(np.array(rhos, dtype=np.float64)) if not ers else np.array(rhos, dtype=np.float64)
215+
216+
self.n_studies = n_studies
217+
self.n = float(len(studies))
218+
219+
if not ers:
220+
self.vi = np.array([(1./float(n-3)) for n in self.n_studies], dtype=np.float64)
221+
self.devs = np.sqrt(self.vi)
222+
else:
223+
self.devs = np.array(ers, dtype=np.float64)
224+
self.vi = self.devs**2.
225+
226+
#self.vi = np.array([(1./float(n-1)) for n in self.n_studies], dtype=np.float64)
227+
228+
self.w = np.array([(1./float(v)) for v in self.vi], dtype=np.float64)
229+
mu_bar = np.sum(a*b for a,b in zip(self.w, self.effects))/np.sum(self.w)
230+
self.Q = np.sum(a*b for a,b in zip(self.w, [(x - mu_bar)**2 for x in self.effects]))
231+
self.Qtest = 2.*(1 - sts.chi2.cdf(np.abs(self.Q), len(self.effects)-1))
232+
233+
H = np.sqrt(self.Q/(self.n - 1))
234+
self.I2 = np.max([0., (self.Q-(len(self.vi)-1))/float(self.Q)])
235+
self.t2_PM, self.t2PM_conv = paule_mandel_tau(self.effects, self.vi)
236+
self.t2_DL = ((self.Q - self.n + 1) / self.scaling( self.w )) if (self.Q > (self.n-1)) else 0.
237+
238+
if self.HET == "PM":
239+
self.W = [(1./float(v+self.t2_PM)) for v in self.vi]
240+
elif self.HET.startswith("FIX"):
241+
self.W = [(1./float(v)) for v in self.vi]
242+
else:
243+
self.W = [(1./float(v+self.t2_DL)) for v in self.vi]
244+
245+
print("Weights: ", " ".join(list(map(str, self.W))))
246+
self.RE = np.sum(self.W*self.effects)/float(np.sum(self.W))
247+
self.RE_Var = 1./float(np.sum(self.W))
248+
249+
print("Random/Fixed Effect model main coefficient: ", self.RE)
250+
print("First round of computation of variance led to: ", self.RE_Var, end="\n")
251+
252+
print("The effect variance: ", self.RE_Var)
253+
self.stdErr = np.sqrt(self.RE_Var)
254+
self.Zscore = self.RE/self.stdErr
255+
print("Meta-analysis Zeta score: ", self.Zscore)
256+
self.Pval = 2.*(1 - sts.norm.cdf(np.abs(self.Zscore)))
257+
print("Meta-analysis p value: ", self.Pval)
258+
self.conf_int = [self.RE - 1.96*self.stdErr, self.RE + 1.96*self.stdErr]
259+
260+
#self.result = self.nice_shape(True)
261+
262+
def scaling(self, W):
263+
C = np.sum(W) - (np.sum([w**2 for w in W])/float(np.sum(W)))
264+
return C
265+
266+
def tau_squared_DL(self, Q, df, C):
267+
return (Q-df)/float(C) if (Q>df) else 0.
268+
269+
def pretty_one_feat_print(self):
270+
_,fdr = fdrcorrection(self.study_pvalues)
271+
pp = { \
272+
"effect": list(self.effects) + [self.RE],
273+
"se": list(self.devs) + [self.stdErr],
274+
"p-val": list(self.study_pvalues) + [self.Pval],
275+
"q-val": list(fdr) + [self.Pval],
276+
"n_s": self.n_studies + [np.sum(self.n_studies)],
277+
"response": self.response_var,
278+
}
279+
return pd.DataFrame(pp, index=list(self.study_names) + ["summary"])
280+
281+
def pretty_print(self):
282+
NS = {}
283+
for rho,std,P,study in zip(self.effects, self.devs, self.study_pvalues, self.studies):
284+
NS[study + "_Correlation"] = np.tanh(rho)
285+
NS[study + "_Pvalue"] = P
286+
NS[study + "_SE"] = np.tanh(std)
287+
288+
NS["RE_Correlation"] = np.tanh(self.RE)
289+
NS["RE_Pvalue"] = self.Pval
290+
NS["RE_stdErr"] = np.tanh(self.stdErr)
291+
NS["RE_conf_int"] = ";".join(list(map(str, [np.tanh(c) for c in self.conf_int])))
292+
NS["RE_Var"] = np.tanh(self.RE_Var)
293+
NS["Zscore"] = self.Zscore
294+
NS["Tau2_DL"] = self.t2_DL
295+
NS["Tau2_PM"] = self.t2_PM
296+
NS["I2"] = self.I2
297+
NS["Q"] = self.Qtest
298+
NS = pd.DataFrame(NS, index=[self.responseName])
299+
return NS
300+
301+
## DISCLAIMER ##
302+
## FOLLOWING CODE FOR PAULE MANDEL TAU WAS TAKEN DIRECTLY FROM statsmodels library
303+
## I DON T OWN THIS CODE
304+
def paule_mandel_tau(eff, var_eff, tau2_start=0, atol=1e-5, maxiter=50):
305+
tau2 = tau2_start
306+
k = eff.shape[0]
307+
converged = False
308+
for i in range(maxiter):
309+
w = 1 / (var_eff + tau2)
310+
m = w.dot(eff) / w.sum(0)
311+
resid_sq = (eff - m)**2
312+
q_w = w.dot(resid_sq)
313+
# estimating equation
314+
ee = q_w - (k - 1)
315+
if ee < 0:
316+
tau2 = 0
317+
converged = 0
318+
break
319+
if np.allclose(ee, 0, atol=atol):
320+
converged = True
321+
break
322+
# update tau2
323+
delta = ee / (w**2).dot(resid_sq)
324+
tau2 += delta
325+
return tau2, converged
326+

0 commit comments

Comments
 (0)