'test_diff(groups, rho=0)
Test for difference between survival curves
Parameters
groups: A list of the values for exog to test for difference.
tests the null hypothesis that the survival curves for all
values of exog in groups are equal
rho: compute the test statistic with weight S(t)^rho, where
S(t) is the pooled estimate for the Kaplan-Meier survival function.
If rho = 0, this is the logrank test, if rho = 0, this is the
Peto and Peto modification to the Gehan-Wilcoxon test.
weight: User specified function that accepts as its sole arguement
an array of times, and returns an array of weights for each time
to be used in the test
Returns
An array whose zeroth element is the chi-square test statistic for
the global null hypothesis, that all survival curves are equal,
the index one element is degrees of freedom for the test, and the
index two element is the p-value for the test.
Examples
>>> import statsmodels.api as sm
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from statsmodels.sandbox.survival2 import KaplanMeier
>>> dta = sm.datasets.strikes.load()
>>> dta = dta.values()[-1]
>>> censoring = np.ones_like(dta[:,0])
>>> censoring[dta[:,0] > 80] = 0
>>> dta = np.c_[dta,censoring]
>>> km = KaplanMeier(dta,0,exog=1,censoring=2)
>>> km.fit()
Test for difference of survival curves
>>> log_rank = km3.test_diff([0.0645,-0.03957])
The zeroth element of log_rank is the chi-square test statistic
for the difference between the survival curves using the log rank test
for exog = 0.0645 and exog = -0.03957, the index one element
is the degrees of freedom for the test, and the index two element
is the p-value for the test
>>> wilcoxon = km.test_diff([0.0645,-0.03957], rho=1)
wilcoxon is the equivalent information as log_rank, but for the
Peto and Peto modification to the Gehan-Wilcoxon test.
User specified weight functions
>>> log_rank = km3.test_diff([0.0645,-0.03957], weight=np.ones_like)
This is equivalent to the log rank test
More than two groups
>>> log_rank = km.test_diff([0.0645,-0.03957,0.01138])
The test can be performed with arbitrarily many groups, so long as
they are all in the column exog'

def test_diff(self, groups, rho=None, weight=None):

    groups = np.asarray(groups)
    if (self.exog == None):
        raise ValueError('Need    an    exogenous    variable    for    logrank    test')
    elif np.in1d(groups, self.groups).all():
        data = self.data[np.in1d(self.data[:, self.exog], groups)]
        t = data[:, self.endog].astype(float).astype(int)
        tind = np.unique(t)
        NK = []
        N = []
        D = []
        Z = []
        if ((rho != None) and (weight != None)):
            raise ValueError('Must    use    either    rho    or    weights,    not    both')
        elif (rho != None):
            s = KaplanMeier(data, self.endog, censoring=self.censoring)
            s.fit()
            s = (s.results[0][0] ** rho)
            s = np.r_[(1, s[:(-1)])]
        elif (weight != None):
            s = weight(tind)
        else:
            s = np.ones_like(tind)
        if (self.censoring == None):
            for g in groups:
                dk = np.bincount(t[(data[:, self.exog] == g)])
                d = np.bincount(t)
                if (np.max(tind) != len(dk)):
                    dif = ((np.max(tind) - len(dk)) + 1)
                    dk = np.r_[(dk, ([0] * dif))]
                dk = dk[:, list(tind)]
                d = d[:, list(tind)]
                dk = dk.astype(float)
                d = d.astype(float)
                dkSum = np.cumsum(dk)
                dSum = np.cumsum(d)
                dkSum = np.r_[(0, dkSum)]
                dSum = np.r_[(0, dSum)]
                nk = (len(data[(data[:, self.exog] == g)]) - dkSum[:(-1)])
                n = (len(data) - dSum[:(-1)])
                d = d[(n > 1)]
                dk = dk[(n > 1)]
                nk = nk[(n > 1)]
                n = n[(n > 1)]
                s = s[(n > 1)]
                ek = ((nk * d) / n)
                Z.append(np.sum((s * (dk - ek))))
                NK.append(nk)
                N.append(n)
                D.append(d)
        else:
            for g in groups:
                censoring = data[:, self.censoring].astype(float).astype(int)
                reverseCensoring = ((-1) * (censoring - 1))
                censored = np.bincount(t, reverseCensoring)
                ck = np.bincount(t[(data[:, self.exog] == g)], reverseCensoring[(data[:, self.exog] == g)])
                dk = np.bincount(t[(data[:, self.exog] == g)], censoring[(data[:, self.exog] == g)])
                d = np.bincount(t, censoring)
                if (np.max(tind) != len(dk)):
                    dif = ((np.max(tind) - len(dk)) + 1)
                    dk = np.r_[(dk, ([0] * dif))]
                    ck = np.r_[(ck, ([0] * dif))]
                dk = dk[:, list(tind)]
                ck = ck[:, list(tind)]
                d = d[:, list(tind)]
                dk = dk.astype(float)
                d = d.astype(float)
                ck = ck.astype(float)
                dkSum = np.cumsum(dk)
                dSum = np.cumsum(d)
                ck = np.cumsum(ck)
                ck = np.r_[(0, ck)]
                dkSum = np.r_[(0, dkSum)]
                dSum = np.r_[(0, dSum)]
                censored = censored[:, list(tind)]
                censored = censored.astype(float)
                censoredSum = np.cumsum(censored)
                censoredSum = np.r_[(0, censoredSum)]
                nk = ((len(data[(data[:, self.exog] == g)]) - dkSum[:(-1)]) - ck[:(-1)])
                n = ((len(data) - dSum[:(-1)]) - censoredSum[:(-1)])
                d = d[(n > 1)]
                dk = dk[(n > 1)]
                nk = nk[(n > 1)]
                n = n[(n > 1)]
                s = s[(n > 1)]
                ek = ((nk * d) / n)
                Z.append(np.sum((s * (dk - ek))))
                NK.append(nk)
                N.append(n)
                D.append(d)
        Z = np.array(Z)
        N = np.array(N)
        D = np.array(D)
        NK = np.array(NK)
        sigma = ((-1) * np.dot(((((NK / N) * ((N - D) / (N - 1))) * D) * np.array(([(s ** 2)] * len(D)))), np.transpose((NK / N))))
        np.fill_diagonal(sigma, np.diagonal(np.dot(((((NK / N) * ((N - D) / (N - 1))) * D) * np.array(([(s ** 2)] * len(D)))), np.transpose((1 - (NK / N))))))
        chisq = np.dot(np.transpose(Z), np.dot(la.pinv(sigma), Z))
        df = (len(groups) - 1)
        return np.array([chisq, df, stats.chi2.sf(chisq, df)])
    else:
        raise ValueError('groups    must    be    in    column    exog')
