retrodesign.py, from Gelman & Carlin (2014), Beyond Power Calculations.

retrodesign.py

import numpy as np
import scipy.stats as st
import warnings

def retrodesign(theta, se, alpha=.05 , df=1, n_sims = 10000):
    """
    Implements the "retrodesign" function from Gelman \& Carlin (2014).

    This computes the power, Sign error rate, and exaggeration ratio.

    Arguments
    ----------
    theta   :   float
                true effect size
    se      :   float
                standard error of effect estimator sampling distribution
    alpha   :   float
                error rate for the analysis
    df      :   int
                degrees of freedom in the estimate
    n_sims  :   int
                number of simulations to conduct to estimate the exaggeration ratio
    """
    if not np.isfinite(df) and df > 0:
        warn("SciPy's 't' distribution does not support infinite degrees of freedom. Setting sufficiently large")
        df = 1000000
    if df < 1:
        raise ValueError('Degrees of freedom is less than 1.')
    z = st.t.ppf(1 - alpha/2, df=df)
    upper = 1 - st.t.cdf(z - theta/se, df=df)
    lower = st.t.cdf(-z - theta/se, df=df)
    power = upper + lower
    type_S = lower / power

    simulations = theta + se*st.t.rvs(df=df, size=n_sims)
    sig_filter= np.abs(simulations) > se*z

    sig_mags = np.abs(simulations[sig_filter])
    exaggeration = sig_mags.mean() / theta

    return power, type_S, exaggeration