Recall the bell-shaped (standard) normal distribution with mean 0 and variance 1.

Although continuous in nature, it can still be a useful approximation for many discrete random variables formed from sums and means.

Suppose \(X_1,\ldots,X_n\) are a random sample from a Bernoulli distribution with \(Pr(X_i=1)=1-Pr(X_i=0)=\pi\) so that \(E(X_i)=\pi\) and \(V(X_i)=\pi(1-\pi)\). By the CLT,

\(\dfrac{\overline{X}-\pi}{\sqrt{\pi(1-\pi)/n}} \)

has an approximate standard normal distribution if \(n\) is large. "Large" in this context usually means the counts of 1s and 0s (successes and failures) should be at least 5, although some authors suggest at least 10 to be more conservative. Note that since the \(X_i\) take on the values 1 and 0, the sample mean \(\overline{X}\) is just the sample proportion of 1s (successes).

The chi-square distribution with \(\nu\) degrees of freedom can be defined as the sum of the squares of \(\nu\) independent standard normal random variables. In particular, if \(Z\) is standard normal, then \(Z^2\) is chi-square with one degree of freedom. For the approximation above, we have (with \(Y=\sum_i X_i\)) that 

\( \left(\dfrac{\overline{X}-\pi}{\sqrt{\pi(1-\pi)/n}}\right)^2 = \left(\dfrac{Y-n\pi}{\sqrt{n\pi(1-\pi)}}\right)^2 \)

is approximately chi-square with one degree of freedom, provided \(n\) is large. The advantage of working with a chi-square distribution is that it allows us to generalize readily to multinomial data when more than two outcomes are possible.

For the discussion here, we assume that \(X_1,\ldots,X_n\) are a random sample from the Bernoulli (or binomial with \(n=1\)) distribution with success probability \(\pi\) so that, equivalently, \(Y=\sum_i X_i\) is binomial with \(n\) trials and success probability \(\pi\). In either case, the MLE of \(\pi\) is \(\hat{\pi}=\overline{X}=Y/n\), with \(E(\hat{\pi})=\pi\) and \(V(\hat{\pi})=\pi(1-\pi)/n\), and the quantity

\(\dfrac{\overline{X}-\mu}{\sigma/\sqrt{n}} \)

is approximately standard normal for large \(n\). The following approaches make use of this. 

If wish to test the null hypothesis \(H_0\colon \pi=\pi_0\) versus \(H_a\colon \pi\ne\pi_0\) for some specified value \(\pi_0\), the Wald test statistic uses the MLE for \(\pi\) in the variance expression:

For large \(n\), \(Z_w\) is approximately standard normal (the MLE is very close to \(\pi\) when \(n\) is large), and we can use standard normal quantiles for thresholds to determine statistical significance. That is, at significance level \(\alpha\), we would reject \(H_0\) if \(|Z_w|\ge z_{\alpha/2}\), the upper \(\alpha/2\) quantile. The reason for the absolute value is because we want to reject \(H_0\) if the MLE significantly differs from \(\pi_0\) in either direction.

Note that \(\pi_0\) will not be rejected if

\(\displaystyle -z_{\alpha/2}<\dfrac{\hat{\pi}-\pi_0}{\sqrt{\hat{\pi}(1-\hat{\pi})/n}} < z_{\alpha/2} \)

Rearranging slightly, we can write the equivalent result

\(\displaystyle \hat{\pi} -z_{\alpha/2}\sqrt{\dfrac{\hat{\pi}(1-\hat{\pi}}{n}} < \pi_0 < \hat{\pi} +z_{\alpha/2}\sqrt{\dfrac{\hat{\pi}(1-\hat{\pi}}{n}} \)

The limits are the familiar \((1-\alpha)100\%\) confidence interval for \(\pi\), which we refer to as the Wald interval:

\(\displaystyle \hat{\pi} \pm z_{\alpha/2}\sqrt{\dfrac{\hat{\pi}(1-\hat{\pi}}{n}} \  \text{(Wald Interval)}\)

The most common choice for \(\alpha\) is 0.05, which gives 95% confidence and multiplier \(z_{.025}=1.96\). We should also mention here that one-sided confidence intervals can be constructed from one-sided tests, although they're less common.

Another popular way to test \(H_0\colon \pi=\pi_0\) is with the Score test statistic:

Considering that a hypothesis test proceeds by assuming the null hypothesis is true until significant evidence shows otherwise, it makes sense to use \(\pi_0\) in place of \(\pi\) in both the mean and variance of \(\hat{\pi}\). The score test thus rejects \(\pi_0\) when \(|Z_s|\ge z_{\alpha/2}\), or equivalently, will not reject \(\pi_0\) when

\(\displaystyle -z_{\alpha/2}<\dfrac{\hat{\pi}-\pi_0}{\sqrt{\pi_0(1-\pi_0)/n}} < z_{\alpha/2} \)

Carrying out the score test is straightforward enough when a particular value \(\pi_0\) is to be tested. Constructing a confidence interval as the set of all \(\pi_0\) that would not be rejected requires a bit more work. Specifically, the score confidence interval limits are roots of the equation \(|Z_s|- z_{\alpha/2}=0\), which is quadratic with respect to \(\pi_0\) and can be solved with the quadratic formula. The full expressions are in the R file below, but the center of this interval is particularly noteworthy (we'll let \(z\) stand for \(z_{\alpha/2}\) for convenience):

\( \displaystyle \dfrac{\hat{\pi}+z^2/2n}{1+z^2/n} \)

Note that this center is close to the MLE \(\hat{\pi}\), but it is pushed slightly toward \(1/2\), depending on the confidence and sample size. This helps the interval's coverage probability when the sample size is small, particularly when \(\pi\) is close to 0 or 1. Recall that we're using a normal approximation for \(\hat{\pi}\) with mean \(\pi\) and variance \(\pi(1-\pi)/n\). The images below illustrate how this approximation depends on values of \(\pi\) and \(n\). 

Problems arise when too much of the normal curve falls outside the (0,1) boundaries allowed for \(\pi\). In the first case, the sample size is small, but \(\pi=0.3\) is far enough away from the boundary that the normal approximation is still useful, whereas in the second case, the normal approximation is quite poor. The larger sample size in the third case decreases the variance enough to offset the small \(\pi\) value.

In practice, the results of a poor normal approximation tend to be intervals that include values outside the range of (0,1), which we know cannot apply to \(\pi\). The score interval performs better than the Wald in these situations because it shifts the center of the interval closer to 0.5. Compare the score interval limits (below) to those of the Wald when applied to the smartphone data.

\( (0.1455, 0.5190) \)

ci = function(y,n,conf)
{   pi.hat = y/n
    z = qnorm(1-(1-conf)/2)
    wald = pi.hat+c(-1,1)*z*sqrt(pi.hat*(1-pi.hat)/n)
    score = (pi.hat+z^2/2/n+c(-1,1)*z*sqrt(pi.hat*(1-pi.hat)/n+z^2/4/n^2))/(1+z^2/n)
    cbind(wald, score) }
ci(6,20,.95)

Our third approach to binomial inference follows the same idea of inverting a test statistic to construct a confidence interval but utilizes the likelihood ratio test (LRT) for the binomial parameter. Recall the likelihood function for \(Y\sim Bin(n,\pi)\):

\(\displaystyle \ell(\pi)= {n\choose y}\pi^y(1 - \pi)^{(n-y)}\)

The LRT statistic for \(H_0:\pi=\pi_0\) versus \(H_a:\pi\ne\pi_0\) is

\(\displaystyle G^2=2\log\dfrac{\ell(\hat{\pi})}{\ell(\pi_0)} = 2\left(y\log\dfrac{\hat{\pi}}{\pi_0}+(n-y)\log\dfrac{1-\hat{\pi}}{1-\pi_0}\right) \)

For large \(n\), \(G^2\) is approximately chi-square with one degree of freedom, and \(\pi_0\) will be rejected if \(G^2\ge \chi^2_{1,\alpha}\). Like the Wald and score test statistics, the LRT statistic is essentially a measure of disagreement between the sample estimate and the hypothesized value for \(\pi\). Larger values indicate more disagreement and more evidence to reject \(H_0\). And we can likewise construct a confidence interval as the set of all values of \(\pi_0\) that would not be rejected. Unfortunately, we must resort to numerical approximation for these limits. Here are the results for the smartphone data.

\( (0.1319, 0.5165) \)

Like the score interval the limits for the LRT interval are centered at a value closer to 0.5, compared with the Wald limits, which are centered at the MLE.

library('rootSolve')
ci = function(y,n,conf)
{   pi.hat = y/n
    z = qnorm(1-(1-conf)/2)
    wald = pi.hat+c(-1,1)*z*sqrt(pi.hat*(1-pi.hat)/n)
    score = (pi.hat+z^2/2/n+c(-1,1)*z*sqrt(pi.hat*(1-pi.hat)/n+z^2/4/n^2))/(1+z^2/n)
    loglik = function(p) 2*(y*log(pi.hat/p)+(n-y)*log((1-pi.hat)/(1-p)))-z^2
    lrt = uniroot.all(loglik,c(0.01,0.99))
    cbind(wald,score,lrt) }
ci(6,20,.95)	

Following up on our brief introduction to this extremely useful distribution, we go into more detail here in preparation for the goodness-of-fit test coming up. Recall that the multinomial distribution generalizes the binomial to accommodate more than two categories. For example, what if the respondents in a survey had three choices:

1. I feel optimistic.
2. I don't feel optimistic.
3. I'm not sure.

If we separately count the number of respondents answering each of these and collect them in a vector, we can use the multinomial distribution to model the behavior of this vector.

A goodness-of-fit test, in general, refers to measuring how well do the observed data correspond to the fitted (assumed) model. We will use this concept throughout the course as a way of checking the model fit. Like in linear regression, in essence, the goodness-of-fit test compares the observed values to the expected (fitted or predicted) values.

A goodness-of-fit statistic tests the following hypothesis:

\(H_0\colon\) the model \(M_0\) fits

vs.

\(H_A\colon\) the model \(M_0\) does not fit (or, some other model \(M_A\) fits)

Most often the observed data represent the fit of the saturated model, the most complex model possible with the given data. Thus, most often the alternative hypothesis \(\left(H_A\right)\) will represent the saturated model \(M_A\) which fits perfectly because each observation has a separate parameter. Later in the course, we will see that \(M_A\) could be a model other than the saturated one. Let us now consider the simplest example of the goodness-of-fit test with categorical data.

In the setting for one-way tables, we measure how well an observed variable X corresponds to a \(Mult\left(n, \pi\right)\) model for some vector of cell probabilities, \(\pi\). We will consider two cases:

1. when vector \(\pi\) is known, and
2. when vector \(\pi\) is unknown.

In other words, we assume that under the null hypothesis data come from a \(Mult\left(n, \pi\right)\) distribution, and we test whether that model fits against the fit of the saturated model. The rationale behind any model fitting is the assumption that a complex mechanism of data generation may be represented by a simpler model. The goodness-of-fit test is applied to corroborate our assumption.

Consider our dice example from Lesson 1. We want to test the hypothesis that there is an equal probability of six faces by comparing the observed frequencies to those expected under the assumed model: \(X \sim Multi(n = 30, \pi_0)\), where \(\pi_0=(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)\). We can think of this as simultaneously testing that the probability in each cell is being equal or not to a specified value:

\(H_0:\pi =\pi_0\)

where the alternative hypothesis is that any of these elements differ from the null value. There are two statistics available for this test.

An easy way to remember it is

\(X^2=\sum\limits_{j=1}^k \dfrac{(O_j-E_j)^2}{E_j}\)

where \(O_j = X_j\) is the observed count in cell \(j\), and \(E_j=E(X_j)=n\pi_{0j}\) is the expected count in cell \(j\) under the assumption that null hypothesis is true.

In some texts, \(G^2\) is also called the likelihood-ratio test (LRT) statistic, for comparing the loglikelihoods \(L_0\) and \(L_1\) of two models under \(H_0\) (reduced model) and \(H_A\) (full model), respectively:

Note that \(X^2\) and \(G^2\) are both functions of the observed data \(X\) and a vector of probabilities \(\pi_0\). For this reason, we will sometimes write them as \(X^2\left(x, \pi_0\right)\) and \(G^2\left(x, \pi_0\right)\), respectively; when there is no ambiguity, however, we will simply use \(X^2\) and \(G^2\). We will be dealing with these statistics throughout the course in the analysis of 2-way and \(k\)-way tables and when assessing the fit of log-linear and logistic regression models.

\(X^2\) and \(G^2\) both measure how closely the model, in this case \(Mult\left(n,\pi_0\right)\) "fits" the observed data. And both have an approximate chi-square distribution with \(k-1\) degrees of freedom when \(H_0\) is true. This allows us to use the chi-square distribution to find critical values and \(p\)-values for establishing statistical significance.

• If the sample proportions \(\hat{\pi}_j\) (i.e., saturated model) are exactly equal to the model's \(\pi_{0j}\) for cells \(j = 1, 2, \dots, k,\) then \(O_j = E_j\) for all \(j\), and both \(X^2\) and \(G^2\) will be zero. That is, the model fits perfectly.

• If the sample proportions \(\hat{\pi}_j\) deviate from the \(\pi_{0j}\)s, then \(X^2\) and \(G^2\) are both positive. Large values of \(X^2\) and \(G^2\) mean that the data do not agree well with the assumed/proposed model \(M_0\).

How large is the discrepancy between the two proposed models? The previous analysis provides a summary of the overall difference between them, but if we want to know more specifically where these differences are coming from, cell-specific residuals can be inspected for relevant clues. There are two types of residuals we will consider: Pearson and deviance residuals.

is called the Pearson residual for cell \(j\), and it compares the observed with the expected counts. The sign (positive or negative) indicates whether the observed frequency in cell \(j\) is higher or lower than the value implied under the null model, and the magnitude indicates the degree of departure. When data do not fit the null model, examination of the Pearson residuals often helps to diagnose where the model has failed.

How large should a "typical" value of \(r_j\) be? Recall that the expectation of a \(\chi^2\) random variable is its degrees of freedom. Thus, if a model is correct, \(E(X^2) \approx E(G^2) \approx k − 1\), and the typical size of a single \(r_{j}^{2}\) is \(\dfrac{k − 1}{k}\). Thus, if the absolute value, \(|r_j|\), is much larger than \(\sqrt{(k-1)/k}\)—say, 2.0 or more—then the model does not appear to fit well for cell \(j\).

Although not as intuitively as the \(X^2\) statistic, the deviance statistic \(G^2=\sum\limits_{j=1}^k d^2_j\) can be regarded as the sum of squared deviance residuals,

\(d_j=\sqrt{\left|2X_j\log\dfrac{X_j}{n\pi_{0j}}\right|}\times \text{sign}(X_j-n\pi_{0j})\)

The sign function can take three values:

• -1 if \((X_j - n\pi_{0j} ) < 0\),
• 0 if \((X_j- n\pi_{0j} ) = 0\), or
• 1 if \((X_j- n\pi_{0j}) > 0\).

When the expected counts \(n\pi_{0j}\) are all fairly large (much greater than 5) the deviance and Pearson residuals resemble each other quite closely.

If an \(X_j\) is zero and all \(\pi_{0j}\)s are positive, then the Pearson \(X^2\) can be calculated without any problems, but there is a problem in computing the deviance, \(G^2\); if \(X_j = 0\) then the deviance residual is undefined, and if we use the standard formula,

\(G^2=2\sum\limits_{j=1}^k X_j\log\dfrac{X_j}{n\pi_{0j}}\)

an error will result. But if we write the deviance as

\(G^2=2\log\dfrac{L(\pi_0;X)}{L(\hat{\pi};X)}=2\log\prod\limits_{j=1}^k \left(\dfrac{X_j/n}{\pi_{0j}}\right)^{X_j}\)

Now, a cell with \(X_j=0\) contributes 1 to the product and may be ignored. Thus, we may calculate the deviance statistic as

\(G^2=2\sum\limits_{i:X_j>0} X_j\log\dfrac{X_j}{n\pi_{0j}}\)

Alternatively, we can set the deviance residuals to zero for cells with \(X_j=0\) and take \(G^2= \sum_j d_j^2\) as before. But if we do that, \(d_j = 0\) should not be interpreted as "the model fits well in cell \(j\)". The fit could be quite poor, especially if \(E_j\) is large.

If any element of vector \(\pi_0\) is zero, then \(X^2\) and \(G^2\) both break down. Simple ways to avoid such problems include putting a very small mass in all the cells, including the zero-cell(s) (e.g., \(0.05\)) or removing or combining cells and re-doing the tests.

Structural zeros and sampling zeros

Note that an observed zero in a cell does not necessarily mean that there cannot be any observation in that cell. A sampling zero is an observed zero count in a cell for which a positive value is possible. An example of this would have occurred in our die experiment if by chance the number 1 had not shown up in any of the 30 rolls. Alternatively, consider an example of categorizing male and female patients on whether they have carcinoma. Since only females can have ovarian cancer, and only males can have prostate cancer, the following cells are said to have structural zeros, and the underlying cell probabilities \(\pi_j=0\). To handle structural zeros, an adjustment to the degrees-of-freedom is required.

For many statistical models, we do not know the vector of probabilities \(\pi\) a priori, but can only specify it up to some unknown parameters. More specifically, the cell proportions may be known functions of one or more other unknown parameters. For example, suppose that a gene is either dominant (A) or recessive (a), and the overall proportion of dominant genes in the population is \(p\). If we assume mating is random (i.e. members of the population choose their mates in a manner that is completely unrelated to this gene), then the three possible genotypes—AA, Aa, and aa—should occur in the so-called Hardy-Weinberg proportions:

This is equivalent to saying that the number of dominant genes that an individual has (0, 1, or 2) is distributed as \(Bin(2,p)\), where the parameter \(p\) is not specified. We have to first estimate it in order to obtain the expected counts under this model.

In such an example, the null hypothesis specifies some constraint on the parameters but stops short of providing their values specifically. In more general notation, the model specifies that:

where \(g_1, \ldots, g_k\) are known functions, but the parameter \(\theta\) is unknown (e.g., \(p\) in the example above). Let \(S_0\) denote the set of all \(\pi\) that satisfy these constraints for some parameter \(\theta\). We want to test

\(H_0 \colon \pi \in S_0\) versus \(H_A \colon \pi \in S\)

where \(S\) denotes the probability simplex (the space) of all possible values of \(\pi\). (Notice that \(S\) is a \((k-1)\)-dimensional space, but the dimension of \(S_0\) is the number of free parameters in \(\theta\).) The method for conducting this test is as follows:

1. Estimate \(\theta\) by an efficient method (e.g., maximum likelihood), and denote the estimate by \(\hat{\theta}\).
2. Calculate (estimated) expected counts under the null hypothesis by \(n\hat{\pi}_{0j}\), where \(\hat{\pi}_{0j} = g_j(\hat{\theta})\), for \(j=0,\ldots,k\). We may still denote these as \(E_j\), even though they involve an estimated parameter.
3. Calculate the goodness-of-fit statistics \(X^2(x,\hat{\pi}_0)\) and \(G^2(x,\hat{\pi}_0)\) using the usual formulas from the \(O_j\) and the \(E_j\).

That is, once the \(E_j\) are obtained from the estimated parameter(s), we use the usual formulas:

\(X^2=\sum\limits_j \dfrac{(O_j-E_j)^2}{E_j}\) and \(G^2=2\sum\limits_j O_j \log\dfrac{O_j}{E_j}\)

If \(X^2\) and \(G^2\) are calculated as described above, then under the null hypothesis, both are approximately  \(\chi^{2}_{\nu}\), where \(\nu\) equals the number of unknown parameters under the alternative hypothesis minus the number of unknown parameters under the null hypothesis, \(\nu = (k − 1) − d\) where \(d = dim(\theta)\), i.e., the number of parameters that required estimation.

In this lesson, we made a distinction among three types of large-sample hypothesis tests and related confidence intervals used in the analysis of categorical data: Wald, Likelihood-Ratio, and Score. The Wald test is the most widely used one. For example, we will see when we fit a logistic regression model that the z-statistic and the confidence intervals for the regression parameter estimates are Wald CIs.

Lesson 2 also introduced the goodness-of-fit test and the corresponding residuals. The distinction was made between the goodness-of-fit tests with known versus unknown population parameters. This idea of comparing the assumed model (H₀) to a distribution of the observed data, and assessing how close they fit, will be used throughout the course, as will these test statistics. For example, in the next lesson, we will see that the chi-square test of independence is just a goodness-of-fit test where the assumed model (H₀) of the independence of two random variables is compared to a saturated model, that is to the observed data.