# $\chi^2$--test - Comparing a categorical variable across the values of another categorical variable {#chi-sq}

Table (dataChiSq) The dataset we consider in this chapter. The two variables are both categorical.

| Martian | Nose colour | Antenna |
|---------|-------------|---------|
| 1       | blue        | yes     | 
| 2       | red         | yes     | 
| 120     | blue        | no      |

Here we are interested in the relationship between two categorical variables, e.g. Martian nose-color (red or blue) and the presence or absence of of antennae. The question we are interested in here is whether the presence of antennae somehow has an effect on the nose colour, i.e. whether these two variables are associated. From Figure \@ref(fig:twogrouptechiques) we know that the test appropriate for this situation is  the $\chi^2$ test (say "kai-squared").
                
As we will learn later (see \@ref(sec:fishers)),  our sample of 20 Martians from sites I + II is too small for this type of analysis, so for the purpose of this chapter we'll imagine that we have made a second trip to Mars and have observed the outcome of the two variables on 120 Martians. Table \@ref(tab:dataChiSq) shows the resulting data table. Arranged in a different form, the data can be looked at in a more convenient way as in shown in Table \@ref(contigencytab). 

```{r echo = FALSE, eval = TRUE}
martian_chi <- read.delim("data/martianchi.txt")
martian <- read.delim("data/martian.txt")
```

contigencytab: Observed contingency table for Martian colour and presence/absence of antennae in our entire Martian sample - 120 observations in total, taken from both sites.

_Observed table of counts_ XX

<!-- |             | Antennae |      |       | -->
| Nose Colour | No       | Yes  | Total |
|-------------|----------|------|-------|
| Red         | $45$     | $25$ | $70$  |
| Blue        | $22$     | $28$ | $50$  |
| Total       | $67$     | $53$ | $120$ |

This sort of table is called a _2x2--contingency table_. Examining it gives you a first hint whether the two variables are associated. However, obviously we are not done yet, so as a next step we rewrite the contingency table in terms of proportions of antennae  (yes or no) in the two nose colour categories:

_Observed table of row proportions_ To get this table, we divided each row in Table \ref{contigencytab} by the corresponding row total, i.e. 70 red--nosed, 50 blue--nosed Martians and 120 Martians in total.

<!-- |             | Antennae  |           |            | -->
| Nose Colour | No        | Yes       | Total      |
|-------------|-----------|-----------|------------|
| Red         | $64.29\%$ | $35.71\%$ | $100.00\%$ |
| Blue        | $44.00\%$ | $55.00\%$ | $100.00\%$ |
| Total       | $55.83\%$ | $44.17\%$ | $100.00\%$ |

It seems that having _no_ antennae is more common in the red--nosed Martians than in the blue--nosed ones. Is this observation a genuine pattern or just due to chance?! 

This is what we are trying to find out with a $\chi^2$--test. The corresponding null hypothesis can be written as $H_0:$ "there is no association between the two variables" ^[Note that once again the null hypothesis corresponds to the state of boredeom.]. To test whether we can reject $H_0$, we need a test statistic and a null distribution. As usual the test statistic is of the form

$$\frac{\text{signal}}{\text{noise}} = \frac{\text{(observed value - hypothesised value)} }{\text{scale factor}}.$$

As you can guess, the "observed value" corresponds our original counts in the contingency table. So here, we have four observed values! As a next step, we need a way to derive the respective hypothesised values. 

Luckily, the hypothesised values are quite easy to calculate: if there is no association between the two variables, then the nose colour doesn't matter when we look at the proportions of Martians with/without antennae. So our best guess are the proportions that we calculated based on the totals in Table \@ref(tab:ChiPropor), i.e. 55.83% no antennae, 44.17% antennae. Under the premise of $H_0$, these proportions would also hold for the red and the blue--nosed Martians, i.e. out of the 70 red--nosed Martians, 55.83% shouldn't have antennae and 44.17% should. The same consideration applies for the 50 blue--nosed Martians. Table \@ref(tab:ChiCalc) shows the corresponding calculations.

_Calculating expected values_ For example, the expected number in the first row and the first column is calculated by multiplying the expected proportion of Martians without antennae, i.e. 55.83%, with the the known number of red--nosed Martians, i.e. 70. The 55.83% corresponds to the overall proportion of Martians without antennae, i.e. 67/120.

<!-- |             | Antennae                 |                          |       | -->
| Nose Colour | No                       | Yes                      | Total |
|-------------|--------------------------|--------------------------|-------|
| Red         | $\frac{70\times67}{120}$ | $\frac{70\times53}{120}$ | $70$  |
| Blue        | $\frac{50\times67}{120}$ | $\frac{50\times53}{120}$ | $50$  |
| Total       | $67$                     | $53$                     | $120$ |

As a result, we get the following cell numbers which quantify what we'd expect when our null hypothesis of no association between colour and having antennae were true:

_Expected table of counts_ The _expected_ contingency table for Martian nose colour and presence/absence of antennae. The values are calculated as in Table \@ref(tab:ChiCalc)

<!-- |             | Antennae |          |       | -->
| Nose Colour | No       | Yes      | Total |
|-------------|----------|----------|-------|
| Red         | $39$ | $31$ | $70$ |
| Blue        | $28$ | $22$ | $50$ |
| Total       | $67$ | $53$ | $120$ |

In order to arrive at a formula that is generally applicable to any data situation, let's rewrite the steps that we have taken so far in a general form. We start with a general contingency table of observed counts in which we give each value a name according to its position in the table, i.e. $n_{11}$, $n_{12}$, $n_{21}$ and $n_{22}$ cf. Table \@ref(tab:ChiGeneral). Also let $e_{11}$, $e_{12}$, $e_{21}$ and $e_{22}$ be the corresponding expected values. 


_General observed table of counts_  The two indices in the cells represent the row and the column that we are looking at, e.e. $n_{11}$ represents the count in the first row and the first column. $n_{1r}$ denotes the total count in row 1, i.e. $n_{1r} = n_{11} + n_{12}$. $n_{c1}$ denotes the total count in column 1, i.e. $n_{c1} = n_{11} + n_{21}$. As before, $n$ is the total number which is the sum, i.e. $n = n_{11} + n_{12} + n_{21} + n_{22}$.

<!-- |             | Antennae |          |       | -->
| Nose Colour | No       | Yes      | Total |
|-------------|----------|----------|-------|
| Red         | $n_{11}$ | $n_{12}$ | $n_{1r}$ |
| Blue        | $n_{21}$ | $n_{22}$ | $n_{2r}$ |
| Total       | $n_{c1}$ | $n_{c2}$ | $n$ |

With this notation we can derive the following components of the test statistic



**observed value:**
The observed cell count for row $i$ and column $j$ is $n_{ij}$, e.g. $n_{11}=45$.

**hypothesised value:**
The expected cell count under the null hypothesis of no association is $e_{ij} = n_{ir}\times n_{cj}/n$, i.e.  $e_{ij}= \frac{70\times67}{120}=39$.

**scale factor:**
For reasons that we won't get into, the scale factor is the same as the expected count for each cell.

The signal--to--noise ratio in the first cell is therefore

$$\frac{\text{signal}}{\text{noise}} = \frac{\text{(observed value - hypothesised value)} }{\text{scale factor}} = \frac{(n_{11} - e_{11})^2}{e_{11}} = \frac{(45-39)^2}{39}=0.92.$$

The ratios for the other three cells can be calculated in the same way. The test statistic we are interested in should obviously combine the signal--to--noise ratios of all four cells, and this is done by just summing them up, i.e. This is a scary looking formula, but it only means that you need to calculate the sum of all four cell counts, i.e. for each of the two rows ($i$=1 and $i$=2) and each of the two columns ($j$=1 and $j$=2).

$$\chi^2 = \sum_{i = 1}^2 \sum_{j = 1}^2 \frac{(n_{ij} - e_{ij})^2}{e_{ij}}.$$

This is our test statistic, called the $\chi^2$--statistic. For our example we get:

$$\chi^2 = \frac{(45-39)^2}{39} + \frac{(25-31)^2}{31} + \frac{(22-28)^2}{28} + \frac{(28-22)^2}{22}$$

$$ = \frac{36}{39} + \frac{36}{31} + \frac{36}{28} + \frac{36}{22}$$
 
$$ =  `r ((45-39)^2)/39 +
       ((25-31)^2)/31 +
       ((22-28)^2)/28 +
       ((28-22)^2)/22` $$
       
Statistics software, such as R, applies a _Yate's correction_. Using this correction, our test statistic is calculated as 

$$\frac{(|45-39|-0.5)^2}{39} + \frac{(|25-31|-0.5)^2}{31} + \frac{(|22-28|-0.5)^2}{28} + \frac{(|28-22|-0.5)^2}{22} 
= `r round((((abs(45-39)-0.5)^2)/39 +
         ((abs(25-31)-0.5)^2)/31 +
         ((abs(22-28)-0.5)^2)/28 +
         ((abs(28-22)-0.5)^2)/22), 2)`.$$

The calculations in our text didn't take this into account, and also we rounded to the nearest whole number. that was just to make things easier to follow. Using more precise values and Yate's correction yields `r round((((abs(45-39.08333)-0.5)^2)/39.08333 + ((abs(25-30.91667)-0.5)^2)/30.91667 + ((abs(22-27.91667)-0.5)^2)/27.91667 + ((abs(28-22.08333)-0.5)^2)/22.08333), 2)`.

```{r echo = FALSE}

# Original observations
M <- as.table(rbind(c(10, 24), c(15, 1)))
dimnames(M) <- list(Nose = c("Blue","Red"),
                    Antennae = c("Yes","No"))

M2 <- as.table(rbind(c(50, 90), c(50, 10)))
dimnames(M2) <- list(Nose = c("Blue","Red"),
                     Antennae = c("Yes","No"))

Xsq1 <- chisq.test(M2, rescale.p = TRUE)

# Xsq1 <- chisq.test(table(martian$Nose, martian.df$Antennae))

myChi <- chisq.test(martian_chi$Nose.Colour,martian_chi$Antennae)

```

Recall that the test statistic tells us how far away our observed results are from the hypothetical results that we'd get if $H_0$ were true. Under $H_0$ it follows a distribution, and in this case it's the $\chi^2$--distribution. Like the $t$--distribution, it depends on some degrees of freedom, and for $2 \times 2$--contingency tables, this number is always equal to 1. We can therefore write

$$ \chi^2 \stackrel{H_0}{\sim} \chi^2_1.$$

Figure \@ref(fig:chiDist) shows this distribution. The red dashed line shows the critical value that your observed statistic needs to pass when you would like to reject $H_0$ with a 5% significance level, and the (Yates corrected) observed value of the test statistic is shown in black. As we have seen before, the $p$--value corresponds to the area under the curve to the right of the observed value of the test statistic. It is shown in red. Our observed value of the test statistic is bigger than the critical value, we therefore know the $p$--value will be smaller than 5%. The exact value is $p = 0.04$.

<!-- $\Sexpr{1 - pchisq(q = myChi$statistic,df = 1)}$. -->

ONLY ONE GIRE WITH BOTH CRITICAL VALUE AND STATIATIC!!!


The $\chi^2_1$ distribution. The critical value that corresponds to a 5\%  significance level is shown in red. Our observed value of the test statistic is shown in black. Note that we show the test statistic that has been corrected according to Yate (which is explained in the marginnote). The $p$--value corresponds to the area under the curve shown in red.
```{r chi-sq-distribution, echo = FALSE}
x <- seq(0, 8, length = 1000)
dist1 <- dchisq(x, df = 1)
dist2 <- dchisq(x, df = 2)
dist3 <- dchisq(x, df = 3)
dist4 <- dchisq(x, df = 4)
dist6 <- dchisq(x, df = 6)
dist9 <- dchisq(x, df = 9)

df.chix <- data.frame(x, get(ls(pattern = "^dist")))
df.chix <- data.frame(x, dist1, dist2, dist3, dist4, dist6, dist9)
df.chix <- reshape2::melt(df.chix, id.var = "x")
levels(df.chix$variable) <- c("1", "2", "3", "4", "6", "9")

df.chix <- df.chix[df.chix$variable == "1",]

# ggplot(df.chix, aes( x = x, y = value, col = variable)) +
 # geom_line() +

# chi.base <- 

ggplot(df.chix, aes( x = x, y = value)) +
  geom_line(col = DarkYellow) +

  geom_vline(xintercept = qchisq(p = .95, df = 1), col = DarkRed, linetype = 2) + # 3.841459
  # annotate("text", x = qchisq(p = .95, df = 1) - 0.1, y = 0.4, label = "Critical\nvalue", hjust=1) +


  scale_x_continuous("$X^2$", expand = c(0,0), breaks = seq(0,6,2)) +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  scale_colour_manual("Degrees of Freedom", values = rev(brewer.pal(9, "Reds")[4:9])) +

  coord_cartesian(xlim = c(0,6), ylim = c(0,0.5)) +

  # ggtitle("The $\\chi^2_1$ Distribution") +
  theme(axis.line = element_blank(),
        legend.position = c(0.8, 0.8),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black")
  )

# chi.base

myChi <- chisq.test(martian_chi$Nose.Colour,martian_chi$Antennae)

```

```{exercise, name="Defining statistical hypotheses"}
How would you report this result in a paper?\\
```

```{exercise, name="Defining statistical hypotheses"}
Why does the $\chi^2$--distribution have low values around zero and falls off quite quickly for larger values?\\
```

What's next? We have seen this scenario before. When we position our test statistic as the critical value on the $\chi^2$ distribution, the area under the curve beyond our test statisic is the probability of observing our results or something more extreme. As you can see from Figure \@ref(fig:chiDistwithP), our $\chi^2$ is to the right of the critical value. There is a probability of only $`r 1 - pchisq(q = myChi$statistic,df = 1)`$ observations like our or more extreme being observed. Which is exactly the $p$--value we get from a $\chi^2$ test! Note that the precise value, as calculated as per Yate's version, above.

<!-- %  -->
<!-- % \begin{figure} -->
<!-- % <<chi-sq-distribution-withP, echo = FALSE, eval = TRUE, fig.width=6, fig.height=3>>= -->
<!-- %  -->
<!-- % chi.base + -->
<!-- %  -->
<!-- %   geom_ribbon(data=subset(df.chix, df.chix$x > myChi$statistic), aes(x=x,ymax=value),ymin=0,fill=MidRed) + -->
<!-- %   geom_vline(xintercept = myChi$statistic, col = "black", linetype = 1) -->
<!-- %   # annotate("text", x = myChi$statistic + 0.1, y = 0.4, label = "test\nstatistic", hjust=0, vjust = 1) + -->
<!-- %  -->
<!-- %  -->
<!-- % @ -->
<!-- % \label{chiDist=withP} -->
<!-- % \end{figure} -->

One more test under our belt! But before we move on, two final notes on the $\chi^2$--test:

1. The $\chi^2$--test is a large--sample test that requires $n$, i.e. the total counts, to be large. Once again, we encounter the question _how large is large enough?_. For the  $\chi^2$--test, a good rule-of-thumb is that each cell in the contingency table of the expected counts should be at least 5. If your sample size is not large enough, there is a solution that we'll cover later (see \@ref(sec:fishers)).
2. Like the $t$--test, the $\chi^2$--test also gets more precise with larger sample size! 

## Large sample sizes {largechi}

For example, consider the table below as a starting point for the analysis which has the same expected proportions as our example but the observed values multiplied by 10:

_Observed table of counts_

| Nose Colour | No    | Yes   | Total  |
|-------------|-------|-------|--------|
| Red         | $450$ | $250$ | $700$  |
| Blue        | $220$ | $280$ | $500$  |
| Total       | $550$ | $650$ | $1200$ |

The resulting $\chi^2$--statistic would be 47.85 (if calcuated with Yate's correction) which would give a $p$--value of $4.59 \times 10^{-12}$. Like with all the other tests, this means that with larger sample size, you have a bigger chance to detect smaller deviations from $H_0$, i.e. you have a bigger power!