Args:
            full_mod: statsmodels GLM to replace MLE model. For when 'predictor'
                is not in the fixed part of the model.
            points: number of points to use for LL approximation.

def _get_slope_stats(self, exog, predictor, sd_corr, full_mod=None,
                         points=4):
        '''
        Args:
            full_mod: statsmodels GLM to replace MLE model. For when 'predictor'
                is not in the fixed part of the model.
            points: number of points to use for LL approximation.
        '''

        if full_mod is None:
            full_mod = self.mle

        # figure out which column of exog to drop for the null model
        keeps = [i for i, x in enumerate(list(exog.columns))
                 if not np.array_equal(predictor, exog[x].values.flatten())]
        i = [x for x in range(exog.shape[1]) if x not in keeps][0]

        # get log-likelihood values from beta=0 to beta=MLE
        values = np.linspace(0., full_mod.params[i], points)
        # if there are multiple predictors, use statsmodels to optimize the LL
        if keeps:
            null = [GLM(endog=self.model.y.data, exog=exog,
                           family=self.model.family.smfamily()).fit_constrained(
                                str(exog.columns[i])+'='+str(val),
                                start_params=full_mod.params.values)
                    for val in values[:-1]]
            null = np.append(null, full_mod)
            ll = np.array([x.llf for x in null])
        # if just a single predictor, use statsmodels to evaluate the LL
        else:
            null = [self.model.family.smfamily().loglike(
                np.squeeze(self.model.y.data), val*predictor)
                for val in values[:-1]]
            ll = np.append(null, full_mod.llf)

        # compute params of quartic approximatino to log-likelihood
        # c: intercept, d: shift parameter
        # a: quartic coefficient, b: quadratic coefficient

        c, d = ll[-1], -(full_mod.params[i].item())
        X = np.array([(values+d)**4,
                       (values+d)**2]).T
        a, b = np.squeeze(
            np.linalg.multi_dot([
                np.linalg.inv(np.dot(X.T, X)),
                X.T,
                (ll[:, None] - c)
            ])
        )

        # m, v: mean and variance of beta distribution of correlations
        # p, q: corresponding shape parameters of beta distribution
        m = .5
        v = sd_corr**2/4
        p = m*(m*(1-m)/v - 1)
        q = (1-m)*(m*(1-m)/v - 1)

        # function to return central moments of rescaled beta distribution
        def moment(k): return (2*p/(p+q))**k * hyp2f1(p, -k, p+q, (p+q)/p)

        # evaluate the derivatives of beta = f(correlation).
        # dict 'point' gives points about which to Taylor expand. We want to
        # expand about the mean (generally 0), but some of the derivatives
        # do not exist at 0. Evaluating at a point very close to 0 (e.g., .001)
        # generally gives good results, but the higher order the expansion, the
        # further from 0 we need to evaluate the derivatives, or they blow up.
        point = dict(zip(range(1, 14), 2**np.linspace(-1, 5, 13)/100))
        vals = dict(a=a, b=b, n=len(self.model.y.data), r=point[self.taylor])
        _deriv = [eval(x, globals(), vals) for x in self.deriv]

        # compute and return the approximate SD
        def term(i, j):
            return 1/np.math.factorial(i) * 1/np.math.factorial(j) \
                * _deriv[i] * _deriv[j] \
                * (moment(i+j) - moment(i)*moment(j))
        terms = [term(i, j) for i in range(1, self.taylor+1)
                 for j in range(1, self.taylor+1)]
        return np.array(terms).sum()**.5