# -*- coding: utf-8 -*-
"""
Created on Thu Apr 2 14:34:25 2020
Author: Josef Perktold
License: BSD-3
"""
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.base import HolderTuple
[docs]class CombineResults(object):
"""Results from combined estimate of means or effect sizes
This currently includes intermediate results that might be removed
"""
def __init__(self, **kwds):
self.__dict__.update(kwds)
self._ini_keys = list(kwds.keys())
self.df_resid = self.k - 1
# TODO: move to property ?
self.sd_eff_w_fe_hksj = np.sqrt(self.var_hksj_fe)
self.sd_eff_w_re_hksj = np.sqrt(self.var_hksj_re)
# explained variance measures
self.h2 = self.q / (self.k - 1)
self.i2 = 1 - 1 / self.h2
# memoize ci_samples
self.cache_ci = {}
[docs] def conf_int_samples(self, alpha=0.05, use_t=None, nobs=None,
ci_func=None):
"""confidence intervals for the effect size estimate of samples
Additional information needs to be provided for confidence intervals
that are not based on normal distribution using available variance.
This is likely to change in future.
Parameters
----------
alpha : float in (0, 1)
Significance level for confidence interval. Nominal coverage is
``1 - alpha``.
use_t : None or bool
If use_t is None, then the attribute `use_t` determines whether
normal or t-distribution is used for confidence intervals.
Specifying use_t overrides the attribute.
If use_t is false, then confidence intervals are based on the
normal distribution. If it is true, then the t-distribution is
used.
nobs : None or float
Number of observations used for degrees of freedom computation.
Only used if use_t is true.
ci_func : None or callable
User provided function to compute confidence intervals.
This is not used yet and will allow using non-standard confidence
intervals.
Returns
-------
ci_eff : tuple of ndarrays
Tuple (ci_low, ci_upp) with confidence interval computed for each
sample.
Notes
-----
CombineResults currently only has information from the combine_effects
function, which does not provide details about individual samples.
"""
# this is a bit messy, we don't have enough information about
# computing conf_int already in results for other than normal
# TODO: maybe there is a better
if (alpha, use_t) in self.cache_ci:
return self.cache_ci[(alpha, use_t)]
if use_t is None:
use_t = self.use_t
if ci_func is not None:
kwds = {"use_t": use_t} if use_t is not None else {}
ci_eff = ci_func(alpha=alpha, **kwds)
self.ci_sample_distr = "ci_func"
else:
if use_t is False:
crit = stats.norm.isf(alpha / 2)
self.ci_sample_distr = "normal"
else:
if nobs is not None:
df_resid = nobs - 1
crit = stats.t.isf(alpha / 2, df_resid)
self.ci_sample_distr = "t"
else:
msg = ("`use_t=True` requires `nobs` for each sample "
"or `ci_func`. Using normal distribution for "
"confidence interval of individual samples.")
import warnings
warnings.warn(msg)
crit = stats.norm.isf(alpha / 2)
self.ci_sample_distr = "normal"
# sgn = np.asarray([-1, 1])
# ci_eff = self.eff + sgn * crit * self.sd_eff
ci_low = self.eff - crit * self.sd_eff
ci_upp = self.eff + crit * self.sd_eff
ci_eff = (ci_low, ci_upp)
# if (alpha, use_t) not in self.cache_ci: # not needed
self.cache_ci[(alpha, use_t)] = ci_eff
return ci_eff
[docs] def conf_int(self, alpha=0.05, use_t=None):
"""confidence interval for the overall mean estimate
Parameters
----------
alpha : float in (0, 1)
Significance level for confidence interval. Nominal coverage is
``1 - alpha``.
use_t : None or bool
If use_t is None, then the attribute `use_t` determines whether
normal or t-distribution is used for confidence intervals.
Specifying use_t overrides the attribute.
If use_t is false, then confidence intervals are based on the
normal distribution. If it is true, then the t-distribution is
used.
Returns
-------
ci_eff_fe : tuple of floats
Confidence interval for mean effects size based on fixed effects
model with scale=1.
ci_eff_re : tuple of floats
Confidence interval for mean effects size based on random effects
model with scale=1
ci_eff_fe_wls : tuple of floats
Confidence interval for mean effects size based on fixed effects
model with estimated scale corresponding to WLS, ie. HKSJ.
ci_eff_re_wls : tuple of floats
Confidence interval for mean effects size based on random effects
model with estimated scale corresponding to WLS, ie. HKSJ.
If random effects method is fully iterated, i.e. Paule-Mandel, then
the estimated scale is 1.
"""
if use_t is None:
use_t = self.use_t
if use_t is False:
crit = stats.norm.isf(alpha / 2)
else:
crit = stats.t.isf(alpha / 2, self.df_resid)
sgn = np.asarray([-1, 1])
m_fe = self.mean_effect_fe
m_re = self.mean_effect_re
ci_eff_fe = m_fe + sgn * crit * self.sd_eff_w_fe
ci_eff_re = m_re + sgn * crit * self.sd_eff_w_re
ci_eff_fe_wls = m_fe + sgn * crit * np.sqrt(self.var_hksj_fe)
ci_eff_re_wls = m_re + sgn * crit * np.sqrt(self.var_hksj_re)
return ci_eff_fe, ci_eff_re, ci_eff_fe_wls, ci_eff_re_wls
[docs] def test_homogeneity(self):
"""Test whether the means of all samples are the same
currently no options, test uses chisquare distribution
default might change depending on `use_t`
Returns
-------
res : HolderTuple instance
The results include the following attributes:
- statistic : float
Test statistic, ``q`` in meta-analysis, this is the
pearson_chi2 statistic for the fixed effects model.
- pvalue : float
P-value based on chisquare distribution.
- df : float
Degrees of freedom, equal to number of studies or samples
minus 1.
"""
pvalue = stats.chi2.sf(self.q, self.k - 1)
res = HolderTuple(statistic=self.q,
pvalue=pvalue,
df=self.k - 1,
distr="chi2")
return res
[docs] def summary_array(self, alpha=0.05, use_t=None):
"""Create array with sample statistics and mean estimates
Parameters
----------
alpha : float in (0, 1)
Significance level for confidence interval. Nominal coverage is
``1 - alpha``.
use_t : None or bool
If use_t is None, then the attribute `use_t` determines whether
normal or t-distribution is used for confidence intervals.
Specifying use_t overrides the attribute.
If use_t is false, then confidence intervals are based on the
normal distribution. If it is true, then the t-distribution is
used.
Returns
-------
res : ndarray
Array with columns
['eff', "sd_eff", "ci_low", "ci_upp", "w_fe","w_re"].
Rows include statistics for samples and estimates of overall mean.
column_names : list of str
The names for the columns, used when creating summary DataFrame.
"""
ci_low, ci_upp = self.conf_int_samples(alpha=alpha, use_t=use_t)
res = np.column_stack([self.eff, self.sd_eff,
ci_low, ci_upp,
self.weights_rel_fe, self.weights_rel_re])
ci = self.conf_int(alpha=alpha, use_t=use_t)
res_fe = [[self.mean_effect_fe, self.sd_eff_w_fe,
ci[0][0], ci[0][1], 1, np.nan]]
res_re = [[self.mean_effect_re, self.sd_eff_w_re,
ci[1][0], ci[1][1], np.nan, 1]]
res_fe_wls = [[self.mean_effect_fe, self.sd_eff_w_fe_hksj,
ci[2][0], ci[2][1], 1, np.nan]]
res_re_wls = [[self.mean_effect_re, self.sd_eff_w_re_hksj,
ci[3][0], ci[3][1], np.nan, 1]]
res = np.concatenate([res, res_fe, res_re, res_fe_wls, res_re_wls],
axis=0)
column_names = ['eff', "sd_eff", "ci_low", "ci_upp", "w_fe", "w_re"]
return res, column_names
[docs] def summary_frame(self, alpha=0.05, use_t=None):
"""Create DataFrame with sample statistics and mean estimates
Parameters
----------
alpha : float in (0, 1)
Significance level for confidence interval. Nominal coverage is
``1 - alpha``.
use_t : None or bool
If use_t is None, then the attribute `use_t` determines whether
normal or t-distribution is used for confidence intervals.
Specifying use_t overrides the attribute.
If use_t is false, then confidence intervals are based on the
normal distribution. If it is true, then the t-distribution is
used.
Returns
-------
res : DataFrame
pandas DataFrame instance with columns
['eff', "sd_eff", "ci_low", "ci_upp", "w_fe","w_re"].
Rows include statistics for samples and estimates of overall mean.
"""
if use_t is None:
use_t = self.use_t
labels = (list(self.row_names) +
["fixed effect", "random effect",
"fixed effect wls", "random effect wls"])
res, col_names = self.summary_array(alpha=alpha, use_t=use_t)
results = pd.DataFrame(res, index=labels, columns=col_names)
return results
[docs] def plot_forest(self, ax=None, **kwds):
"""Forest plot with means and confidence intervals
Parameters
----------
ax : None or matplotlib axis instance
If ax is provided, then the plot will be added to it.
kwds : optional keyword arguments
Keywords are forwarded to the dot_plot function that creates the
plot.
Returns
-------
fig : Matplotlib figure instance
See Also
--------
dot_plot
"""
from statsmodels.graphics.dotplots import dot_plot
res_df = self.summary_frame()
hw = np.abs(res_df[["ci_low", "ci_upp"]] - res_df[["eff"]].values)
fig = dot_plot(points=res_df["eff"], intervals=hw,
lines=res_df.index, line_order=res_df.index, **kwds)
return fig
[docs]def effectsize_smd(mean1, sd1, nobs1, mean2, sd2, nobs2):
"""effect sizes for mean difference for use in meta-analysis
mean1, sd1, nobs1 are for treatment
mean2, sd2, nobs2 are for control
Effect sizes are computed for the mean difference ``mean1 - mean2``
standardized by an estimate of the within variance.
This does not have option yet.
It uses standardized mean difference with bias correction as effect size.
This currently does not use np.asarray, all computations are possible in
pandas.
Parameters
----------
mean1 : array
mean of second sample, treatment groups
sd1 : array
standard deviation of residuals in treatment groups, within
nobs1 : array
number of observations in treatment groups
mean2, sd2, nobs2 : arrays
mean, standard deviation and number of observations of control groups
Returns
-------
smd_bc : array
bias corrected estimate of standardized mean difference
var_smdbc : array
estimate of variance of smd_bc
Notes
-----
Status: API will still change. This is currently intended for support of
meta-analysis.
References
----------
Borenstein, Michael. 2009. Introduction to Meta-Analysis.
Chichester: Wiley.
Chen, Ding-Geng, and Karl E. Peace. 2013. Applied Meta-Analysis with R.
Chapman & Hall/CRC Biostatistics Series.
Boca Raton: CRC Press/Taylor & Francis Group.
"""
# TODO: not used yet, design and options ?
# k = len(mean1)
# if row_names is None:
# row_names = list(range(k))
# crit = stats.norm.isf(alpha / 2)
# var_diff_uneq = sd1**2 / nobs1 + sd2**2 / nobs2
var_diff = (sd1**2 * (nobs1 - 1) +
sd2**2 * (nobs2 - 1)) / (nobs1 + nobs2 - 2)
sd_diff = np.sqrt(var_diff)
nobs = nobs1 + nobs2
bias_correction = 1 - 3 / (4 * nobs - 9)
smd = (mean1 - mean2) / sd_diff
smd_bc = bias_correction * smd
var_smdbc = nobs / nobs1 / nobs2 + smd_bc**2 / 2 / (nobs - 3.94)
return smd_bc, var_smdbc
[docs]def effectsize_2proportions(count1, nobs1, count2, nobs2, statistic="diff",
zero_correction=None, zero_kwds=None):
"""Effects sizes for two sample binomial proportions
Parameters
----------
count1, nobs1, count2, nobs2 : array_like
data for two samples
statistic : {"diff", "odds-ratio", "risk-ratio", "arcsine"}
statistic for the comparison of two proportions
Effect sizes for "odds-ratio" and "risk-ratio" are in logarithm.
zero_correction : {None, float, "tac", "clip"}
Some statistics are not finite when zero counts are in the data.
The options to remove zeros are:
* float : if zero_correction is a single float, then it will be added
to all count (cells) if the sample has any zeros.
* "tac" : treatment arm continuity correction see Ruecker et al 2009,
section 3.2
* "clip" : clip proportions without adding a value to all cells
The clip bounds can be set with zero_kwds["clip_bounds"]
zero_kwds : dict
additional options to handle zero counts
"clip_bounds" tuple, default (1e-6, 1 - 1e-6) if zero_correction="clip"
other options not yet implemented
Returns
-------
effect size : array
Effect size for each sample.
var_es : array
Estimate of variance of the effect size
Notes
-----
Status: API is experimental, Options for zero handling is incomplete.
The names for ``statistics`` keyword can be shortened to "rd", "rr", "or"
and "as".
The statistics are defined as:
- risk difference = p1 - p2
- log risk ratio = log(p1 / p2)
- log odds_ratio = log(p1 / (1 - p1) * (1 - p2) / p2)
- arcsine-sqrt = arcsin(sqrt(p1)) - arcsin(sqrt(p2))
where p1 and p2 are the estimated proportions in sample 1 (treatment) and
sample 2 (control).
log-odds-ratio and log-risk-ratio can be transformed back to ``or`` and
`rr` using `exp` function.
See Also
--------
statsmodels.stats.contingency_tables
"""
if zero_correction is None:
cc1 = cc2 = 0
elif zero_correction == "tac":
# treatment arm continuity correction Ruecker et al 2009, section 3.2
nobs_t = nobs1 + nobs2
cc1 = nobs2 / nobs_t
cc2 = nobs1 / nobs_t
elif zero_correction == "clip":
clip_bounds = zero_kwds.get("clip_bounds", (1e-6, 1 - 1e-6))
cc1 = cc2 = 0
elif zero_correction:
# TODO: check is float_like
cc1 = cc2 = zero_correction
else:
msg = "zero_correction not recognized or supported"
raise NotImplementedError(msg)
zero_mask1 = (count1 == 0) | (count1 == nobs1)
zero_mask2 = (count2 == 0) | (count2 == nobs2)
zmask = np.logical_or(zero_mask1, zero_mask2)
n1 = nobs1 + (cc1 + cc2) * zmask
n2 = nobs2 + (cc1 + cc2) * zmask
p1 = (count1 + cc1) / (n1)
p2 = (count2 + cc2) / (n2)
if zero_correction == "clip":
p1 = np.clip(p1, *clip_bounds)
p2 = np.clip(p2, *clip_bounds)
if statistic in ["diff", "rd"]:
rd = p1 - p2
rd_var = p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2
eff = rd
var_eff = rd_var
elif statistic in ["risk-ratio", "rr"]:
# rr = p1 / p2
log_rr = np.log(p1) - np.log(p2)
log_rr_var = (1 - p1) / p1 / n1 + (1 - p2) / p2 / n2
eff = log_rr
var_eff = log_rr_var
elif statistic in ["odds-ratio", "or"]:
# or_ = p1 / (1 - p1) * (1 - p2) / p2
log_or = np.log(p1) - np.log(1 - p1) - np.log(p2) + np.log(1 - p2)
log_or_var = 1 / (p1 * (1 - p1) * n1) + 1 / (p2 * (1 - p2) * n2)
eff = log_or
var_eff = log_or_var
elif statistic in ["arcsine", "arcsin", "as"]:
as_ = np.arcsin(np.sqrt(p1)) - np.arcsin(np.sqrt(p2))
as_var = (1 / n1 + 1 / n2) / 4
eff = as_
var_eff = as_var
else:
msg = 'statistic not recognized, use one of "rd", "rr", "or", "as"'
raise NotImplementedError(msg)
return eff, var_eff
[docs]def combine_effects(effect, variance, method_re="iterated", row_names=None,
use_t=False, alpha=0.05, **kwds):
"""combining effect sizes for effect sizes using meta-analysis
This currently does not use np.asarray, all computations are possible in
pandas.
Parameters
----------
effect : array
mean of effect size measure for all samples
variance : array
variance of mean or effect size measure for all samples
method_re : {"iterated", "chi2"}
method that is use to compute the between random effects variance
"iterated" or "pm" uses Paule and Mandel method to iteratively
estimate the random effects variance. Options for the iteration can
be provided in the ``kwds``
"chi2" or "dl" uses DerSimonian and Laird one-step estimator.
row_names : list of strings (optional)
names for samples or studies, will be included in results summary and
table.
alpha : float in (0, 1)
significance level, default is 0.05, for the confidence intervals
Returns
-------
results : CombineResults
Contains estimation results and intermediate statistics, and includes
a method to return a summary table.
Statistics from intermediate calculations might be removed at a later
time.
Notes
-----
Status: Basic functionality is verified, mainly compared to R metafor
package. However, API might still change.
This computes both fixed effects and random effects estimates. The
random effects results depend on the method to estimate the RE variance.
Scale estimate
In fixed effects models and in random effects models without fully
iterated random effects variance, the model will in general not account
for all residual variance. Traditional meta-analysis uses a fixed
scale equal to 1, that might not produce test statistics and
confidence intervals with the correct size. Estimating the scale to account
for residual variance often improves the small sample properties of
inference and confidence intervals.
This adjustment to the standard errors is often referred to as HKSJ
method based attributed to Hartung and Knapp and Sidik and Jonkman.
However, this is equivalent to estimating the scale in WLS.
The results instance includes both, fixed scale and estimated scale
versions of standard errors and confidence intervals.
References
----------
Borenstein, Michael. 2009. Introduction to Meta-Analysis.
Chichester: Wiley.
Chen, Ding-Geng, and Karl E. Peace. 2013. Applied Meta-Analysis with R.
Chapman & Hall/CRC Biostatistics Series.
Boca Raton: CRC Press/Taylor & Francis Group.
"""
k = len(effect)
if row_names is None:
row_names = list(range(k))
crit = stats.norm.isf(alpha / 2)
# alias for initial version
eff = effect
var_eff = variance
sd_eff = np.sqrt(var_eff)
# fixed effects computation
weights_fe = 1 / var_eff # no bias correction ?
w_total_fe = weights_fe.sum(0)
weights_rel_fe = weights_fe / w_total_fe
eff_w_fe = weights_rel_fe * eff
mean_effect_fe = eff_w_fe.sum()
var_eff_w_fe = 1 / w_total_fe
sd_eff_w_fe = np.sqrt(var_eff_w_fe)
# random effects computation
q = (weights_fe * eff**2).sum(0)
q -= (weights_fe * eff).sum()**2 / w_total_fe
df = k - 1
if method_re.lower() in ["iterated", "pm"]:
tau2, _ = _fit_tau_iterative(eff, var_eff, **kwds)
elif method_re.lower() in ["chi2", "dl"]:
c = w_total_fe - (weights_fe**2).sum() / w_total_fe
tau2 = (q - df) / c
else:
raise ValueError('method_re should be "iterated" or "chi2"')
weights_re = 1 / (var_eff + tau2) # no bias_correction ?
w_total_re = weights_re.sum(0)
weights_rel_re = weights_re / weights_re.sum(0)
eff_w_re = weights_rel_re * eff
mean_effect_re = eff_w_re.sum()
var_eff_w_re = 1 / w_total_re
sd_eff_w_re = np.sqrt(var_eff_w_re)
# ci_low_eff_re = mean_effect_re - crit * sd_eff_w_re
# ci_upp_eff_re = mean_effect_re + crit * sd_eff_w_re
scale_hksj_re = (weights_re * (eff - mean_effect_re)**2).sum() / df
scale_hksj_fe = (weights_fe * (eff - mean_effect_fe)**2).sum() / df
var_hksj_re = (weights_rel_re * (eff - mean_effect_re)**2).sum() / df
var_hksj_fe = (weights_rel_fe * (eff - mean_effect_fe)**2).sum() / df
res = CombineResults(**locals())
return res
[docs]def _fit_tau_iterative(eff, var_eff, tau2_start=0, atol=1e-5, maxiter=50):
"""Paule-Mandel iterative estimate of between random effect variance
implementation follows DerSimonian and Kacker 2007 Appendix 8
see also Kacker 2004
Parameters
----------
eff : ndarray
effect sizes
var_eff : ndarray
variance of effect sizes
tau2_start : float
starting value for iteration
atol : float, default: 1e-5
convergence tolerance for absolute value of estimating equation
maxiter : int
maximum number of iterations
Returns
-------
tau2 : float
estimate of random effects variance tau squared
converged : bool
True if iteration has converged.
"""
tau2 = tau2_start
k = eff.shape[0]
converged = False
for i in range(maxiter):
w = 1 / (var_eff + tau2)
m = w.dot(eff) / w.sum(0)
resid_sq = (eff - m)**2
q_w = w.dot(resid_sq)
# estimating equation
ee = q_w - (k - 1)
if ee < 0:
tau2 = 0
converged = 0
break
if np.allclose(ee, 0, atol=atol):
converged = True
break
# update tau2
delta = ee / (w**2).dot(resid_sq)
tau2 += delta
return tau2, converged
[docs]def _fit_tau_mm(eff, var_eff, weights):
"""one-step method of moment estimate of between random effect variance
implementation follows Kacker 2004 and DerSimonian and Kacker 2007 eq. 6
Parameters
----------
eff : ndarray
effect sizes
var_eff : ndarray
variance of effect sizes
weights : ndarray
weights for estimating overall weighted mean
Returns
-------
tau2 : float
estimate of random effects variance tau squared
"""
w = weights
m = w.dot(eff) / w.sum(0)
resid_sq = (eff - m)**2
q_w = w.dot(resid_sq)
w_t = w.sum()
expect = w.dot(var_eff) - (w**2).dot(var_eff) / w_t
denom = w_t - (w**2).sum() / w_t
# moment estimate from estimating equation
tau2 = (q_w - expect) / denom
return tau2
[docs]def _fit_tau_iter_mm(eff, var_eff, tau2_start=0, atol=1e-5, maxiter=50):
"""iterated method of moment estimate of between random effect variance
This repeatedly estimates tau, updating weights in each iteration
see two-step estimators in DerSimonian and Kacker 2007
Parameters
----------
eff : ndarray
effect sizes
var_eff : ndarray
variance of effect sizes
tau2_start : float
starting value for iteration
atol : float, default: 1e-5
convergence tolerance for change in tau2 estimate between iterations
maxiter : int
maximum number of iterations
Returns
-------
tau2 : float
estimate of random effects variance tau squared
converged : bool
True if iteration has converged.
"""
tau2 = tau2_start
converged = False
for _ in range(maxiter):
w = 1 / (var_eff + tau2)
tau2_new = _fit_tau_mm(eff, var_eff, w)
tau2_new = max(0, tau2_new)
delta = tau2_new - tau2
if np.allclose(delta, 0, atol=atol):
converged = True
break
tau2 = tau2_new
return tau2, converged