Classical tests, linear models and their extensions for the analysis of 2 × 2 contingency tables

Ecologists and evolutionary biologists are regularly tasked with the comparison of binary data across groups. There is, however, some discussion in the biostatistics literature about the best methodology for the analysis of data comprising binary explanatory and response variables forming a 2 × 2 contingency table.We assess several methodologies for the analysis of 2 × 2 contingency tables using a simulation scheme of different sample sizes with outcomes evenly or unevenly distributed between groups. Specifically, we assess the commonly recommended logistic (generalised linear model [GLM]) regression analysis, the classical Pearson chi‐squared test and four conventional alternatives (Yates' correction, Fisher's exact, exact unconditional and mid‐p), as well as the widely discouraged linear model (LM) regression.We found that both LM and GLM analyses provided unbiased estimates of the difference in proportions between groups. LM and GLM analyses also provided accurate standard errors and confidence intervals when the experimental design was balanced. When the experimental design was unbalanced, sample size was small, and one of the two groups had a probability close to 1 or 0, LM analysis could substantially over‐ or under‐represent statistical uncertainty. For null hypothesis significance testing, the performance of the chi‐squared test and LM analysis were almost identical. Across all scenarios, both had high power to detect non‐null effects and reject false positives. By contrast, the GLM analysis was underpowered when using z‐based p‐values, in particular when one of the two groups had a probability near 1 or 0. The GLM using the LRT had better power to detect non‐null results.Our simulation results suggest that, wherever a chi‐squared test would be recommended, a linear regression is a suitable alternative for the analysis of 2 × 2 contingency table data. When researchers opt for more sophisticated procedures, we provide R functions to calculate the standard error of a difference between two probabilities from a Bernoulli GLM output using the delta method. We also explore approaches to compliment GLM analysis of 2 × 2 contingency tables with credible intervals on the probability scale. These additional operations should support researchers to make valid assessments of both statistical and practical significances.

Keywords: 2 × 2 contingency table; chi‐squared test; linear models; logistic GLMs; uncertainty estimates

From sex differences in survival rates (e.g. Teder & Kaasik, [44]) to taxonomic bias in publications (e.g. Rosenthal et al., [37]), assessing categorical differences in discrete variables is central to many questions in ecology and evolution. A legitimate concern among biologists is how to correctly analyse this type of data. In particular, what statistical tests produce an unbiased estimate of group differences, what assumptions need to be satisfied, and how should one interpret and present the results?

Here, we consider a simple but commonly encountered situation where two samples of binary data are compared. Such data are often presented in a 2 × 2 table, also called a contingency table since the key question is whether there is a contingency between the row and column variables. For example, given a set of sampled rodents, an individual's sex and whether they showed symptoms of some specific disease may be recorded. The natural set of questions to then ask are whether sex and disease status are linked and, if so, what is the strength and direction of the effect. We would also generally like to quantify how confident we can be in our answers to these questions.

Given such a 2 × 2 contingency table of count data, the Pearson chi‐squared test (Pearson, [34]) is widely viewed as a reasonable way to test for differences between groups (Albert, [4]; Altman & Krzywinski, [5]; Crawley, [12]; Dytham, [14]; Fagerland et al., [15]; Seltman, [42]; Whitlock & Schluter, [48]). Variations on the chi‐squared test depending on data structure are also commonly discussed. For example, to account for small sample sizes, Yates' correction for continuity (Yates, [50]) or Fisher's exact test (Fisher, [17]) are sometimes recommended (Altman & Krzywinski, [5]; Crawley, [12]) although not without some dissention (Fagerland et al., [15]; Ruxton & Neuhäuser, [39]). The use of Fisher's exact test, in particular, is controversial given that both margins of the contingency table are rarely naturally fixed (see e.g. Agresti, [1]; Berkson, [9]; Kempthorne, [25]), leading to unnecessarily conservative estimates. The Fisher mid‐p (Lancaster, [26]) test is therefore often suggested to better represent evidence against the null hypothesis (Fagerland et al., [15]; Hwang & Yang, [22]; Routledge, [38]). Alternatively, Lydersen et al. ([28]) recommend the exact unconditional test (Barnard, [7]; Boschloo, [11]) with a Berger and Boos correction (Berger & Boos, [8]) as the 'gold standard' for testing association in 2 × 2 tables.

There are, however, two major disadvantages of using the chi‐squared test (or its alternatives and refinements discussed above) to analyse contingency tables. First, these test statistics only assess the distribution of the response variable under the null hypothesis that there is no difference between groups. They do not provide any information on the strength or direction of association among variables. This is a particularly relevant criticism as ecologists and evolutionary biologists shift away from principally relying on null hypothesis significance testing (NHST) and p‐values (Stephens et al., [43]). Second, these test statistics are somewhat limited in application. When testing for an association between two categorical variables, both the explanatory and response variables must be dichotomous. While this type of data is very common in ecology and evolution, it limits the expansion of the analysis to more complex datasets. For example, if the hypothetical rodent study above also recorded the mass of each individual, this continuously distributed factor cannot be easily incorporated into any of the methods discussed above.

Both these limitations of the chi‐squared test can be overcome by instead adopting regression analyses. Not only does a regression analysis provide more informative insights into the relationship between explanatory and response variables (describing strength and direction of effects), the precision of an estimate is also easily assessed. Results can thus be more intuitive, that is, offering an effect size and uncertainty (mean ± SE). A regression analysis is also more versatile and can be expanded beyond the 2 × 2 contingency table. Using this single methodology, a biologist can investigate combinations of quantitative and categorical explanatory variables hypothesised to have an important influence on the response variable. Despite these obvious advantages, there seems to be some confusion and discussion in the literature about the appropriateness of regression analyses to assess the association between binary variables.

Foremost, biologists are widely discouraged from using linear regressions to analyse contingency tables. The most commonly cited cause of concern is that, for the analysis of data comprising binary response and categorical explanatory variables, the assumptions of normality and homogeneity of variance are violated (Kaplan, [24]; Seltman, [42]; Tutz, [45]). Using a linear regression is thus regarded by some as 'unsatisfactory' (Tutz, [45]) at best and 'completely unreliable' (Seltman, [42]) at worst. By contrast, the logistic regression is largely accepted as an appropriate method to estimate probabilities of categorical explanatory variables, compensating for the aforementioned problems with linear regressions (Dunn & Smyth, [13]; Fagerland et al., [15]; Lever et al., [27]; Orme & Combs‐Orme, [33]; Ramos et al., [36]; Tutz, [45]). Nonetheless, a clear advantage of linear over logistic regression is the often more straightforward interpretation of model coefficients (without requiring conversion).

When assessing the possible differences between groups in dichotomous variables, biologists tend to ask about differences in probabilities between groups. In the linear regression of a binomial outcome on a group variable, a one‐unit increase in the explanatory x increases the conditional expectation of the response y by β units. By contrast, the coefficients of logistic regressions are given in log‐odds, which will often be a less intuitive unit for researchers asking biological questions. Even after converting the coefficient for the explanatory variable back to an ordinary scale, this odds ratio does not represent a constant increase or decrease in the response y given the explanatory variable x. This can make results difficult to contextualise (Gallis & Turner, [18]; Halvorson et al., [20]), even more so by non‐experts who might be unfamiliar with odds as a statistical measure of the probability of one outcome versus the other (Grant, [19]; Schwartz et al., [41]).

Here, we assess the appropriateness of different statistical tests for the analysis of data comprising binary explanatory and response variables (both coded as 0,1). In contrast to many previous works on this subject, we focus on different goals of statistical analysis beyond p‐values, namely estimation and statements about uncertainty in differences in probabilities. Odds, odds ratios and log‐odds ratios are all valid ways of presenting results and are well treated in many works. We focus specifically on differences in probabilities, in extension of typical considerations, because this is the scale most relevant and intuitive to biologists and the scale that will frequently map onto how biologists formulate research questions. Using a simulation scheme of different sample sizes with n outcomes evenly or unevenly distributed between groups, we assess the biasedness of inferences of (a) the difference in proportions between groups and (b) estimates of statistical uncertainty (standard errors and confidence intervals) when using linear and logistic regressions. We also assessed and compared (c) the false‐positive rate and (d) the consistency of p‐values for the Pearson chi‐squared test and its alternatives, linear regression and logistic regression.

Each simulation generated data for two groups. These groups could represent control vs. treatment, male vs. female, etc. Sample size in each group was determined according to four schemes. In the first, the total sample size of $n_t$ was divided evenly between the two groups ( $n_t$ /2). In the second, $n_t$ was apportioned between groups following a binomial distribution with probability 0.5. In the third scenario, the groups had sample sizes 0.8 $n_t$ and 0.2 $n_t$ . Finally, the total of $n_t$ samples was partitioned between groups following a binomial distribution with probability 0.2 for group x = 0. When data were distributed randomly between groups, if a group had fewer than two observations, we modified the values of x so that there were at least two cases in each group. This was done on the grounds that most researchers would not look to explore contingency between groups when one of the groups had such a small sample size.

We then simulated data to represent a Bernoulli response variable for all the observations in both groups. This variable $y$ took values of 1 and 0, which could represent outcomes such as survived vs. died, mated vs. unmated, etc. Observations from group x = 0 were assigned a probability $p_0$ of success (i.e. of y = 1, rather than y = 0), and the n values for y were drawn as independent samples from a Bernoulli distribution with probability $p_0$ . The probability of success for observations in group x = 1 was defined by  $p_1=p_0+\delta$ , and values of y for individuals in group x = 1 were drawn from a Bernoulli distribution with probability of success of $p_1$ . Formally, the simulation of the y data can be described according to1 $$y_i~\text{Bernoulli}\left(p_0+\delta x_i\right),$$ where $y_i$ represents success or failure (e.g. survived vs. died) coded as 1 and 0, respectively, with $x_i$ coding the group membership (e.g. male vs. female; 0, 1) of individual i. $p_0$ is the probability of success in the first group (x = 0). $\delta$ is the difference in success probability between the two groups.

For each simulation scheme, we simulated a range of sample sizes of $n_t$ (the total number of samples across both groups) from [10, 20, 30, 50, 70, 100]. We considered two sets of values of $p_0$ in combination with $\delta$ . First, we considered a value of $p_0$ of 0.5 and values of $\delta$ between −0.5 and +0.5 in increments of 0.1. Second, we considered a value of $p_0$ of 0.1 and values of $\delta$ from −0.1 to +0.9, also in increments of 0.1.

For every simulated dataset of values of x and y, we conducted a range of analyses focused on estimation of the parameter $\delta$ , that is, of the difference in the underlying probability of success (of y = 1) in the two groups coded x = 0 and x = 1. First, we conducted linear model (LM) analyses of the dependence of y on x. We fitted the model2 $$y_i={\alpha }_{\mathrm{LM}}+{\beta }_{\mathrm{LM}}{x}_i+e_i,$$ where $y_i$ is the outcome (0, 1) for individual i, from group $x_i$ (0, 1). ${\alpha }_{\mathrm{LM}}$ is the intercept and ${\beta }_{\mathrm{LM}}$ is the slope of the linear regression of y on x, which estimate $p_0$ and $\delta$ , respectively. We use the standard ordinary least squares (OLS) standard error (SE) of ${\beta }_{\mathrm{LM}}$ as the standard error of $\delta$ , and use a Wald‐type confidence interval (estimate ± 1.96 SEs) as our LM‐based confidence interval (CI) for $\delta$ . We use the t‐test‐based p‐value for ${\beta }_{\mathrm{LM}}$ from the summary.lm() function in base R as our LM‐based p‐value for $\delta$ .

Secondly, we conducted a generalised linear model (GLM) analysis with a logit link function and a Bernoulli response, also known as a logistic regression. This model was constructed as3a $${\eta }_i={\alpha }_{\mathrm{GLM}}+{\beta }_{\mathrm{GLM}}{x}_i$$ 3b $$y_i~\text{Bernoulli}\left(g^{-1}\left({\eta }_i\right)\right),$$ where the data are modelled as a linear function on the scale of a latent variable $\eta$ , related to expected probabilities via the inverse link function $g^{-1}\left(\right)$ (i.e. the inverse of the logit link function), defined $g^{-1}\left(a\right)=\frac{e^a}{1+e^a}$ . Bernoulli sampling of the data $y_i$ was performed given the probabilities $g^{-1}\left({\eta }_i\right)$ . As such a GLM‐based estimator of $\delta$ is4 $${\widehat{\delta }}_{\mathrm{GLM}}=g^{-1}\left({\alpha }_{\mathrm{GLM}}+{\beta }_{\mathrm{GLM}}\right)-g^{-1}\left({\alpha }_{\mathrm{GLM}}\right).$$

${\widehat{\delta }}_{\mathrm{GLM}}$ is thus a somewhat complex quantity, depending not only on the GLM parameter that describes the difference between the groups ( ${\beta }_{\mathrm{GLM}}$ ) but also on the model intercept ( ${\alpha }_{\mathrm{GLM}}$ ). This renders its interpretation more difficult than the estimate from the LM analysis. While biological inferences on the logit data scale are possible, it would be useful if estimates of differences in probabilities (as in Equation 4) were used more widely. Being able to put uncertainty on the probability scale would be useful as well. A standard error for the ${\widehat{\delta }}_{\mathrm{GLM}}$ estimator, constructed by the delta method (Lynch & Walsh, [29]; Ver Hoef, [47]), is5 $$\mathrm{SE}\left[{\widehat{\delta }}_{\mathrm{GLM}}\right]\approx \sqrt{{\left[\right]}^T{\mathrm{\Sigma }}_{\mathit{\alpha \beta },\mathrm{GLM}}\left[\right]},$$ where $\raisebox{1ex}{$\partial {\widehat{\delta }}_{\mathrm{GLM}}$}\!\left/ \!\raisebox{-1ex}{$\partial {\widehat{\alpha }}_{\mathrm{GLM}}$}\right.$ and $\raisebox{1ex}{$\partial {\widehat{\delta }}_{\mathrm{GLM}}$}\!\left/ \!\raisebox{-1ex}{$\partial {\widehat{\beta }}_{\mathrm{GLM}}$}\right.$ are the derivatives of ${\widehat{\delta }}_{\mathrm{GLM}}$ in Equation 4 with respect to the logistic intercept and contrast estimated by the GLM (as specified in Equation 3a,b), evaluated at their estimated values. T is the transpose operator.  ${\mathrm{\Sigma }}_{\mathit{\alpha \beta },\mathrm{GLM}}$ is the covariance matrix of the GLM's parameter estimates. The square roots of the diagonal elements are the standard errors of the logistic intercept and contrast terms.

It is relatively straightforward to generate a CI from a binomial GLM on the probability scale for a conditional probability (e.g. for one group or the other, in the present context). One would typically generate a Wald‐type CI on the linear predictor scale and transform its upper and lower limits to the probability scale with the inverse logit function. However, this strategy is not generally possible for differences in conditional probabilities. Consequently, similarly to how SEs for differences on the probability scale are not generally considered—useful as they would be—CIs for differences or data‐scale effects in GLMs are not routinely generated by biologists.

There are many methods for setting CIs, and our aim was not to provide a comprehensive overview. We did, however, investigate the utility of the three most commonly accepted ways of constructing CIs for statistical models (Alan, [3]; Fagerland et al., [15]; Newcombe, [32]). First, we tested the simple probability scale Wald‐type CIs constructed as 1.96 times the delta‐SE (Equation 5) above and below the estimated difference (Equation 4). Second, we assessed the Wilson or score CIs (Wilson, [49]), making use of an adapted z2stat function from the R package PropCIs for our calculations (Scherer, [40]). Finally, we used profile likelihood‐based CIs (Venzon & Moolgavkar, [46]) of $\delta$ . To generate the profile likelihood CI, we fixed $\delta$ to a range of values between [−0.99 and 0.99] and found the value of $p_0$ that maximised the likelihood of observing the data, given the fixed value of $\delta$ . We then recorded the upper and lower values of $\delta$ for which twice the log‐likelihood difference from the unconstrained model was less than the 95% quantile of a chi‐squared distribution with one degree of freedom. We assessed significance of the ${\beta }_{\mathrm{GLM}}$ using both the standard z‐test output from the summary.glm() function in base R and the likelihood ratio test (LRT). To implement the LRT, we generated a p‐value under the assumption that twice the difference in the log‐likelihood of the GLM in Equation 3a,b with the log‐likelihood of an intercept‐only GLM is chi‐squared distributed with one degree of freedom.

Finally, we conducted four classical contingency table analyses on the 2 × 2 contingency tables of x and y values. These analyses do not estimate the value of $\delta$ , but can be interpreted as tests of statistical significance of $\delta$ , that is, of whether its value differs from a (null) hypothetical value of zero (no difference between $p_0$ and $p_1$ ). The first two such tests were chi‐squared tests, without (Pearson, [34]) and with Yates' continuity correction (Yates, [50]). The third test was a Fisher's exact test (Fisher, [17]). In addition to the exact test's ordinary p‐value, we calculated the mid‐p (Lancaster, [26]) using the R package epitools (Aragon, [6]). The final test was an exact unconditional test (Barnard, [7]; Boschloo, [11]) with the Berger and Boos correction (Berger & Boos, [8]), implemented using the R package exact2x2 (Fay & Hunsberger, [16]).

We first assessed results of linear and logistic regression analyses in terms of biasedness. Bias is the difference between the average value of an estimate and the true value of the quantity it is estimating. Formally,6 $$\text{bias}\left[\widehat{\delta }\right]=E\left[\widehat{\delta }\right]-\delta ,$$ where $\delta$ is the true value of the difference in probabilities and $E\left[\widehat{\delta }\right]$ is the expected value of the estimator of the difference. As such, in unbiased analyses, the average value of the estimator $\delta$ applied to replicate simulated datasets would not differ from the true value in those simulations.

While bias is typically defined as a difference between the expected value of an estimator and the true value of the estimand (as in Equation 6), we expressed (un)biasedness in a proportional sense for our assessment of different methods for generating standard errors of differences in probability. Specifically, we divided the average SE for any given method by the empirical SD of the estimator,7 $${\text{bias}}_{\text{proportional}}\left[\widehat{\delta }\right]=\frac{E\left[\mathrm{SE}\left[\widehat{\delta }\right]\right]}{\mathrm{SD}\left[\widehat{\delta }\right]},$$ where $E\left[\mathrm{SE}\left[\widehat{\delta }\right]\right]$ is the average standard error across replicate simulations, and $\mathrm{SD}\left[\widehat{\delta }\right]$ is the standard deviation of the estimated differences in probabilities across replicate simulations in any given scenario. As such, a value of one for the proportional bias measure for the standard error indicates ideal behaviour.

Similarly, we assessed confidence intervals by considering their coverage properties. For a given scenario, a 95% CI has the correct or unbiased (nominal) coverage when estimates fall within its bounds in 95% of replicate simulations.

All analyses were done in R (R Core Team, [35]). Code for the simulations, analyses, and figures presented in this manuscript can be found on GitHub and Zenodo (Nagel et al., [30]; https://github.com/rebebba/ProbUncertainty_MSCode).

All results were very similar among the four data‐generating schemes. In other words, whether the distribution of the predictor variable was fixed or random, or whether the data were balanced or not between the two groups, most results were very similar. Therefore, except when notable differences occur, we report on results for the simulations with fixed and equal proportions of the data between the two groups. Results for the other three scenarios are presented in the Supplementary Materials.

The two main analyses that generate estimates of $\delta$ are the LM analysis, which estimates $\delta$ directly, and the GLM analysis, from which an estimator of $\delta$ can be recovered using Equation 4. Both the LM and the GLM returned unbiased estimates of $\delta$ in the sense that the average value of the estimator was equal to the true value (Figure 1; also see Figures A1–A3) across all parameter value combinations that we considered.

Providing that observations were reasonably evenly divided among groups and that sample sizes were not extremely small (i.e. >5 observations per group and $n_t$  > 10) standard errors of $\delta$ from LM and GLM analyses were good reflections of the uncertainty in the estimation process (Figure 2; also see Figure A4). However, when observations were extremely unevenly distributed between groups, standard errors from the LM analysis could substantially over‐ or under‐represent statistical uncertainty (Figure 3; also see Figure A5).

For the LM analysis, coverage properties largely reflected the behaviour of standard errors. Performance was good providing a reasonable sample size ( $n_t$  ≥ 20), and only deteriorated appreciably when data were highly unequal across groups and there were large differences in the true probabilities between groups (Figure 4a,e; also see corresponding plots in Figures A6–A8). The Wald‐type CIs we generated based on the delta‐SEs from the GLM analysis were generally conservative (Figure 4b,f), indicating greater uncertainty in the GLM estimates of differences between groups than was actually achieved. The score (Figure 4c,g) and profile likelihood (Figure 4d,h) methods for generating CIs for $\delta$ performed reasonably well across all scenarios that we simulated (also see corresponding plots in Figures A6–A8).

Under a balanced experimental design, the chi‐squared test, mid‐p, exact unconditional and LM analysis had a type 1 error rate very close to the conventionally accepted 5% threshold (Figure 5 at $\delta$  = 0; also see corresponding plots in A10). Yates' correction, Fisher's exact and the GLM analysis with z‐test p‐values were generally most conservative. When $p_0$  = 0.1, observations were very unevenly distributed between groups, and sample sizes were small ( $n_t$  ≤ 20) the type 1 error rate for both the chi‐squared test and LM analysis was closer to 10% (Figure A9). Across most scenarios, the highest type 1 error rate was realised using a GLM analysis with p‐values generated using the LRT (Figure 5 at $\delta$  = 0; also see corresponding plots in Figures A9–A11).

Across the range of scenarios considered, the chi‐squared test, LM analysis and GLM analysis with p‐values generated using the LRT generally gave the best performance, with the highest power to detect non‐null effects (Figure 5; also see corresponding plots in Figures A9–A11). The poorest performance was realised in the GLM analysis with default p‐values (z‐test based) from the summary.glm() function in base R.

The underlying p‐values generated from the linear models and chi‐squared tests were nearly equivalent, both across the full dataset and at the α = 0.05 level (Figure 6; also see Figures A12–A14).

Many biological phenomena have binomial outcomes: diagnostic tests may be passed or failed, individuals may survive or die, sites may or may not be occupied. Consequently, many biological questions boil down to investigating differences in probabilities. Here, we have investigated a range of approaches to characterise differences between groups. Our primary focus was on the estimation of the difference in probability between two groups and the ability to make reasonable statements about uncertainty in those differences. To this end, we first considered the performance of LMs and Bernoulli GLMs to produce unbiased estimates and meaningful standard errors and confidence intervals. Secondly, given that widely recommended approaches for the analysis of binary data (e.g. GLMs) do not directly generate standard errors or confidence intervals for differences in probabilities, we developed and tested methods to do so. Finally, we compared model‐based approaches and classical tests for 2 × 2 contingency table analysis in terms of the performance of p‐values associated with the null hypothesis of no difference between two groups. We discuss the main observations from these exercises in terms of what they mean for practicing scientists, and how analysing differences in probabilities might fit into the process of learning practical statistics.

A post‐doc is looking over a study area with the long‐term PI on their new project:

• New post‐doc Do the unicorns in the valley have higher survival than those on the hills?
• Experienced PI Hmm, yes. Their log‐odds of survival probably differs by 1.5.
• New post‐doc As much as 1.5 on the log‐odds scale?
• Experienced PI Indeed; maybe as much as 2!

A conversation that definitely never happened.

Almost certainly, the underlying interests of both participants in this fanciful conversation would have been about survival probability and the difference in those probabilities between habitats—not about differences in log‐odds ratios. This is not to discount the fact that the log‐odds are a statistically convenient scale, and not without biological importance. Rather, we think it likely that most biologists prefer to consider their data and results in terms of probabilities, finding them a more intuitive and illustratable statistic (Gallis & Turner, [18]; Halvorson et al., [20]). That, in the above example, the post‐doc made a judgement about biological effect size and the PI made a statement about uncertainty, both on the log‐odds ratio scale, is pretty improbable. We therefore consider how best to support this inclination and focus on the estimation of differences in probabilities between groups and what we can say about uncertainty in those estimates.

First, we focus on the linear regression, which gives estimated differences in probabilities and its standard errors directly as a model output. While many biostatistics texts give the impression that using a linear model to analyse binomial data would be dangerously naïve (Kaplan, [24]; Lever et al., [27]; Seltman, [42]; Tutz, [45]), we found that a linear model with a binary outcome (0, 1) regressed on a binary explanatory variable (0, 1) gave unbiased estimates of differences in probabilities. More specifically, the LM accurately estimated the true difference between the means of two groups (Figure 1). Standard errors from the LM did not, however, perfectly reflect the true uncertainty in estimated differences (Figures 2 and 3). Basic Wald‐type confidence intervals (estimate ± 1.96 SEs) performed well under a broad range of circumstances but were compromised when sample size was small and data distributions were very skewed (Figure 4). Nonetheless, SEs and CIs from the LM analysis were generally of the correct order of magnitude. Furthermore, the conditions under which the LM‐based SEs and CIs performed poorly are the same conditions under which a chi‐squared test is generally discouraged (Crawley, [12]; Fagerland et al., [15]; Kang et al., [23]; Ruxton & Neuhäuser, [39]).

Where linear regression is discouraged for the analysis of contingency tables, logistic regression is widely recommended as an ideal means of modelling binomial data (Dunn & Smyth, [13]; Lever et al., [27]; Orme & Combs‐Orme, [33]; Ramos et al., [36]; Tutz, [45]). While a binomial GLM analysis does not directly return inference of differences in probabilities, the conversion of predictions to the probability scale is reasonably commonplace. Calculating standard measures of uncertainty (SEs and CIs) on the probability scale is, however, less common. For this reason, we investigated a method that is standard in some areas of statistics, but not routinely utilised by biologists, to derive SEs of differences in probabilities: the delta method. We find that this linear approximation of the transformation from the log‐odds scale to the difference in probabilities (Equation 5) performs reasonably well. It gave accurate SEs of differences in probabilities and was outperformed by the OLS‐SEs only in a few circumstances involving the smallest sample sizes, and then only modestly so (Figure 2). The GLM delta‐SEs also performed well in some scenarios where the LM analysis performed poorly, namely when differences between groups were large and one group had a true value very near 0 or 1 (Figure 3).

Similar to SEs of differences in probabilities, generating CIs of differences in probabilities from binomial GLMs is not standard practice in the fields of ecology and evolution. We tested three methods to do so. First, we applied the standard Wald‐type method (estimate ± 1.96 SEs) using the delta‐SEs. This generated CIs that were typically overly conservative, that is, that depicted greater ranges of uncertainty in the estimate of differences than actually occurred (Figure 4). This is undesirable since excessive conservatism can equate to wasting sample size, money and effort. We then tested the score and profile likelihood‐based methods for generating CIs. Neither are directly linked to LM‐ or GLM‐based inferences, but both focus on estimated differences in probabilities and make statements about uncertainty of those differences. Both methods generated CIs with approximately correct coverage properties across the range of scenarios investigated (Figure 4).

While our primary focus was on the estimation of differences in probabilities and the generation of statements about uncertainty to accompany these estimates, we also considered two p‐value‐based assessments of model performance: power to reject false null hypotheses and false‐positive rates in the case of no effect. Given that most classical approaches for analysing contingency tables generate p‐values but not estimates or uncertainty statements about differences in probabilities, this also allowed us to compare regression model outputs with four commonly recommended alternatives and refinements to the chi‐squared test. We offer some general conclusions drawn from this comparison.

The power and false‐positive rates from the LM analysis and the Pearson chi‐squared test were nearly identical (Figures 5 and 6). This is inevitable, given that the underlying p‐values generated by the two approaches are also nearly identical. The p‐values from the LM analysis assume that estimation errors follow a t‐distribution, which rapidly converges on a normal distribution when sample size ≥10 (Seltman, [42]). The p‐values from the chi‐squared test arise from a normal approximation to deviations from the null model of equal probabilities (Lydersen et al., [28]).

The GLM‐based approach was the most powerful NHST when H0 was false; it was not outperformed under any combination of simulation parameters assessed (Figure 5). Also, the GLM false‐positive rate most closely reflected the nominal rate at α = 0.05, except in some cases with the smallest sample size. However, these properties are only true for GLM p‐values generated from the LRT. The default statistical hypothesis test for GLMs, the z‐test (i.e. the standard summary.glm() output in base R), can severely lack power when one or both groups have observed rates near 0 or 1; it can also be excessively conservative under the null hypothesis. This is a known phenomenon (Bolker et al., [10]). It is thus widely recommended to apply the LRT to binomial GLMs, especially when some cells of an experiment have very high or low probabilities (Agresti, [2]; Hauck & Donner, [21]). It is worth noting, however, that the LRT is not perfect. It gives modestly inflated type 1 error rates (i.e. p‐values <0.05 when the true difference between groups is 0) in a range of circumstances and is sometimes excessively conservative.

All chi‐squared test refinements and alternatives (Yates' correction, Fisher's exact test, mid‐P and exact unconditional; see Section 2.3) produced results either more conservative than or comparable to the LM and Pearson chi‐squared analyses (Figure 5). More specifically, they were excessively conservative when the null hypothesis was true and generally less powerful when the null hypothesis was false. Being excessively conservative or aiming to never conduct a test that might have a false positive rate slightly above the stated rate may initially seem laudable. However, conservatism when the null hypothesis is true is inevitably linked to lower power when the null hypothesis is false, and it could be argued that employing tests that underestimate the power and precision of an estimate is an inefficient use of resources.

Broadly speaking, across the range of circumstances that we simulated, p‐values from the LM analysis (or very nearly equivalently, from the Pearson chi‐squared test) were more often near the nominal value of α = 0.05 than any other test, including the often‐recommended GLM‐LRT approach. However, marginal differences between tests under most simulated circumstances suggest that adopting the GLM‐LRT approach on this basis risks adding unnecessary confusion to the analysis of contingency tables. While rigorously avoiding excessive false positives is of course desirable, it can be taken to extremes. These results also highlight the problems with expressing results only in terms of p‐values (Stephens et al., [43]): None of the methodologies we assessed proved perfect for NHSTing under all simulated conditions.

The LM analysis robustly estimated differences of proportions between groups. Measures of statistical uncertainty (SEs and associated CIs) were largely correct representations of the true uncertainty, although somewhat compromised when total sample size was small and strongly unbalanced between the groups. However, these are the same conditions under which a chi‐squared test is generally discouraged (Crawley, [12]; Fagerland et al., [15]; Kang et al., [23]; Ruxton & Neuhäuser, [39]). As such, the LM analysis of 2 × 2 contingency table data is equivalent to (in power and performance under the null) or superior to (easily generated estimates and uncertainty statements) the chi‐squared analysis. This suggests at least three broad benefits to more widespread application of the LM approach to 2 × 2 contingency table analysis.

First, instruction in biological statistics may be well served by shifting to a LM‐based approach. However, given the ubiquity of advice that standard linear regressions are unsuited to data such as arises in a 2 × 2 contingency table analysis, biostatistics instructors face a quandary. One option is to teach LMs for everything except contingency tables and then cover all statistical procedures typically handled in introductory statistics (e.g. the chi‐squared test and its alternatives). A second option has instructors teaching GLM analyses. This would undoubtedly benefit students, but most curricula do not afford enough time to teach GLMs well and this leads to greater problems down the road. By one estimation, 58% of GLMs in ecology and evolution are inappropriate in some way (Bolker et al., [10]). Alternatively, our results show that the typical advice is wrong and LMs are reasonable alternatives to the chi‐squared test. This should allow basic LMs to be used for model‐based teaching covering the full range of classical tests (chi‐squared, t‐test, simple and multiple regression, ANOVA, ANCOVA, etc.).

Second, our results suggest that simple LMs provide a good balance between robustness and appropriateness of the methods, and ease of application and interpretation. LMs produce direct parameter estimates and unbiased uncertainty statements of differences in probabilities without common errors in application or interpretation. While the modern biologist may opt for more sophisticated procedures (e.g. GLM), confidence that more straightforward methods are indeed quite robust is useful. Furthermore, researchers can take comfort in the fact that results from those who prefer simpler methods (e.g. LM) are likely reliable. This view may be particularly useful to meta‐analysts, as inclusion of results from relatively simple analyses in synthetic works may be perfectly justified.

Third, while GLM analysis can generally be used to good effect for basic 2 × 2 contingency table analyses, it is only with substantial additional effort that one can extract the basic biological information that the LM analysis generates directly. To our surprise, no general biostatistics source that we consulted discussed how to generate a standard error for the difference in probability between two groups from a fitted GLM analysis. To this end, for those interested in implementing the linear approximation of the sampling error in the difference in probabilities (delta method SE) in the case of the 2 × 2 contingency table analysis, we provide a function in the form of a GitHub R package named ProbUncertainty (Nagel et al., [31]; https://github.com/rebebba/ProbUncertainty). We also provide functions to calculate the score and profile likelihood CIs.

We have confirmed the consensus in the biostatistics literature that GLMs are well suited to the analysis of 2 × 2 contingency table data. However, our results highlight that substantial care is needed for their application and interpretation. Little consideration has previously been given to the fact that the GLM does not readily output uncertainty in differences between groups in terms of probabilities, despite the fact that such information might correspond best to how biologists formulate research questions. We have illustrated ways to rectify this issue. Perhaps more surprising to many readers, we found that the direct inference of differences in probability and associated uncertainly (SEs and CIs) generated by LM analyses perform well in all circumstances where the basic chi‐squared test can be recommended. This realisation, that the LM analysis is generally reasonable and directly yields the kind of information most people want to know, should be helpful when designing instruction in introductory biostatistics courses and could be broadly useful to researchers. Finally, while statisticians have been industrious in inventing alternatives to the basic Pearson chi‐squared test, we found no circumstances where the commonly recommended alternatives unambiguously outperformed the classical test (or the LM analysis).

Michael B. Morrissey and Graeme D. Ruxton conceived the ideas and designed the methodology. Michael B. Morrissey and Rebecca Nagel drafted the simulation code. Rebecca Nagel ran the final versions of all simulations, made final versions of all graphs and led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Rebecca Nagel was supported by a German Research Foundation grant (515410943). Michael B. Morrissey was funded by a Royal Society University Research Fellowship.

The authors declare no competing interests.

The peer review history for this article is available at https://www.webofscience.com/api/gateway/wos/peer‐review/10.1111/2041‐210X.14318.

All data were simulated. Code for the simulations, analyses and figures presented in this manuscript can be found on GitHub (https://github.com/rebebba/ProbUncertainty%5fMSCode) and Zenodo (https://zenodo.org/doi/10.5281/zenodo.10807611). The code and a worked example for the introduced R package can also be found on GitHub (https://github.com/rebebba/ProbUncertainty) and Zenodo (https://zenodo.org/doi/10.5281/zenodo.10807613).

GRAPH: Figure A1. As for Figure 1, but with an unbalanced design.Figure A2. As for Figure 1, but with a type 3 experimental design.Figure A3. As for Figure 1, but with an unbalanced type 3 experimental design.Figure A4. As for Figure 2, but with a type 3 experimental design.Figure A5. As for Figure 3, but with a type 3 experimental design.Figure A6. As for Figure 4, but with an unbalanced design.Figure A7. As for Figure 4, but with a type 3 experimental design.Figure A8. As for Figure 4, but with an unbalanced type 3 experimental design.Figure A9. As for Figure 5, but with an unbalanced design.Figure A10. As for Figure 5, but with a type 3 experimental design.Figure A11. As for Figure 5, but with an unbalanced type 3 experimental design.Figure A12. As for Figure 6, but with an unbalanced design.Figure A13. As for Figure 6, but with a type 3 experimental design.Figure A14. As for Figure 6, but with an unbalanced type 3 experimental design.