1/n - Proportional Correction to Tests of Independence

Jan Vrbik

Mathematics Department, Brock University, St. Catharines, Canada.

DOI: 10.4236/am.2025.163013                

Abstract

It is known that both the Pearson and G-test of independence have the same asymptotic distribution, namely ${\chi }^2$ with $\left(K-1\right)\left(M-1\right)$ degrees of freedom, where $K$ and $M$ are the number of levels, respectively, of two attributes. However, when extending the accuracy ot this approximation by $\frac{1}{n}$ -proportional terms, the resulting corrections differ quite dramatically. The purpose of this article is to derive each of these corrections and demonstrate their use in practical application.

Keywords

Test of Independence, Asymptotic Distribution, Chi-Square Test, Contingency Tables, G-Test, Kronecker Product

1. Introduction

When testing whether two nominal-scale attributes (one having $K$ discrete values, the other one $M$ ) are independent, the traditional (Pearson’s) test uses the following test statistic

$U:={\displaystyle \sum }_{i=1}^K{\displaystyle \sum }_{j=1}^M\frac{{\left(X_{i,j}-\frac{X_{i,•}{X}_{•,j}}{n}\right)}^2}{\frac{X_{i,•}{X}_{•,j}}{n}}$ (1)

where $n$ is the total number of observations, $X_{i,j}$ is the number of those resulting in the $i^{\text{th}}$ value of the first attribute and in the $j^{\text{th}}$ value of the second attribute, and The null hypothesis claims that the probability of an observation being of the ${\left(i,j\right)}^{\text{th}}$ type is $p_i{q}_j$ where ${\displaystyle {\sum }_{j=1}^M{p}_i}=1$ and ${\displaystyle {\sum }_{i=1}^K{q}_j}=1$ ; note that the multinomial probability mass function of the $X_{i.j}$ ’s is therefore

$f\left(x\right)=\left(\begin{array}{c}n\\ x_{1,1},x_{1,2},\cdots ,x_{K,M}\end{array}\right){\left(p_1{q}_1\right)}^{x_{1,1}}{\left(p_1{q}_2\right)}^{x_{1,2}}\cdots {\left(p_K{q}_M\right)}^{x_{K,M}}$ (2)

having $K+M-2$ algebraically independent parameters [1].

This model implies that the conditional probability of an observation’s second attribute being equal to $j$ , given that the value of its first attribute is $i$ , is given by $q_j$ (for any $j$ and $i$ ), i.e. that the attributes are independent of each other. The test also assumes that the $p_i$ and $q_j$ probabilities are unknown and are estimated by their maximum-likelihood estimators, namely $\frac{X_{i,•}}{n}$ and $\frac{X_{•,j}}{n}$ respectively. The alternate hypothesis allows the probabilities of the $X_{i,j}$ ’s to have any non-negative values, as long as these add up to 1, thus resulting in $KM-1$ algebraically independent parameters. Failing the null hypothesis makes the above test statistic increase in value; the critical region is always an upper tail of the corresponding distribution.

Another test statistic proposed decades later for testing the same hypothesis is (using the same notation, for future convenience)

(3)

of the so-called G test [2] (note that when $X_{i.j}=0$ , which happens occasionally, makes the corresponding term of the last expression equal to 0). It is well known (something we verify shortly) that, in the $n\to \infty$ limit, both Formula (1) and Formula (3) converge to the same expression, whose asymptotic distribution is chi-squared with $\left(K-1\right)\left(M-1\right)$ degrees of freedom [3].

The main objective of this article is to find $\frac{1}{n}$ -proportional correction to this distribution, separately for each of our two test statistics.

Next comes a brief review of the basic mathematical tool used throughout this article.

2. Kronecker Product

In this section $a,$ $b,$ $c$ and $d$ are matrices of any dimensions (a dimension equal to 1 implies a column or row vector; both dimensions equal to 1 represents a scalar). Kronecker product $a\otimes b$ of two matrices is created by multiplying every element of $a$ by the full matrix $b$ and organizing the resulting blocks into a single matrix [4]. This product is clearly non-commutative, associative and distributive over addition; furthermore (assuming conformable dimensions for each matrix multiplication)

${\left(a\otimes b\right)}^{\text{T}}=a^{\text{T}}\otimes b^{\text{T}}$

${\left(a\otimes b\right)}^{-1}=a^{-1}\otimes b^{-1}$

$a\otimes b\circ c\otimes d=a\cdot c\otimes b\cdot d$

where both $\circ$ and $\cdot$ indicate matrix multiplication; note that $\cdot$ takes precedence over $\otimes$ which, in turn, takes precedence over $\circ$ (this is the reason for the duplicate notation). When $a$ and $b$ are square matrices ($K$ by $K$ and $M$ by $M$ respectively)

$\text{det}\left(a\otimes b\right)=\text{det}{\left(a\right)}^M\text{det}{\left(b\right)}^K$

We also need the following Woodbury’s identity [5]

${\left(a+b\cdot c\right)}^{-1}=a^{-1}-a^{-1}\cdot b\cdot {\left(I+c\cdot a^{-1}\cdot b\right)}^{-1}\cdot c\cdot a^{-1}$

where $I$ stands for a (conformable) identity matrix, and the following Sylvester’s identity [6]

$\text{det}\left(a+b\cdot c\right)=\text{det}\left(a\right)\text{det}\left(I+c\cdot a^{-1}\cdot b\right)$

3. Asymptotic Theory

In this section we investigate the asymptotic distribution of Formula (1) and Formula (3) when $n\to \infty$ , getting the same answer for both of these [7].

3.1. Normal Limit of Formula (2)

To make our task easier, we assume that the sum all probabilities of the Multinomial distribution defined in Formula (2) equals to $r<1$ , and thus allowing for the possibility of an extra outcome whose probability is $1-r$ , denoting its observed total by $W$ . This substantially simplifies subsequent development by removing the singularity of the variance-covariance (V-C) matrix of the $X_{i,j}$ variables, while letting us reach correct conclusions in the $r\to 1$ limit. Correspondingly modifying Formula (2) is easy: extend the bottom row of the multinomial coefficient by $w=n-{\displaystyle {\sum }_{i=1}^K{\displaystyle {\sum }_{j=1}^M{x}_{i,j}}}$ and further multiply by ${\left(1-r\right)}^w$ . Since all variables must add up to $n$ , the resulting $KM+1$ by $KM+1$ V-C matrix is singular, but the singularity disappears when replacing $W$ by $n-{\displaystyle {\sum }_{i=1}^K{\displaystyle {\sum }_{j=1}^M{X}_{i,j}}}$ and considering the distribution of the remaining $X_{i,j}$ ’s only. The corresponding moments are

$\mathbb{E}\left(X_{i.j}\right)=n{p}_i{q}_j$

$\text{Var}\left(X_{i,j}\right)=n{p}_i{q}_j\left(1-p_i{q}_j\right)$

$\text{Cov}\left(X_{i,j},X_{k,m}\right)=-n{p}_i{q}_j{p}_k{q}_m\text{when}i\ne k\text{or}j\ne m$

which remains correct even with the extra $W$ . We then define

$Y_{i,j}:=\frac{X_{i,j}-n{p}_i{q}_j}{\sqrt{n}}$ (4)

and substitute

$x_{i,j}=n{p}_i{q}_j+y_{i,j}\sqrt{n}$

$w=-{\displaystyle \sum }_{i=1}^K{\displaystyle \sum }_{j=1}^M{y}_{i,j}$

into the natural logarithm of Formula (2), after is has been modified by introducing the extra $W$ . Adding ln of the corresponding Jacobian (namely $KM\cdot \ln \sqrt{n}$ , since the distribution of the $Y_{i,j}$ ’s is $KM$ dimensional) and taking the $n\to \infty$ limit of the resulting expression (all limits of this article are routinely done by a computer equipped with Mathematica or a similar software) yields

(5)

This identifies the asymptotic distribution of the $Y_{i,j}$ ’s to be multivariate Normal, with a probability density function (PDF) whose natural logarithm is given by the last expression. Its linear and quadratic moments agree with those of the original discrete distribution, being equal to 0 for the expected value of each $Y_{i,j}$ and to

$p_i{\delta }_{i,k}{q}_j{\delta }_{j,m}-p_i{p}_k{q}_j{q}_m$ (6)

for the covariance between $Y_{i,j}$ and $Y_{k.m}$ (for any combination of indices); the corresponding V-C matrix is thus

$\mathbb{V}=\mathbb{P}\otimes \mathbb{Q}-p\cdot p^{\text{T}}\otimes q\cdot q^{\text{T}}$ (7)

where $p$ and $q$ are column vectors whose elements are the $p_i$ and $q_j$ probabilities (respectively), and $\mathbb{P}$ and $\mathbb{Q}$ are diagonal matrices with the $p_i$ and $q_j$ probabilities (respectively) on the main diagonal. The PDF can then be written, in correspondence with Formula (5), as

$\frac{\text{exp}\left(-\frac{Y^{\text{T}}\circ \mathbb{A}\circ Y}{2}\right)}{{\left(2\pi \right)}^{KM/2}}\sqrt{\det \mathbb{A}}$ (8)

where $Y$ is a column vector whose $KM$ elements are $Y_{i,j}$ (organized in accordance with our Kronecker’s product, i.e. $Y_{1,1},Y_{1,2},\cdots ,Y_{1,M},Y_{2,1},\cdots ,Y_{K,M}$ ),

$\mathbb{A}={\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}+\frac{1_k\cdot 1_k^{\text{T}}\otimes 1_m\cdot 1_m^{\text{T}}}{1-r}={\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}+\frac{1_k\otimes 1_m\circ 1_k^{\text{T}}\otimes 1_m^{\text{T}}}{1-r}$ (9)

and

$\begin{array}{c}\text{det}\left(\mathbb{A}\right)=\text{det}\left({\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}\right)\text{det}\left(1+\frac{1_k^{\text{T}}\otimes 1_m^{\text{T}}\circ \mathbb{P}\otimes \mathbb{Q}\circ 1_k\otimes 1_m}{1-r}\right)\\ =\frac{1}{{\left({{\displaystyle \prod }}_{i=1}^K{p}_i\right)}^M{\left({{\displaystyle \prod }}_{j=1}^M{q}_j\right)}^K\left(1-r\right)}\end{array}$ (10)

the last two implied by Formula (5).

Alternately, they can be derived from Formula (7) using the Woodbury and Silvester identities. Bypassing the details, we only verify that $\mathbb{A}$ is the inverse of $\mathbb{V}$ , by multiplying

$\begin{array}{l}\left({\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}+\frac{1_k\cdot 1_k^{\text{T}}\otimes 1_m\cdot 1_m^{\text{T}}}{1-r}\right)\circ \left(\mathbb{P}\otimes \mathbb{Q}-p\cdot p^{\text{T}}\otimes q\cdot q^{\text{T}}\right)\\ ={\mathbb{I}}_k\otimes {\mathbb{I}}_m-1_k\cdot p^{\text{T}}\otimes 1_m\cdot q^{\text{T}}+\frac{1_k\cdot p^{\text{T}}\otimes 1_m\cdot q^{\text{T}}}{1-r}-\frac{r\left(1_k\cdot p^{\text{T}}\otimes 1_m\cdot q^{\text{T}}\right)}{1-r}\\ ={\mathbb{I}}_k\otimes {\mathbb{I}}_m\end{array}$

where ${\mathbb{I}}_k$ and ${\mathbb{I}}_m$ are $k$ by $k$ and $m$ by $m$ identity matrices respectively.

3.2. Asymptotic Distribution of U

We now take the $n\to \infty$ limit of Formula (1) and of Formula (3), first replacing $X_{i,j}$ with $Y_{i,j}$ by using Formula (4); this results in the same answer in both cases, namely in

${\displaystyle \sum }_{i=1}^K{\displaystyle \sum }_{j=1}^M\frac{Y_{i,j}^2}{p_i{q}_j}-{\displaystyle \sum }_{i=1}^K\frac{Y_{i,•}^2}{p_i}-{\displaystyle \sum }_{j=1}^M\frac{Y_{•,j}^2}{q_i}$ (11)

where $Y_{•,j}={\displaystyle {\sum }_{i=1}^K{Y}_{i,j}}$ and $Y_{i,•}:={\displaystyle {\sum }_{j=1}^M{Y}_{i.j}}$ . To get the corresponding moment generating function (MGF) of $U$ , we need to compute the expected value of $\text{exp}\left(tU\right)$ , i.e.

$\frac{{\displaystyle \int \cdots }{\displaystyle {\int }_{-\infty }^{\infty }\text{exp}\left(-\frac{Y^{\text{T}}\circ \left(\mathbb{A}-2t\mathbb{U}\right)\circ Y}{2}\right)\text{d}Y}}{{\left(2\pi \right)}^{KM/2}}\sqrt{\det \mathbb{A}}$ (12)

where

$\mathbb{U}={\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}-{\mathbb{P}}^{-1}\otimes 1_m\cdot 1_m^{\text{T}}-1_k\cdot 1_k^{\text{T}}\otimes {\mathbb{Q}}^{-1}$ (13)

which is the matrix version of Formula (11); to see this, it helps to re-write Formula (11) as follows

${\displaystyle \sum }_{i=1}^K{\displaystyle \sum }_{k=1}^K{\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{m=1}^M{Y}_{i.j}\left(\frac{{\delta }_{i,k}{\delta }_{j,m}}{p_i{q}_j}-\frac{{\delta }_{i,k}}{p_i}-\frac{{\delta }_{j,m}}{q_j}\right)Y_{k,m}$

To find the result of Formula (12) necessitates computing the determinant of $\mathbb{A}-2t\mathbb{U}$ ; this matrix can be expressed in the following form

$\left(1-2t\right){\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}+\frac{1_k\cdot 1_k^{\text{T}}\otimes 1_m\cdot 1_m^{\text{T}}}{1-r}+2t\left[\begin{array}{ll}{\mathbb{P}}^{-1}\otimes 1_m\hfill & 1_k\otimes {\mathbb{Q}}^{-1}\hfill \end{array}\right]\circ \left[\begin{array}{l}{\mathbb{I}}_k\otimes 1_m^{\text{T}}\hfill \\ 1_k^{\text{T}}\otimes {\mathbb{I}}_m\hfill \end{array}\right]$ (14)

which follows from Formula (9) and Formula (13); the square brackets indicate that the matrix consists of two (and subsequently also of four) blocks.

Pre-multiplying Formula (14) by $\frac{\mathbb{P}\otimes \mathbb{Q}}{1-2t}$ , whose determinant is

$D_1:=\frac{{\left({{\displaystyle \prod }}_{i=1}^K{p}_i\right)}^M{\left({{\displaystyle \prod }}_{j=1}^M{q}_j\right)}^K}{{\left(1-2t\right)}^{KM}}$

leaves us with the task of finding the determinant of

$\begin{array}{c}ℂ:={\mathbb{I}}_k\otimes {\mathbb{I}}_m+\frac{p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}}{\left(1-r\right)\left(1-2t\right)}+\frac{2t}{1-2t}\left[\begin{array}{ll}{\mathbb{I}}_k\otimes q\hfill & p\otimes {\mathbb{I}}_m\hfill \end{array}\right]\circ \left[\begin{array}{l}{\mathbb{I}}_k\otimes 1_m^{\text{T}}\hfill \\ 1_k^{\text{T}}\otimes {\mathbb{I}}_m\hfill \end{array}\right]\\ :=\mathbb{F}+\frac{2t}{1-2t}\mathbb{G}\circ ℍ\end{array}$ (15)

This needs to be done in two steps: for the determinant of $\mathbb{F}$ , we get

$D_2:=\det \mathbb{F}=1+\frac{1_k^{\text{T}}\otimes 1_m^{\text{T}}\circ p\otimes q}{\left(1-r\right)\left(1-2t\right)}=1+\frac{r}{\left(1-r\right)\left(1-2t\right)}=\frac{1-2t\left(1-r\right)}{\left(1-r\right)\left(1-2t\right)}$

while its inverse equals to

${\mathbb{F}}^{-1}={\mathbb{I}}_k\otimes {\mathbb{I}}_m-\frac{p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}}{1-2t\left(1-r\right)}$ (16)

easily verifiable by simple multiplication, utilizing $p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}\circ p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}=r\left(p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}\right)$ .

Sylvester’s identity implies that the determinant of $ℂ$ is a product of $\det \mathbb{F}$ and the determinant of (at this point, we can start replacing $r$ by 1)

$\begin{array}{l}{\mathbb{I}}_{k+m}+\frac{2t}{1-2t}\left[\begin{array}{l}{\mathbb{I}}_k\otimes 1_m^{\text{T}}\hfill \\ 1_k^{\text{T}}\otimes {\mathbb{I}}_m\hfill \end{array}\right]\circ \left({\mathbb{I}}_k\otimes {\mathbb{I}}_m-p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}\right)\circ \left[\begin{array}{ll}{\mathbb{I}}_k\otimes q\hfill & p\otimes {\mathbb{I}}_m\hfill \end{array}\right]\\ ={\mathbb{I}}_{k+m}+\frac{2t}{1-2t}\left[\begin{array}{cc}{\mathbb{I}}_k\otimes 1& p\otimes 1_m^{\text{T}}\\ 1_k^{\text{T}}\otimes q& 1\otimes {\mathbb{I}}_m\end{array}\right]-\frac{2t}{1-2t}\left[\begin{array}{cc}p\cdot 1_k^{\text{T}}\otimes 1& p\otimes 1_m^{\text{T}}\\ 1_k^{\text{T}}\otimes q& 1\otimes q\cdot 1_m^{\text{T}}\end{array}\right]\\ =\frac{1}{1-2t}\left[\begin{array}{cc}\left({\mathbb{I}}_k-2tp\cdot 1_k^{\text{T}}\right)\otimes 1& 0_k\otimes 0_m^{\text{T}}\\ 0_m\otimes 0_k^{\text{T}}& 1\otimes \left({\mathbb{I}}_m-2tq\cdot 1_m^{\text{T}}\right)\end{array}\right]\end{array}$ (17)

where $0$ is a zero column-vector of indicated length and 1 is a scalar. The determinant of the last matrix is

$D_3:=\frac{{\left(1-2t\right)}^2}{{\left(1-2t\right)}^{K+M}}$

since the determinant of each main-diagonal block is equal to $1-2t,$ due to Sylvester’s identity.

In the $r\to 1$ limit Formula (12), whose value is given by $\sqrt{\frac{\det \mathbb{A}}{\det \left(\mathbb{A}-2t\mathbb{U}\right)}}=\sqrt{\frac{\det \mathbb{A}}{D_1{D}_2{D}_3}}$ , then equals to

${\left(\frac{1}{1-2t}\right)}^{\left(K-1\right)\left(M-1\right)/2}$ (18)

proving the well-known result that, to this level of approximation, the distribution of $U$ is ${\chi }^2$ with $\left(K-1\right)\left(M-1\right)$ degrees of freedom.

3.3. Inverting $\mathbb{A}-2t\mathbb{U}$

Making the integrand of Formula (12) into a PDF of a Normal distribution results in

$\frac{\text{exp}\left(-\frac{Y^{\text{T}}\circ \left(\mathbb{A}-2t\mathbb{U}\right)\circ Y}{2}\right)}{{\left(2\pi \right)}^{KM/2}}\sqrt{\det \mathbb{A}-2t\mathbb{U}}$ (19)

whose higher (quartic and hexic in particular) moments are key to finding $\frac{1}{n}$ -proportional corrections to the MGF of each Formula (1) and Formula (3). Furthermore, these moments are polynomial functions of its second moments (a special property of a multivariate Normal distribution), which in turn are simply the elements of ${\left(\mathbb{A}-2t\mathbb{U}\right)}^{-1}$ .

To compute this inverse, we first write Formula (14) as a product of $\left(1-2t\right){\mathbb{P}}^{-1}\otimes {\mathbb{Q}}^{-1}$ (whose inverse is easy) and of $ℂ=\mathbb{F}+\frac{2t}{1-2t}\mathbb{G}\circ ℍ$ , introduced in Formula (15). From Formula (16) we know that

${\mathbb{F}}^{-1}={\mathbb{I}}_k\otimes {\mathbb{I}}_m-p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}$

as $r$ is no longer needed and can be set equal to 1.

To complete the exercise, we need to simplify

${\mathbb{F}}^{-1}-\frac{2t}{1-2t}{\mathbb{F}}^{-1}\circ \mathbb{G}\circ {\left({\mathbb{I}}_{k+m}+\frac{2t}{1-2t}ℍ\circ {\mathbb{F}}^{-1}\circ \mathbb{G}\right)}^{-1}\circ ℍ\circ {\mathbb{F}}^{-1}$

obtained from Woodbury’s identity as a part of inverting $ℂ$ . The last line of Formula (17) has already done exactly that for the matrix in parentheses; the corresponding inverse is

$\left(1-2t\right)\left[\begin{array}{cc}{\mathbb{I}}_k\otimes 1& 0_k\otimes 0_m^{\text{T}}\\ 0_m\otimes 0_k^{\text{T}}& 1\otimes {\mathbb{I}}_m\end{array}\right]+2t\left[\begin{array}{cc}p\cdot 1_k^{\text{T}}\otimes 1& 0_k\otimes 0_m^{\text{T}}\\ 0_m\otimes 0_k^{\text{T}}& 1\otimes q\cdot 1_m^{\text{T}}\end{array}\right]$

easily verifiable by simple matrix multiplication (note that both matrices in the last expression are idempotent).

This inverse still needs to be pre-multiplied by $\left[\begin{array}{ll}{\mathbb{I}}_k\otimes q\hfill & p\otimes {\mathbb{I}}_m\hfill \end{array}\right]$ and post-multiplied by $\left[\begin{array}{l}{\mathbb{I}}_k\otimes 1_m^{\text{T}}\hfill \\ 1_k^{\text{T}}\otimes {\mathbb{I}}_m\hfill \end{array}\right]$ , thus getting

$\left(1-2t\right)\left({\mathbb{I}}_k\otimes q\cdot 1_m^{\text{T}}+p\cdot 1_k^{\text{T}}\otimes {\mathbb{I}}_m\right)+4tp\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}$

and further pre and post-multiplied by ${\mathbb{F}}^{-1}$ , resulting in

$\left(1-2t\right)\left({\mathbb{I}}_k\otimes q\cdot 1_m^{\text{T}}+p\cdot 1_k^{\text{T}}\otimes {\mathbb{I}}_m-2p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}\right)+4t{0}_k\cdot 0_k^{\text{T}}\otimes 0_m\cdot 0_m^{\text{T}}$

Multiplying the last expression by $\frac{2t}{1-2t}$ and subtracting the result from ${\mathbb{F}}^{-1}$ yields ${ℂ}^{-1}$ , namely

${\mathbb{I}}_k\otimes {\mathbb{I}}_m+\left(4t-1\right)p\cdot 1_k^{\text{T}}\otimes q\cdot 1_m^{\text{T}}-2t\left({\mathbb{I}}_k\otimes q\cdot 1_m^{\text{T}}+p\cdot 1_k^{\text{T}}\otimes {\mathbb{I}}_m\right)$

Finally, post-multiplying by $\frac{\mathbb{P}\otimes \mathbb{Q}}{1-2t}$ results in the following expression for ${\left(\mathbb{A}-2t\mathbb{U}\right)}^{-1}$

$\frac{\mathbb{P}\otimes \mathbb{Q}}{1-2t}+\left(\frac{1}{1-2t}-2\right)p\cdot p^{\text{T}}\otimes q\cdot q^{\text{T}}+\left(1-\frac{1}{1-2t}\right)\left(\mathbb{P}\otimes q\cdot q^{\text{T}}+p\cdot p^{\text{T}}\otimes \mathbb{Q}\right)$

and the corresponding formula for computing the desired second moments, namely

$\begin{array}{c}\mathbb{E}\left(Y_{i,j}{Y}_{k,m}\right)=-2p_i{p}_k{q}_j{q}_m+{\delta }_{i.k}{p}_i{q}_j{q}_m+{\delta }_{j,m}{p}_i{p}_k{q}_j\\ +\frac{p_i\left({\delta }_{i.k}-p_k\right)q_j\left({\delta }_{j,m}-q_m\right)}{1-2t}\end{array}$ (20)

(note the expected agreement with Formula (6) when $t=0$ ).

Consequently,

$\begin{array}{l}\mathbb{E}\left(Y_{i,j}{Y}_{k,•}\right)=p_i{q}_j\left({\delta }_{i,k}-p_k\right)\\ \mathbb{E}\left(Y_{i,j}{Y}_{•,m}\right)=p_i{q}_j\left({\delta }_{j,m}-q_m\right)\\ \mathbb{E}\left(Y_{i,•}{Y}_{k,•}\right)=p_i\left({\delta }_{i.k}-p_k\right)\\ \mathbb{E}\left(Y_{•,j}{Y}_{•,m}\right)=q_j\left({\delta }_{j,m}-q_m\right)\\ \mathbb{E}\left(Y_{i,•}{Y}_{•,j}\right)=0\end{array}$ (21)

provide the remaining moments needed to proceed with computing the desired $\frac{1}{n}$ -proportional corrections. Note that fourth-order moments are then found from

$\mathbb{E}\left(abcd\right)=\mathbb{E}\left(ab\right)\mathbb{E}\left(cd\right)+\mathbb{E}\left(ac\right)\mathbb{E}\left(bd\right)+\mathbb{E}\left(ad\right)\mathbb{E}\left(bc\right)$

(with subsequent special cases, such as $\mathbb{E}\left(a^2{b}^2\right)=2\mathbb{E}{\left(ab\right)}^2+\mathbb{E}\left(a^2\right)\mathbb{E}\left(b^2\right)$ etc.), and the sixth-order moments from

$\mathbb{E}\left(abcdef\right)=\mathbb{E}\left(ab\right)\mathbb{E}\left(cd\right)\mathbb{E}\left(ef\right)+\cdots +\mathbb{E}\left(af\right)\mathbb{E}\left(cd\right)\mathbb{E}\left(be\right)$

(where the RHS consists of 15 terms corresponding to all possible ways of pairing the six arguments), usually needed in its special form such as

$\mathbb{E}\left(a^3{b}^3\right)=9\mathbb{E}\left(a^2\right)\mathbb{E}\left(b^2\right)\mathbb{E}\left(ab\right)+6\mathbb{E}{\left(ab\right)}^3$

etc.

4. $\frac{1}{n}$ -Proportional Corrections

To find the desired correction to the ${\chi }^2$ distribution of the Pearson test, we must

include $\frac{1}{\sqrt{n}}$ and $\frac{1}{n}$ -proportional terms in the Formula (13) expansion, thus getting (listing the extra terms only)

$\begin{array}{l}\mathfrak{U}=\cdots +\frac{{{\displaystyle \sum }}_{i=1}^K\frac{Y_{i,•}^3}{p_i^2}+{{\displaystyle \sum }}_{j=1}^M\frac{Y_{•,j}^3}{q_j^2}}{\sqrt{n}}+\frac{{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\left(2\frac{Y_{i,j}{Y}_{i,•}{Y}_{•,j}}{p_i{q}_j}-\frac{Y_{i,j}^2{Y}_{i,•}}{p_i^2{q}_j}-\frac{Y_{i,j}^2{Y}_{•,j}}{p_i{q}_j^2}\right)}{\sqrt{n}}\\ +\frac{{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\left(-2\frac{Y_{i,j}{Y}_{i,•}^2{Y}_{•.j}}{p_i^2{q}_j}-2\frac{Y_{i,j}{Y}_{i,•}{Y}_{•.j}^2}{p_i{q}_j^2}+\frac{Y_{i,j}^2{Y}_{i,•}^2}{p_i^3{q}_j}+\frac{Y_{i,j}^2{Y}_{•.j}^2}{p_i{q}_j^3}+\frac{Y_{i,j}^2{Y}_{i,•}{Y}_{•.j}}{p_i^2{q}_j^2}+\frac{Y_{i,•}^2{Y}_{•.j}^2}{p_i{q}_j}\right)}{n}\\ +\frac{-{{\displaystyle \sum }}_{i=1}^K\frac{Y_{i,•}^4}{p_i^3}-{{\displaystyle \sum }}_{j=1}^M\frac{Y_{•.j}^4}{q_j^3}}{n}\\ :=\frac{u_1}{\sqrt{n}}+\frac{u_2}{\sqrt{n}}+\frac{u_3}{n}+\frac{u_4}{n}\end{array}$

Similarly, we need the extra terms in the expansion of the natural logarithm of the RHS of Formula (2); this time we get

$\begin{array}{c}\mathcal{A}=\cdots +\frac{{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\left(\frac{Y_{i,j}^3}{6p_i^2{q}_j^2}-\frac{Y_{i,j}}{2p_i{q}_j}\right)}{\sqrt{n}}+\frac{1}{12n}\\ -\frac{{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\left(\frac{Y_{i,j}^4}{12p_i^3{q}_j^3}-\frac{Y_{i,j}^2}{4p_i^2{q}_j^2}+\frac{1}{12p_i{q}_j}\right)}{n}\\ :=\frac{a_1}{\sqrt{n}}+\frac{1}{12n}+\frac{a_2}{n}\end{array}$

Combining these two, we expand $\text{exp}\left(\mathcal{A}-2t\mathfrak{U}\right)$ up to and including $\frac{1}{n}$ -proportional terms, getting

$\frac{2\left(u_1^2+u_2^2+2u_1{u}_2\right)t^2-2t\left(u_3+u_4\right)-2t\left(u_1+u_2\right)a_1+\frac{1}{2}{a}_1^2+a_2+\frac{1}{12}}{n}$ (22)

(having discarded $\frac{1}{\sqrt{n}}$ -proportional terms, whose expected values are equal to zero).

What remains to be done is to take the expected value of Formula (22); this is done by first expanding every sum in this expression, then finding the expected value of each term of these expansions, finally followed by carrying out the summation itself (over up to four indices). This results in scores of individual contributions, even though, rather surprisingly, most of them cancelling each other upon subsequent simplification. This computation can be successfully completed only with a computer’s help; the actual Mathematica program required many lines of code and is available upon request.

As an example, we demonstrate how to find $\mathbb{E}\left(u_1{u}_2\right)$ . We start by ignoring the triple summations, multiplying only the first term of $u_1$ by $u_2$ (note that this requires using a different index for $u_1$ ) and expanding the resulting expected value, thus

$\begin{array}{l}\mathbb{E}\left(\frac{Y_{k,•}^3}{p_k^2}\left(2\frac{Y_{i,j}{Y}_{i,•}{Y}_{•,j}}{p_i{q}_j}-\frac{Y_{i,j}^2{Y}_{i,•}}{p_i^2{q}_j}-\frac{Y_{i,j}^2{Y}_{•,j}}{p_i{q}_j^2}\right)\right)\\ =\frac{6\mathbb{E}\left(Y_{i,•}{Y}_{k,•}\right)\mathbb{E}{\left(Y_{i,j}{Y}_{k,•}\right)}^2-3\mathbb{E}\left(Y_{i,•}{Y}_{k,•}\right)\mathbb{E}\left(Y_{k,•}^2\right)\mathbb{E}\left(Y_{i,j}^2\right)-6\mathbb{E}\left(Y_{i,j}{Y}_{i,•}\right)\mathbb{E}\left(Y_{i,j}{Y}_{k,•}\right)\mathbb{E}\left(Y_{k,•}^2\right)}{p_i^2{q}_j{p}_k^2}\\ -\frac{6\mathbb{E}\left(Y_{i,j}{Y}_{•,j}\right)\mathbb{E}\left(Y_{i,j}{Y}_{k,•}\right)\mathbb{E}\left(Y_{k,•}^2\right)}{p_i{q}_j^2{p}_k^2}+\frac{6\mathbb{E}\left(Y_{i,j}{Y}_{•,j}\right)\mathbb{E}\left(Y_{i,•}{Y}_{k,•}\right)\mathbb{E}\left(Y_{k,•}^2\right)}{p_i{q}_j{p}_k^2}\end{array}$

where we have already discarded terms containing a zero factor such as $\mathbb{E}\left(Y_{i,•}{Y}_{•,j}\right)$ . We then evaluate the expected values based on Formula (20) and Formula (21), and complete the $i,j$ and $k$ summations, resulting in

$-12+18K+9K^2-15P_0+\frac{3\left(M-1\right)\left(K^2-P_0\right)}{1-2t}$

where $P_0:={\displaystyle {\sum }_{i=1}^K{p}_i^{-1}}$ . To include the second term of $u_1$ , we simply interchange $K$ and $M$ of the previous answer, replace $P_0$ by $Q_0:={\displaystyle {\sum }_{j=1}^M{q}_j^{-1}}$ , and add the result to the previous answer.

When this is completed for all terms of Formula (22), adding them and simplifying (many terms cancel out, as mentioned already) results in

$\begin{array}{l}\frac{d}{2n}\left(\frac{1}{1-2t}-1\right)\\ +\frac{\tilde{P}\tilde{Q}-2\left(M-1\right)\tilde{P}-2\left(K-1\right)\tilde{Q}-2\left(K+M-3\right)d}{8n}\left(\frac{1}{{\left(1-2t\right)}^2}-\frac{2}{1-2t}+1\right)\\ +\frac{5\tilde{P}\tilde{Q}+2\left(M-1\right)\left(M-2\right)\tilde{P}+2\left(K-1\right)\left(K-2\right)\tilde{Q}+2d\left(K-2\right)\left(M-2\right)}{24n}\\ \times \left(\frac{1}{{\left(1-2t\right)}^3}-\frac{3}{{\left(1-2t\right)}^2}+\frac{3}{1-2t}-1\right)\end{array}$ (23)

where $\tilde{P}:=P_0-K^2,\tilde{Q}:=Q_0-M^2$ and $d:=\left(K-1\right)\left(M-1\right)$ ; note that when $p_i=\frac{1}{K}$ for all $i$ values, $\tilde{P}$ is then equal to zero, and similarly $q_j=\frac{1}{M}$ implies $\tilde{Q}=0$ .

The last expression, further divided by ${\left(1-2t\right)}^d$ , is the actual correction to the MGF of the test statistic Formula (13), implying that the corresponding PDF correction is a linear combination of ${\chi }^2$ distributions; to get an explicit formula for this correction, make the following replacement in Formula (23)

$\left(\frac{1}{1-2t}-1\right)\to \left(\frac{x}{d}-1\right){\chi }_d^2\left(x\right)$

$\left(\frac{1}{{\left(1-2t\right)}^2}-\frac{2}{1-2t}+1\right)\to \left(\frac{x^2}{d\left(d+2\right)}-\frac{2x}{d}+1\right){\chi }_d^2\left(x\right)$

$\left(\frac{1}{{\left(1-2t\right)}^3}-\frac{3}{{\left(1-2t\right)}^2}+\frac{3}{1-2t}-1\right)\to \left(\frac{x^3}{d\left(d+2\right)\left(d+4\right)}-\frac{3x^2}{d\left(d+2\right)}+\frac{3x}{d}-1\right){\chi }_d^2\left(x\right)$

where

${\chi }_d^2\left(x\right)=\frac{x^{d/2-1}\exp \left(-x/2\right)}{2^{d/2}\Gamma \left(d/2\right)}$

is the PDF of the ${\chi }_d^2$ distribution. Note that each of the three corrections integrates to zero (the main reason for keeping them in this form).

These results are in perfect agreement with those obtained, using a different technique, by [8].

G-Test Results

This time we get the following extra terms in the expansion of Formula (3)

$\begin{array}{c}\mathfrak{U}=\cdots +\frac{2{{\displaystyle \sum }}_{i=1}^K\frac{Y_{i,•}^3}{p_i^2}+2{{\displaystyle \sum }}_{i=1}^M\frac{Y_{•,j}^3}{q_j^2}}{3\sqrt{n}}+\frac{-2{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\frac{Y_{i,j}^3}{p_i^2{q}_j^2}}{3\sqrt{n}}\\ +\frac{{{\displaystyle \sum }}_{i=1}^K\frac{Y_{i,•}^4}{p_i^3}+{{\displaystyle \sum }}_{i=1}^M\frac{Y_{•,j}^4}{q_j^3}}{3n}+\frac{-{{\displaystyle \sum }}_{i=1}^K{{\displaystyle \sum }}_{j=1}^M\frac{Y_{i,j}^4}{p_i^3{q}_j^3}}{3n}\\ :=\frac{u_1}{\sqrt{n}}+\frac{u_2}{\sqrt{n}}+\frac{u_3}{n}+\frac{u_4}{n}\end{array}$

We then repeat the steps of the previous section, getting the following (surprisingly simple) correction to the corresponding PDF

$\frac{\tilde{P}\tilde{Q}+\left(M^2-1\right)\tilde{P}+\left(K^2-1\right)\tilde{Q}+\left(K^2-1\right)\left(M^2-1\right)}{12n}\left(\frac{x}{d}-1\right){\chi }_d^2\left(x\right)$

5. Monte-Carlo Verification and Conclusion

To demonstrate the utility of our correction formulas (but also their limitations, in case of the G-test), we generate a million observations from a contingency-table model and compare the results with the ${\chi }^2$ approximation, both without and including the $\frac{1}{n}$ -proportional correction.

Our first example assumes that two attributes have 10 equally likely levels each (for a total of 100 cells), and the number of observations is equal to 150. The results are displayed in Figure 1 for the Pearson test and in Figure 2 for the G-test (the dotted lines show the basic ${\chi }^2$ approximation, the solid lines include the respective corrections). From these graphs it is quite obvious that, in the former case, the correction terms are quite capable of removing the (still rather serious) error, while the G-test is so hopelessly inaccurate that repairing it becomes impossible (the coefficient of our formula has a value bigger than 5, being no longer ‘small’, resulting in making the ‘corrected’ PDF partly negative; a strong indication that the error is too large to be dealt with in this manner). The benefit of our correction is to clearly indicate a presence of an inacceptable error in the usual ${\chi }^2$ approximation (suggesting a solution would require a separate study). We conclude that the G-test should not be used with large number of cells (unless the number of observations is correspondingly large).

Figure 1. Pearson’s test correction (n = 150).

Figure 2. G-test and its error (n = 150).

Nevertheless, we can still confirm that our G-test correction works when the attributes have only 3 levels each (we will still make them equally likely) and $n$ is equal to 35; this is indicated in Figure 3.

Figure 3. G-test corrected (n = 35).

Note that the empirical distribution (represented by a histogram) has some minor irregularities due to the exact distribution’s discrete nature (correcting for these would require a new study). Nevertheless, the graph indicates that the right-tail region is not affected by this phenomenon (this is particularly true for the Pearson test, as seen in Figure 4). We should also mention that, in this same situation and using the same data, the Pearson test is still more accurate than the G-test, both when using the basic ${\chi }^2$ approximation and using its corrected version. In conclusion, we find no reason for recommending the G-test over the Pearson test of independence.

Figure 4. Correcting Pearson’s test (n = 35).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Kendall, M.G. and Stuart, A. (1973) The Advanced Theory of Statistics, Vol. 2, Chapter 33. Harper Publishing Company.
[2]	MacDonald, H.J. (2014) G-Test of Goodness-of-Fit, Handbook of Biological Statistics. 3rd Edition, Sparky House Publishing, 53-58.
[3]	Sokal, R.R. and Rohlf, F.J. (1981) Biometry: The Principles and Practice of Statistics in Biological Research. 2nd Edition, Freeman.
[4]	Wikipedia. Kronecker Product.
[5]	Wikipedia. Woodbury Matrix Identity.
[6]	Wikipedia. Sylvester’s Determinant Identity.
[7]	Serfling, R.J. (1980). Approximation Theorems of Mathematical Statistics. Wiley. https://doi.org/10.1002/9780470316481
[8]	Vrbik, J. (2014) Improving Accuracy of Test of Independence. Advances and Applications in Statistics, 39, 91-94.