Module stikpetP.tests.test_james_owa

Expand source code
import pandas as pd
from scipy.stats import chi2

def ts_james_owa(nomField, scaleField, categories=None, order=2, ddof=2):
    '''
    James One-Way ANOVA
    ----------------------
    Tests if the means (averages) of each category could be the same in the population.
    
    If the p-value is below a pre-defined threshold (usually 0.05), the null hypothesis is rejected, and there are then at least two categories who will have a different mean on the scaleField score in the population.
    
    James (1951) proposed three tests, one for large group sizes, a 'first order test', and a 'second order test'. The later two a significance level (α) is chosen and a critical value is then calculated based on a modification of the chi-square distribution.
    
    The James test statistic value J is the same as the test statistic in Cochran's test, calculated slightly different, but will lead to the same result.
    
    Schneider and Penfield (1997) looked at the Welch, Alexander-Govern and the James test (they ignored the Brown-Forsythe since they found it to perform worse than Welch or James), and concluded: “Under variance heterogeneity, Alexander-Govern’s approximation was not only comparable to the Welch test and the James second-order test but was superior, in certain instances, when coupled with the power results for those tests” (p. 285). 
    
    There are quite some alternatives for this, the stikpet library has Fisher, Welch, James, Box, Scott-Smith, Brown-Forsythe, Alexander-Govern, Mehrotra modified Brown-Forsythe, Hartung-Agac-Makabi, Özdemir-Kurt and Wilcox as options. See the notes from ts_fisher_owa() for some discussion on the differences.
    
    Parameters
    ----------
    nomField : pandas series
        data with categories
    scaleField : pandas series
        data with the scores
    categories : list or dictionary, optional
        the categories to use from catField
    order : {2, 1, 0}, optional
        order of James test to use. 0 will use the large-sample approximation.
    ddof : int, optional
        offset for degrees of freedom. Default is 2.
    
    Returns
    -------
    Dataframe with:
    
    * *n*, the sample size
    * *statistic*, the J test statistic
    * *J critical*, the critical J value
    * *df*, degrees of freedom
    * *p-value*, the p-value (significance)
    * *test*, description of test used
    
    
    Notes
    -----
    The formula used is (James, 1951, p. 324):
    $$J = \\sum_{j=1}^k w_j\\times\\left(\\bar{x}_j - \\bar{y}_w\\right)^2 = \\chi_{Cochran}^2$$
    
    James also wrote the formula in a different way, but with the same result (James, 1951, p. 324):
    $$J = \\sum_{j=1}^k w_j\\times\\bar{x}_j^2 - \\frac{\\left(\\sum_{s=1}^k w_s\\times\\bar{x}_s\\right)^2}{w}$$
    
    With:
    $$\\bar{y}_w = \\sum_{j=1}^k h_j\\times\\bar{x}_j$$
    $$h_j = \\frac{w_j}{w}$$
    $$w = \\sum_{j=1}^k w_j$$
    $$w_j = \\frac{n_j}{s_j^2}$$
    $$s_j^2 = \\frac{\\sum_{i=1}^{n_j}\\left(x_{i,j} - \\bar{x}_j\\right)^2}{n_j - 1}$$
    $$\\bar{x}_j = \\frac{\\sum_{i=1}^{n_j} x_{i,j}}{n_j}$$
    
    **Large Sample (order=0)**
    
    For a large sample the p-value (sig.) can be determined by using a chi-square distribution (James, 1951, p. 324):    
    $$df = k - 1$$
    $$sig. = 1 - \\chi^2\\left(J, df\\right)$$
    
    **First Order (order=1)**
    
    James describes to reject the null hypothesis if \\(J>J_{crit1}\\), where \\(J_{crit1}\\) is determined using:
    $$J_{crit1} = \\chi_{crit}^2\\times\\left(1 + \\frac{3\\times\\chi_{crit}^2 +k+1}{2\\times\\left(k^2-1\\right)}\\times\\lambda\\right)$$
    
    With:
    $$\\lambda = \\sum_{j=1}^k \\frac{\\left(1-h_j\\right)^2}{v_j}$$
    $$\\chi_{crit}^2 = Q\\left(\\chi^2\\left(1 - \\alpha, df\\right)\\right)$$
    $$v_j = n_j - 1$$
    
    This function will do an iterative search to find the approximate p-value.
    
    **Second Order (order=2)**
    
    James describes to reject the null hypothesis if \\(J>J_{crit2}\\), where \\(J_{crit2}\\) is determined using:
    $$J_{crit2} = \\sum_{r=1}^9 a_r$$
    
    With:
    $$ a_1 =  \\chi_{crit}^2 + \\frac{1}{2}\\times\\left(3\\times\\chi_4 + \\chi_2\\right)\\times\\ \\lambda_2 $$
    $$ a_2 =  \\frac{1}{16}\\times\\left(3\\times\\chi_4 + \\chi_2\\right)^2\\times\\left(1 - \\frac{k - 3}{\\chi_{crit}^2}\\right)\\times\\ \\lambda_2^2$$
    $$ a_{3f} = \\frac{1}{2}\\times\\left(3\\times\\chi_4+\\chi_2\\right)$$
    $$ a_{3a} = 8\\times R_{23} - 10\\times R_{22} + 4\\times R_{21} - 6\\times R_{12}^2 + 8\\times R_{12}\\times R_{11} - 4\\times R_{11}^2 $$
    $$ a_{3b} = \\left(2\\times R_{23} - 4\\times R_{22} + 2\\times R_{21} - 2\\times R_{12}^2 + 4\\times R_{12}\\times R_{11} - 2\\times R_{11}^2\\right)\\times\\left(\\chi_2-1\\right) $$
    $$ a_{3c} = \\frac{1}{4}\\times\\left(-R_{12}^2 + 4\\times R_{12}\\times R_{11} - 2\\times R_{12}\\times R_{10} - 4\\times R_{11}^2 + 4\\times R_{11}\\times R_{10} - R_{10}^2\\right)\\times\\left(3\\times\\chi_4 - 2\\times\\chi_2 - 1\\right) $$
    $$ a_3 = a_{3f}\\times\\left(a_{3a} + a_{3b} + a_{3c} \\right) $$
    $$ a_4 = \\left(R_{23} - 3\\times R_{22} + 3\\times R_{21} - R_{20}\\right)\\times\\left(5\\times \\chi_6 + 2\\times\\chi_4 + \\chi_2\\right) $$
    $$ a_5 = \\frac{3}{16}\\times\\left(R_{12}^2 - 4\\times R_{23} + 6\\times R_{22} - 4\\times R_{21} + R_{20}\\right)\\times\\left(35\\times\\chi_8 + 15\\times\\chi_6 + 9\\times\\chi_4 + 5\\times\\chi_2\\right) $$
    $$ a_6 = \\frac{1}{16}\\times\\left(-2\\times R_{22}^2 + 4\\times R_{21} - R_{20} + 2\\times R_{12}\\times R_{10} - 4\\times R_{11}\\times R_{10} + R_{10}^2\\right)\\times\\left(9\\times\\chi_8 - 3\\times\\chi_6 - 5\\times\\chi_4 - \\chi_2\\right) $$
    $$ a_7 = \\frac{1}{16}\\times\\left(-2\\times R_{22}^2 + 4\\times R_{21} - R_20 + 2\\times R_{12}\\times R_{10} - 4\\times R_{11}\\times R_{10} + R_{10}^2\\right)\\times\\left(9\\times\\chi_8 - 3\\times\\chi_6 - 5\\times\\chi_4 - \\chi_2\\right) $$
    $$ a_8 = \\frac{1}{4}\\times\\left(-R_{22} + R_{11}^2\\right)\\times\\left(27\\times\\chi_8 + 3\\times\\chi_6 + \\chi_4 + \\chi_2\\right) $$
    $$ a_9 = \\frac{1}{4}\\times\\left(R_{23} - R_{12}\\times R_{11}\\right)\\times\\left(45\\times\\chi_8 + 9\\times\\chi_6 + 7\\times\\chi_4 + 3\\times\\chi_2\\right) $$
    $$ \\lambda_2 = \\sum_{j=1}^k \\frac{\\left(1 - h_j\\right)^2}{v_j^*} $$
    $$ v_j^* = n_j - 2 $$
    $$ R_{xy} = \\sum_{j=1}^k \\frac{1}{\\left(v_j^*\\right)^x}\\times\\left(\\frac{w_j}{w}\\right)^y $$
    $$ \\chi_{2\\times r} = \\frac{\\left(\\chi_{crit}^2\\right)^r}{\\prod_{i=1}^r \\left(k + 2\\times i - 3\\right)} $$
    $$ \\chi_{crit}^2 = \\chi_{crit}^2\\left(1 - \\alpha, df\\right) $$
    
    This function will do an iterative search to find the approximate p-value.
    
    Note that the formula \\(v_j^* = n_j - 2\\) for the second-order test is based on: "..., we finally obtain, as the approximation of order -2 in the,..." (James, 1951, p. 328). It can als be found in Deshon and Alexander (1994, p. 331). However, other authors use in the calculation, for example Myers (1998, p. 209) and Cribbie et al. (2002, p. 62).
    
    By setting the *variation* the function will use \\(v_j^* = n_j - variation\\)
    
    *Symbols used:*
    
    * \\(x_{i,j}\\), the i-th score in category j
    * \\(k\\), the number of categories
    * \\(n_j\\), the sample size of category j
    * \\(\\bar{x}_j\\), the sample mean of category j
    * \\(s_j^2\\), the sample variance of the scores in category j
    * \\(w_j\\), the weight for category j
    * \\(h_j\\), the adjusted weight for category j
    * \\(df\\), the degrees of freedom. 
    * \\(Q\\left(\\chi^2\\left(\\dots\\right)\\right)\\), the inverse cumulative chi-square distribution.
    
    References
    ----------
    Cribbie, R. A., Fiksenbaum, L., Keselman, H. J., & Wilcox, R. R. (2012). Effect of non-normality on test statistics for one-way independent groups designs. *British Journal of Mathematical and Statistical Psychology, 65*(1), 56–73. doi:10.1111/j.2044-8317.2011.02014.x
    
    Deshon, R. P., & Alexander, R. A. (1994). A generalization of James’s second-order approximation to the test for regression slope equality. *Educational and Psychological Measurement, 54*(2), 328–335. doi:10.1177/0013164494054002007
    
    James, G. S. (1951). The comparison of several groups of observations when the ratios of the population variances are unknown. *Biometrika, 38*(3–4), 324–329. doi:10.1093/biomet/38.3-4.324
    
    Myers, L. (1998). Comparability of the james’ second-order approximation test and the alexander and govern  A  statistic for non-normal heteroscedastic data. *Journal of Statistical Computation and Simulation, 60*(3), 207–222. doi:10.1080/00949659808811888
    
    Schneider, P. J., & Penfield, D. A. (1997). Alexander and Govern’s approximation: Providing an alternative to ANOVA under variance heterogeneity. *The Journal of Experimental Education, 65*(3), 271–286. doi:10.1080/00220973.1997.9943459

    Author
    ------
    Made by P. Stikker
    
    Companion website: https://PeterStatistics.com  
    YouTube channel: https://www.youtube.com/stikpet  
    Donations: https://www.patreon.com/bePatron?u=19398076
    
    '''
    if type(nomField) == list:
        nomField = pd.Series(nomField)
        
    if type(scaleField) == list:
        scaleField = pd.Series(scaleField)
        
    data = pd.concat([nomField, scaleField], axis=1)
    data.columns = ["category", "score"]
    
    #remove unused categories
    if categories is not None:
        data = data[data.category.isin(categories)]
    
    #Remove rows with missing values and reset index
    data = data.dropna()    
    data.reset_index()
    
    #overall n, mean and ss
    n = len(data["category"])
    m = data.score.mean()
    sst = data.score.var()*(n-1)
    
    #sample sizes, variances and means per category
    nj = data.groupby('category').count()
    sj2 = data.groupby('category').var()
    mj = data.groupby('category').mean()
    
    #number of categories
    k = len(mj)
    
    wj = nj / sj2
    w = float(wj.sum())
    hj = wj/w
    yw = float((hj*mj).sum())
    
    Jstat = float((wj*(mj - yw)**2).sum())
    
    if order==0:
        testUsed = "James large-sample approximation"
        
        df = k - 1
        pVal = chi2.sf(Jstat, df)
        res = pd.DataFrame([[n, Jstat, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "df", "p-value", "test"]
    
    elif order==1:
        testUsed = "James first-order"
        vj = nj - 1
        lm = float(((1 - hj)**2/vj).sum())
        
        df = k - 1
        pLow=0
        pHigh=1
        pVal=.05
        nIter = 1
        
        c = chi2.ppf(1-pVal, df)
        jCrit = c*(1+(3*c+k+1)/(2*(k**2-1))*lm)
        jCritOriginal = jCrit
        while jCrit != Jstat and nIter < 100:
            if jCrit < Jstat:
                pHigh = pVal
                pVal = (pLow + pVal)/2
            elif jCrit > Jstat:
                pLow = pVal
                pVal = (pHigh + pVal)/2
                
            c = chi2.ppf(1-pVal, df)
            jCrit = c*(1+(3*c+k+1)/(2*(k**2-1))*lm)
            
            nIter = nIter + 1
            
        res = pd.DataFrame([[n, Jstat, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "df", "p-value", "test"]
    
    elif order==2:
        testUsed = "James second-order"
        vj = nj-ddof
        lm = float(((1 - hj)**2/vj).sum())
        
        #R values
        R10 = float((hj**0/(vj**1)).sum())
        R11 = float((hj**1/(vj**1)).sum())
        R12 = float((hj**2/(vj**1)).sum())
        R20 = float((hj**0/(vj**2)).sum())
        R21 = float((hj**1/(vj**2)).sum())
        R22 = float((hj**2/(vj**2)).sum())
        R23 = float((hj**3/(vj**2)).sum())
        
        #determine denominator for chi-variables
        chi2den = k + 2 * 1 - 3
        chi4den = chi2den * (k + 2 * 2 - 3)
        chi6den = chi4den * (k + 2 * 3 - 3)
        chi8den = chi6den * (k + 2 * 4 - 3)
        
        df = k - 1
        pLow=0
        pHigh=1
        pVal=.05
        nIter = 1
        loop = True
        original = True
        
        while loop:
            #critical chi-square value
            c = chi2.ppf(1-pVal, df)
            
            #chi-variables
            chi2v = c / chi2den
            chi4 = c**2 / chi4den
            chi6 = c**3 / chi6den
            chi8 = c**4 / chi8den

            prt1 = c + 1 / 2 * (3 * chi4 + chi2v) * lm + 1 / 16 * (3 * chi4 + chi2v)**2 * (1 - (k - 3) / c) * lm**2
            prt2a = 1 / 2 * (3 * chi4 + chi2v)
            prt2b = (8 * R23 - 10 * R22 + 4 * R21 - 6 * R12 **2 + 8 * R12 * R11 - 4 * R11 **2)
            prt2c = (2 * R23 - 4 * R22 + 2 * R21 - 2 * R12**2 + 4 * R12 * R11 - 2 * R11**2) * (chi2v - 1)
            prt2d = 1 / 4 * (-1 * R12 **2 + 4 * R12 * R11 - 2 * R12 * R10 - 4 * R11 **2 + 4 * R11 * R10 - R10 **2) * (3 * chi4 - 2 * chi2v - 1)
            prt2 = prt2a * (prt2b + prt2c + prt2d)
            prt3 = (R23 - 3 * R22 + 3 * R21 - R20) * (5 * chi6 + 2 * chi4 + chi2v)
            prt4 = 3 / 16 * (R12**2 - 4 * R23 + 6 * R22 - 4 * R21 + R20) * (35 * chi8 + 15 * chi6 + 9 * chi4 + 5 * chi2v)
            prt5 = 1 / 16 * (-2 * R22 + 4 * R21 - R20 + 2 * R12 * R10 - 4 * R11 * R10 + R10**2) * (9 * chi8 - 3 * chi6 - 5 * chi4 - chi2v)
            prt6 = 1 / 4 * (-1 * R22 + R11**2) * (27 * chi8 + 3 * chi6 + chi4 + chi2v)
            prt7 = 1 / 4 * (R23 - R12 * R11) * (45 * chi8 + 9 * chi6 + 7 * chi4 + 3 * chi2v)

            Jcrit = prt1 + prt2 + prt3 + prt4 + prt5 + prt6 + prt7
            
            if original:
                JcritOr = Jcrit
                original = False
            
            if Jcrit < Jstat:
                pHigh = pVal
                pVal = (pLow + pVal) / 2
            elif Jcrit > Jstat:
                pLow = pVal
                pVal = (pHigh + pVal) / 2

            nIter = nIter + 1
            
            if nIter > 100 or Jcrit==Jstat:
                loop = False
        
        res = pd.DataFrame([[n, Jstat, JcritOr, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "J critical", "df", "p-value", "test"]
        
    
    return res

Functions

def ts_james_owa(nomField, scaleField, categories=None, order=2, ddof=2)

James One-Way ANOVA

Tests if the means (averages) of each category could be the same in the population.

If the p-value is below a pre-defined threshold (usually 0.05), the null hypothesis is rejected, and there are then at least two categories who will have a different mean on the scaleField score in the population.

James (1951) proposed three tests, one for large group sizes, a 'first order test', and a 'second order test'. The later two a significance level (α) is chosen and a critical value is then calculated based on a modification of the chi-square distribution.

The James test statistic value J is the same as the test statistic in Cochran's test, calculated slightly different, but will lead to the same result.

Schneider and Penfield (1997) looked at the Welch, Alexander-Govern and the James test (they ignored the Brown-Forsythe since they found it to perform worse than Welch or James), and concluded: “Under variance heterogeneity, Alexander-Govern’s approximation was not only comparable to the Welch test and the James second-order test but was superior, in certain instances, when coupled with the power results for those tests” (p. 285).

There are quite some alternatives for this, the stikpet library has Fisher, Welch, James, Box, Scott-Smith, Brown-Forsythe, Alexander-Govern, Mehrotra modified Brown-Forsythe, Hartung-Agac-Makabi, Özdemir-Kurt and Wilcox as options. See the notes from ts_fisher_owa() for some discussion on the differences.

Parameters

nomField : pandas series
data with categories
scaleField : pandas series
data with the scores
categories : list or dictionary, optional
the categories to use from catField
order : {2, 1, 0}, optional
order of James test to use. 0 will use the large-sample approximation.
ddof : int, optional
offset for degrees of freedom. Default is 2.

Returns

Dataframe with:
 
  • n, the sample size
  • statistic, the J test statistic
  • J critical, the critical J value
  • df, degrees of freedom
  • p-value, the p-value (significance)
  • test, description of test used

Notes

The formula used is (James, 1951, p. 324): J = \sum_{j=1}^k w_j\times\left(\bar{x}_j - \bar{y}_w\right)^2 = \chi_{Cochran}^2

James also wrote the formula in a different way, but with the same result (James, 1951, p. 324): J = \sum_{j=1}^k w_j\times\bar{x}_j^2 - \frac{\left(\sum_{s=1}^k w_s\times\bar{x}_s\right)^2}{w}

With: \bar{y}_w = \sum_{j=1}^k h_j\times\bar{x}_j h_j = \frac{w_j}{w} w = \sum_{j=1}^k w_j w_j = \frac{n_j}{s_j^2} s_j^2 = \frac{\sum_{i=1}^{n_j}\left(x_{i,j} - \bar{x}_j\right)^2}{n_j - 1} \bar{x}_j = \frac{\sum_{i=1}^{n_j} x_{i,j}}{n_j}

Large Sample (order=0)

For a large sample the p-value (sig.) can be determined by using a chi-square distribution (James, 1951, p. 324):
df = k - 1 sig. = 1 - \chi^2\left(J, df\right)

First Order (order=1)

James describes to reject the null hypothesis if J>J_{crit1}, where J_{crit1} is determined using: J_{crit1} = \chi_{crit}^2\times\left(1 + \frac{3\times\chi_{crit}^2 +k+1}{2\times\left(k^2-1\right)}\times\lambda\right)

With: \lambda = \sum_{j=1}^k \frac{\left(1-h_j\right)^2}{v_j} \chi_{crit}^2 = Q\left(\chi^2\left(1 - \alpha, df\right)\right) v_j = n_j - 1

This function will do an iterative search to find the approximate p-value.

Second Order (order=2)

James describes to reject the null hypothesis if J>J_{crit2}, where J_{crit2} is determined using: J_{crit2} = \sum_{r=1}^9 a_r

With: a_1 = \chi_{crit}^2 + \frac{1}{2}\times\left(3\times\chi_4 + \chi_2\right)\times\ \lambda_2 a_2 = \frac{1}{16}\times\left(3\times\chi_4 + \chi_2\right)^2\times\left(1 - \frac{k - 3}{\chi_{crit}^2}\right)\times\ \lambda_2^2 a_{3f} = \frac{1}{2}\times\left(3\times\chi_4+\chi_2\right) a_{3a} = 8\times R_{23} - 10\times R_{22} + 4\times R_{21} - 6\times R_{12}^2 + 8\times R_{12}\times R_{11} - 4\times R_{11}^2 a_{3b} = \left(2\times R_{23} - 4\times R_{22} + 2\times R_{21} - 2\times R_{12}^2 + 4\times R_{12}\times R_{11} - 2\times R_{11}^2\right)\times\left(\chi_2-1\right) a_{3c} = \frac{1}{4}\times\left(-R_{12}^2 + 4\times R_{12}\times R_{11} - 2\times R_{12}\times R_{10} - 4\times R_{11}^2 + 4\times R_{11}\times R_{10} - R_{10}^2\right)\times\left(3\times\chi_4 - 2\times\chi_2 - 1\right) a_3 = a_{3f}\times\left(a_{3a} + a_{3b} + a_{3c} \right) a_4 = \left(R_{23} - 3\times R_{22} + 3\times R_{21} - R_{20}\right)\times\left(5\times \chi_6 + 2\times\chi_4 + \chi_2\right) a_5 = \frac{3}{16}\times\left(R_{12}^2 - 4\times R_{23} + 6\times R_{22} - 4\times R_{21} + R_{20}\right)\times\left(35\times\chi_8 + 15\times\chi_6 + 9\times\chi_4 + 5\times\chi_2\right) a_6 = \frac{1}{16}\times\left(-2\times R_{22}^2 + 4\times R_{21} - R_{20} + 2\times R_{12}\times R_{10} - 4\times R_{11}\times R_{10} + R_{10}^2\right)\times\left(9\times\chi_8 - 3\times\chi_6 - 5\times\chi_4 - \chi_2\right) a_7 = \frac{1}{16}\times\left(-2\times R_{22}^2 + 4\times R_{21} - R_20 + 2\times R_{12}\times R_{10} - 4\times R_{11}\times R_{10} + R_{10}^2\right)\times\left(9\times\chi_8 - 3\times\chi_6 - 5\times\chi_4 - \chi_2\right) a_8 = \frac{1}{4}\times\left(-R_{22} + R_{11}^2\right)\times\left(27\times\chi_8 + 3\times\chi_6 + \chi_4 + \chi_2\right) a_9 = \frac{1}{4}\times\left(R_{23} - R_{12}\times R_{11}\right)\times\left(45\times\chi_8 + 9\times\chi_6 + 7\times\chi_4 + 3\times\chi_2\right) \lambda_2 = \sum_{j=1}^k \frac{\left(1 - h_j\right)^2}{v_j^*} v_j^* = n_j - 2 R_{xy} = \sum_{j=1}^k \frac{1}{\left(v_j^*\right)^x}\times\left(\frac{w_j}{w}\right)^y \chi_{2\times r} = \frac{\left(\chi_{crit}^2\right)^r}{\prod_{i=1}^r \left(k + 2\times i - 3\right)} \chi_{crit}^2 = \chi_{crit}^2\left(1 - \alpha, df\right)

This function will do an iterative search to find the approximate p-value.

Note that the formula v_j^* = n_j - 2 for the second-order test is based on: "…, we finally obtain, as the approximation of order -2 in the,…" (James, 1951, p. 328). It can als be found in Deshon and Alexander (1994, p. 331). However, other authors use in the calculation, for example Myers (1998, p. 209) and Cribbie et al. (2002, p. 62).

By setting the variation the function will use v_j^* = n_j - variation

Symbols used:

  • x_{i,j}, the i-th score in category j
  • k, the number of categories
  • n_j, the sample size of category j
  • \bar{x}_j, the sample mean of category j
  • s_j^2, the sample variance of the scores in category j
  • w_j, the weight for category j
  • h_j, the adjusted weight for category j
  • df, the degrees of freedom.
  • Q\left(\chi^2\left(\dots\right)\right), the inverse cumulative chi-square distribution.

References

Cribbie, R. A., Fiksenbaum, L., Keselman, H. J., & Wilcox, R. R. (2012). Effect of non-normality on test statistics for one-way independent groups designs. British Journal of Mathematical and Statistical Psychology, 65(1), 56–73. doi:10.1111/j.2044-8317.2011.02014.x

Deshon, R. P., & Alexander, R. A. (1994). A generalization of James’s second-order approximation to the test for regression slope equality. Educational and Psychological Measurement, 54(2), 328–335. doi:10.1177/0013164494054002007

James, G. S. (1951). The comparison of several groups of observations when the ratios of the population variances are unknown. Biometrika, 38(3–4), 324–329. doi:10.1093/biomet/38.3-4.324

Myers, L. (1998). Comparability of the james’ second-order approximation test and the alexander and govern A statistic for non-normal heteroscedastic data. Journal of Statistical Computation and Simulation, 60(3), 207–222. doi:10.1080/00949659808811888

Schneider, P. J., & Penfield, D. A. (1997). Alexander and Govern’s approximation: Providing an alternative to ANOVA under variance heterogeneity. The Journal of Experimental Education, 65(3), 271–286. doi:10.1080/00220973.1997.9943459

Author

Made by P. Stikker

Companion website: https://PeterStatistics.com
YouTube channel: https://www.youtube.com/stikpet
Donations: https://www.patreon.com/bePatron?u=19398076

Expand source code
def ts_james_owa(nomField, scaleField, categories=None, order=2, ddof=2):
    '''
    James One-Way ANOVA
    ----------------------
    Tests if the means (averages) of each category could be the same in the population.
    
    If the p-value is below a pre-defined threshold (usually 0.05), the null hypothesis is rejected, and there are then at least two categories who will have a different mean on the scaleField score in the population.
    
    James (1951) proposed three tests, one for large group sizes, a 'first order test', and a 'second order test'. The later two a significance level (α) is chosen and a critical value is then calculated based on a modification of the chi-square distribution.
    
    The James test statistic value J is the same as the test statistic in Cochran's test, calculated slightly different, but will lead to the same result.
    
    Schneider and Penfield (1997) looked at the Welch, Alexander-Govern and the James test (they ignored the Brown-Forsythe since they found it to perform worse than Welch or James), and concluded: “Under variance heterogeneity, Alexander-Govern’s approximation was not only comparable to the Welch test and the James second-order test but was superior, in certain instances, when coupled with the power results for those tests” (p. 285). 
    
    There are quite some alternatives for this, the stikpet library has Fisher, Welch, James, Box, Scott-Smith, Brown-Forsythe, Alexander-Govern, Mehrotra modified Brown-Forsythe, Hartung-Agac-Makabi, Özdemir-Kurt and Wilcox as options. See the notes from ts_fisher_owa() for some discussion on the differences.
    
    Parameters
    ----------
    nomField : pandas series
        data with categories
    scaleField : pandas series
        data with the scores
    categories : list or dictionary, optional
        the categories to use from catField
    order : {2, 1, 0}, optional
        order of James test to use. 0 will use the large-sample approximation.
    ddof : int, optional
        offset for degrees of freedom. Default is 2.
    
    Returns
    -------
    Dataframe with:
    
    * *n*, the sample size
    * *statistic*, the J test statistic
    * *J critical*, the critical J value
    * *df*, degrees of freedom
    * *p-value*, the p-value (significance)
    * *test*, description of test used
    
    
    Notes
    -----
    The formula used is (James, 1951, p. 324):
    $$J = \\sum_{j=1}^k w_j\\times\\left(\\bar{x}_j - \\bar{y}_w\\right)^2 = \\chi_{Cochran}^2$$
    
    James also wrote the formula in a different way, but with the same result (James, 1951, p. 324):
    $$J = \\sum_{j=1}^k w_j\\times\\bar{x}_j^2 - \\frac{\\left(\\sum_{s=1}^k w_s\\times\\bar{x}_s\\right)^2}{w}$$
    
    With:
    $$\\bar{y}_w = \\sum_{j=1}^k h_j\\times\\bar{x}_j$$
    $$h_j = \\frac{w_j}{w}$$
    $$w = \\sum_{j=1}^k w_j$$
    $$w_j = \\frac{n_j}{s_j^2}$$
    $$s_j^2 = \\frac{\\sum_{i=1}^{n_j}\\left(x_{i,j} - \\bar{x}_j\\right)^2}{n_j - 1}$$
    $$\\bar{x}_j = \\frac{\\sum_{i=1}^{n_j} x_{i,j}}{n_j}$$
    
    **Large Sample (order=0)**
    
    For a large sample the p-value (sig.) can be determined by using a chi-square distribution (James, 1951, p. 324):    
    $$df = k - 1$$
    $$sig. = 1 - \\chi^2\\left(J, df\\right)$$
    
    **First Order (order=1)**
    
    James describes to reject the null hypothesis if \\(J>J_{crit1}\\), where \\(J_{crit1}\\) is determined using:
    $$J_{crit1} = \\chi_{crit}^2\\times\\left(1 + \\frac{3\\times\\chi_{crit}^2 +k+1}{2\\times\\left(k^2-1\\right)}\\times\\lambda\\right)$$
    
    With:
    $$\\lambda = \\sum_{j=1}^k \\frac{\\left(1-h_j\\right)^2}{v_j}$$
    $$\\chi_{crit}^2 = Q\\left(\\chi^2\\left(1 - \\alpha, df\\right)\\right)$$
    $$v_j = n_j - 1$$
    
    This function will do an iterative search to find the approximate p-value.
    
    **Second Order (order=2)**
    
    James describes to reject the null hypothesis if \\(J>J_{crit2}\\), where \\(J_{crit2}\\) is determined using:
    $$J_{crit2} = \\sum_{r=1}^9 a_r$$
    
    With:
    $$ a_1 =  \\chi_{crit}^2 + \\frac{1}{2}\\times\\left(3\\times\\chi_4 + \\chi_2\\right)\\times\\ \\lambda_2 $$
    $$ a_2 =  \\frac{1}{16}\\times\\left(3\\times\\chi_4 + \\chi_2\\right)^2\\times\\left(1 - \\frac{k - 3}{\\chi_{crit}^2}\\right)\\times\\ \\lambda_2^2$$
    $$ a_{3f} = \\frac{1}{2}\\times\\left(3\\times\\chi_4+\\chi_2\\right)$$
    $$ a_{3a} = 8\\times R_{23} - 10\\times R_{22} + 4\\times R_{21} - 6\\times R_{12}^2 + 8\\times R_{12}\\times R_{11} - 4\\times R_{11}^2 $$
    $$ a_{3b} = \\left(2\\times R_{23} - 4\\times R_{22} + 2\\times R_{21} - 2\\times R_{12}^2 + 4\\times R_{12}\\times R_{11} - 2\\times R_{11}^2\\right)\\times\\left(\\chi_2-1\\right) $$
    $$ a_{3c} = \\frac{1}{4}\\times\\left(-R_{12}^2 + 4\\times R_{12}\\times R_{11} - 2\\times R_{12}\\times R_{10} - 4\\times R_{11}^2 + 4\\times R_{11}\\times R_{10} - R_{10}^2\\right)\\times\\left(3\\times\\chi_4 - 2\\times\\chi_2 - 1\\right) $$
    $$ a_3 = a_{3f}\\times\\left(a_{3a} + a_{3b} + a_{3c} \\right) $$
    $$ a_4 = \\left(R_{23} - 3\\times R_{22} + 3\\times R_{21} - R_{20}\\right)\\times\\left(5\\times \\chi_6 + 2\\times\\chi_4 + \\chi_2\\right) $$
    $$ a_5 = \\frac{3}{16}\\times\\left(R_{12}^2 - 4\\times R_{23} + 6\\times R_{22} - 4\\times R_{21} + R_{20}\\right)\\times\\left(35\\times\\chi_8 + 15\\times\\chi_6 + 9\\times\\chi_4 + 5\\times\\chi_2\\right) $$
    $$ a_6 = \\frac{1}{16}\\times\\left(-2\\times R_{22}^2 + 4\\times R_{21} - R_{20} + 2\\times R_{12}\\times R_{10} - 4\\times R_{11}\\times R_{10} + R_{10}^2\\right)\\times\\left(9\\times\\chi_8 - 3\\times\\chi_6 - 5\\times\\chi_4 - \\chi_2\\right) $$
    $$ a_7 = \\frac{1}{16}\\times\\left(-2\\times R_{22}^2 + 4\\times R_{21} - R_20 + 2\\times R_{12}\\times R_{10} - 4\\times R_{11}\\times R_{10} + R_{10}^2\\right)\\times\\left(9\\times\\chi_8 - 3\\times\\chi_6 - 5\\times\\chi_4 - \\chi_2\\right) $$
    $$ a_8 = \\frac{1}{4}\\times\\left(-R_{22} + R_{11}^2\\right)\\times\\left(27\\times\\chi_8 + 3\\times\\chi_6 + \\chi_4 + \\chi_2\\right) $$
    $$ a_9 = \\frac{1}{4}\\times\\left(R_{23} - R_{12}\\times R_{11}\\right)\\times\\left(45\\times\\chi_8 + 9\\times\\chi_6 + 7\\times\\chi_4 + 3\\times\\chi_2\\right) $$
    $$ \\lambda_2 = \\sum_{j=1}^k \\frac{\\left(1 - h_j\\right)^2}{v_j^*} $$
    $$ v_j^* = n_j - 2 $$
    $$ R_{xy} = \\sum_{j=1}^k \\frac{1}{\\left(v_j^*\\right)^x}\\times\\left(\\frac{w_j}{w}\\right)^y $$
    $$ \\chi_{2\\times r} = \\frac{\\left(\\chi_{crit}^2\\right)^r}{\\prod_{i=1}^r \\left(k + 2\\times i - 3\\right)} $$
    $$ \\chi_{crit}^2 = \\chi_{crit}^2\\left(1 - \\alpha, df\\right) $$
    
    This function will do an iterative search to find the approximate p-value.
    
    Note that the formula \\(v_j^* = n_j - 2\\) for the second-order test is based on: "..., we finally obtain, as the approximation of order -2 in the,..." (James, 1951, p. 328). It can als be found in Deshon and Alexander (1994, p. 331). However, other authors use in the calculation, for example Myers (1998, p. 209) and Cribbie et al. (2002, p. 62).
    
    By setting the *variation* the function will use \\(v_j^* = n_j - variation\\)
    
    *Symbols used:*
    
    * \\(x_{i,j}\\), the i-th score in category j
    * \\(k\\), the number of categories
    * \\(n_j\\), the sample size of category j
    * \\(\\bar{x}_j\\), the sample mean of category j
    * \\(s_j^2\\), the sample variance of the scores in category j
    * \\(w_j\\), the weight for category j
    * \\(h_j\\), the adjusted weight for category j
    * \\(df\\), the degrees of freedom. 
    * \\(Q\\left(\\chi^2\\left(\\dots\\right)\\right)\\), the inverse cumulative chi-square distribution.
    
    References
    ----------
    Cribbie, R. A., Fiksenbaum, L., Keselman, H. J., & Wilcox, R. R. (2012). Effect of non-normality on test statistics for one-way independent groups designs. *British Journal of Mathematical and Statistical Psychology, 65*(1), 56–73. doi:10.1111/j.2044-8317.2011.02014.x
    
    Deshon, R. P., & Alexander, R. A. (1994). A generalization of James’s second-order approximation to the test for regression slope equality. *Educational and Psychological Measurement, 54*(2), 328–335. doi:10.1177/0013164494054002007
    
    James, G. S. (1951). The comparison of several groups of observations when the ratios of the population variances are unknown. *Biometrika, 38*(3–4), 324–329. doi:10.1093/biomet/38.3-4.324
    
    Myers, L. (1998). Comparability of the james’ second-order approximation test and the alexander and govern  A  statistic for non-normal heteroscedastic data. *Journal of Statistical Computation and Simulation, 60*(3), 207–222. doi:10.1080/00949659808811888
    
    Schneider, P. J., & Penfield, D. A. (1997). Alexander and Govern’s approximation: Providing an alternative to ANOVA under variance heterogeneity. *The Journal of Experimental Education, 65*(3), 271–286. doi:10.1080/00220973.1997.9943459

    Author
    ------
    Made by P. Stikker
    
    Companion website: https://PeterStatistics.com  
    YouTube channel: https://www.youtube.com/stikpet  
    Donations: https://www.patreon.com/bePatron?u=19398076
    
    '''
    if type(nomField) == list:
        nomField = pd.Series(nomField)
        
    if type(scaleField) == list:
        scaleField = pd.Series(scaleField)
        
    data = pd.concat([nomField, scaleField], axis=1)
    data.columns = ["category", "score"]
    
    #remove unused categories
    if categories is not None:
        data = data[data.category.isin(categories)]
    
    #Remove rows with missing values and reset index
    data = data.dropna()    
    data.reset_index()
    
    #overall n, mean and ss
    n = len(data["category"])
    m = data.score.mean()
    sst = data.score.var()*(n-1)
    
    #sample sizes, variances and means per category
    nj = data.groupby('category').count()
    sj2 = data.groupby('category').var()
    mj = data.groupby('category').mean()
    
    #number of categories
    k = len(mj)
    
    wj = nj / sj2
    w = float(wj.sum())
    hj = wj/w
    yw = float((hj*mj).sum())
    
    Jstat = float((wj*(mj - yw)**2).sum())
    
    if order==0:
        testUsed = "James large-sample approximation"
        
        df = k - 1
        pVal = chi2.sf(Jstat, df)
        res = pd.DataFrame([[n, Jstat, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "df", "p-value", "test"]
    
    elif order==1:
        testUsed = "James first-order"
        vj = nj - 1
        lm = float(((1 - hj)**2/vj).sum())
        
        df = k - 1
        pLow=0
        pHigh=1
        pVal=.05
        nIter = 1
        
        c = chi2.ppf(1-pVal, df)
        jCrit = c*(1+(3*c+k+1)/(2*(k**2-1))*lm)
        jCritOriginal = jCrit
        while jCrit != Jstat and nIter < 100:
            if jCrit < Jstat:
                pHigh = pVal
                pVal = (pLow + pVal)/2
            elif jCrit > Jstat:
                pLow = pVal
                pVal = (pHigh + pVal)/2
                
            c = chi2.ppf(1-pVal, df)
            jCrit = c*(1+(3*c+k+1)/(2*(k**2-1))*lm)
            
            nIter = nIter + 1
            
        res = pd.DataFrame([[n, Jstat, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "df", "p-value", "test"]
    
    elif order==2:
        testUsed = "James second-order"
        vj = nj-ddof
        lm = float(((1 - hj)**2/vj).sum())
        
        #R values
        R10 = float((hj**0/(vj**1)).sum())
        R11 = float((hj**1/(vj**1)).sum())
        R12 = float((hj**2/(vj**1)).sum())
        R20 = float((hj**0/(vj**2)).sum())
        R21 = float((hj**1/(vj**2)).sum())
        R22 = float((hj**2/(vj**2)).sum())
        R23 = float((hj**3/(vj**2)).sum())
        
        #determine denominator for chi-variables
        chi2den = k + 2 * 1 - 3
        chi4den = chi2den * (k + 2 * 2 - 3)
        chi6den = chi4den * (k + 2 * 3 - 3)
        chi8den = chi6den * (k + 2 * 4 - 3)
        
        df = k - 1
        pLow=0
        pHigh=1
        pVal=.05
        nIter = 1
        loop = True
        original = True
        
        while loop:
            #critical chi-square value
            c = chi2.ppf(1-pVal, df)
            
            #chi-variables
            chi2v = c / chi2den
            chi4 = c**2 / chi4den
            chi6 = c**3 / chi6den
            chi8 = c**4 / chi8den

            prt1 = c + 1 / 2 * (3 * chi4 + chi2v) * lm + 1 / 16 * (3 * chi4 + chi2v)**2 * (1 - (k - 3) / c) * lm**2
            prt2a = 1 / 2 * (3 * chi4 + chi2v)
            prt2b = (8 * R23 - 10 * R22 + 4 * R21 - 6 * R12 **2 + 8 * R12 * R11 - 4 * R11 **2)
            prt2c = (2 * R23 - 4 * R22 + 2 * R21 - 2 * R12**2 + 4 * R12 * R11 - 2 * R11**2) * (chi2v - 1)
            prt2d = 1 / 4 * (-1 * R12 **2 + 4 * R12 * R11 - 2 * R12 * R10 - 4 * R11 **2 + 4 * R11 * R10 - R10 **2) * (3 * chi4 - 2 * chi2v - 1)
            prt2 = prt2a * (prt2b + prt2c + prt2d)
            prt3 = (R23 - 3 * R22 + 3 * R21 - R20) * (5 * chi6 + 2 * chi4 + chi2v)
            prt4 = 3 / 16 * (R12**2 - 4 * R23 + 6 * R22 - 4 * R21 + R20) * (35 * chi8 + 15 * chi6 + 9 * chi4 + 5 * chi2v)
            prt5 = 1 / 16 * (-2 * R22 + 4 * R21 - R20 + 2 * R12 * R10 - 4 * R11 * R10 + R10**2) * (9 * chi8 - 3 * chi6 - 5 * chi4 - chi2v)
            prt6 = 1 / 4 * (-1 * R22 + R11**2) * (27 * chi8 + 3 * chi6 + chi4 + chi2v)
            prt7 = 1 / 4 * (R23 - R12 * R11) * (45 * chi8 + 9 * chi6 + 7 * chi4 + 3 * chi2v)

            Jcrit = prt1 + prt2 + prt3 + prt4 + prt5 + prt6 + prt7
            
            if original:
                JcritOr = Jcrit
                original = False
            
            if Jcrit < Jstat:
                pHigh = pVal
                pVal = (pLow + pVal) / 2
            elif Jcrit > Jstat:
                pLow = pVal
                pVal = (pHigh + pVal) / 2

            nIter = nIter + 1
            
            if nIter > 100 or Jcrit==Jstat:
                loop = False
        
        res = pd.DataFrame([[n, Jstat, JcritOr, df, pVal, testUsed]])
        res.columns = ["n", "statistic", "J critical", "df", "p-value", "test"]
        
    
    return res