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:listordictionary, 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=19398076Expand 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