Doing a chi-square independence test in Mathematica

Question

I have a matrix consisting of the following data:

Needs["HypothesisTesting`"]
data = {{11, 206}, {32, 1374}}

How would I make a chi-square independence test with a level of 5% of significance out of this? When I try the following:

In[9]:= ChiSquarePValue[data, 0.05]

I get this a very little value out of it.

Out[9]= OneSidedPValue -> 1.89247*10^-303

However, the correct answer would be something around 0.017 or 1.7%.

Answer 1

When conducting a Chi-square test on a two-way table, you want to create and inspect the following. (The calculations are made in a way that generalizes to any $r$ by $c$ table.)

1. The data:

   data = {{11, 206}, {32, 1374}}; 
   rc =  {{"Row 1", "Row 2"}, {"Column 1", "Column 2"}};
   TableForm[data, TableHeadings -> rc]
   

2. The fit. This is obtained from the row and column sums:

   fit = Outer[Times, Plus @@@ data, Plus @@@ Transpose[data]] / Plus @@ Flatten[data];
   TableForm[fit // N, TableHeadings -> rc]
   

3. The residuals, equal to the differences between the data and the fit:

   residual = data - fit;
   TableForm[residual // N, TableHeadings -> rc]
   

   (For a 2 by 2 table, all residuals will have equal size.)

4. The squared residuals, scaled by the reciprocal of the fit. Where these are substantially larger than $1$ in absolute value, they signal bad fits:

   χ2array = residual^2 / fit;
   TableForm[χ2array // N, TableHeadings -> rc]
   

   In this example, the entry for row 1, column 1 has a value of 4.8, suggesting a (slightly) bad fit there. The other entries are all small, indicating decent to excellent fits. (Appropriately signed square roots of these values are normally considered "residuals", but it's not really necessary to do this extra computation.)

5. The sum of these scaled squared residuals. This is the chi-squared statistic, $\chi^2$. It is an overall measure of how well the fit matches the data.

   χ2 = Plus @@ Flatten[χ2array];
   χ2 // N
   

   Here, it equals 5.68632.

6. The p-value. This assesses the chance that a chi-squared random variable could attain a value of $\chi^2$ or larger. As a preliminary step, we need to compute the "degrees of freedom" (df) of the statistic.

   df = Length[Flatten[data] ] - Length[data] - Length[Transpose[data]] + 1;
   1 - CDF[ChiSquareDistribution[df], χ2] // N
   

   This calculation returns 0.0170977, or about 1.7%, based on one degree of freedom.

In a full report of the test, all these results would be presented. In an abbreviated report, only df, $\chi^2$, and the p-value would be given (as computed in steps 5 and 6).

Finally, to conduct the test at the 5% level, one would remark that the p-value is less than 5%. Because you have generated these intermediate results and inspected the residual tables (steps 3 and 4), you might remark that this low p-value appears to be due solely to a lack of fit in the first row and column. You might urge some caution in the interpretation because (looking at the data and fit tables, steps 1 and 2) you notice the counts in this cell (11 and 5.75) are small. In fact, you might elect to confirm your result with a permutation test or, when applicable, Fisher's Exact Test.

As a double-check--because these calculations have been coded from scratch--you might compare the results to a calculation with other software, such as R; e.g.,

chisq.test(matrix(c(11,32,206,1374), nrow=2), correct=FALSE)

This produces the abbreviated summary:

X-squared = 5.6863, df = 1, p-value = 0.01710

By saving the result of this calculation you can inspect the auxiliary material and compare it with the Mathematica calculations. Unlike R (or almost any other statistical program), Mathematica will present exact results: simply remove the "// N" bits from the code. (It can be surprising how many people have gotten into trouble by over-rounding intermediate results in their statistical calculations; using exact arithmetic avoids that problem.)

Answer 2

Or use Poisson regression via GeneralizedLinearModelFit. Use dummy variables to indicate rows and columns.

data = {{1, 1, 11}, {1, 0, 206}, {0, 1, 32}, {0, 0, 1374}};
glm = GeneralizedLinearModelFit[data, {x1, x2}, {x1, x2},ExponentialFamily -> "Poisson"];
testStatistic = glm["PearsonChiSquare"]

This gives you the test statistic. Then get the P-value:

df = glm["ResidualDegreesOfFreedom"];
Needs["HypothesisTesting`"];
ChiSquarePValue[testStatistic, df]

You can also get the expected values:

glm["PredictedResponse"]

Much more is in the documentation for GeneralizedLinearModelFit.

Answer 3

As the manual clearly indicates, ChiSquarePValue expects as parameters x and df, x being a cumulative probability beyond x and df a number indicating the degrees of freedom. Your parameters are a matrix and an alpha level. Obviously not what is required.

In this case you have to calculate the $\chi^2$ value from the table. This Wikipedia lemma shows you how to do that. The degrees of freedom is the number of columns times the number of rows in the table both minus 1 (result here: df=1).

I don't have Mathematica at hand, but a manual calculation gets me a $\chi^2$ value of 5.68, which according to Wolfram|Alpha, corresponds to a one-sided p-value of 0.017.

Answer 4

Summary of whuber's answer with J. M's comment with another example data of 3x2 matrix about

data = {{58, 35}, {11, 25}, {10, 23}}; 
rc = {{"agree", "neutral", "disgaree"}, {"men", "women"}};
TableForm[data, TableHeadings -> rc]
fit = Outer[Times, Total[data, {2}], Total[data]]/Total[data, 2];
TableForm[fit // N, TableHeadings -> rc]
residual = data - fit;
TableForm[residual // N, TableHeadings -> rc]
\[Chi]2array = residual^2/fit;
TableForm[\[Chi]2array // N, TableHeadings -> rc]
\[Chi]2 = Total[\[Chi]2array, 2];
\[Chi]2 // N
df = (Times @@ # - Total[#] + 1) &[Dimensions[data]]; 
1 - CDF[ChiSquareDistribution[df], \[Chi]2] // N