---
title: "Stat 5421 Notes: Review of Intro Stats with Some Extras"
author: "Charles J. Geyer"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
    html_document:
        md_extensions: -tex_math_single_backslash
        mathjax: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML
        number_sections: true
        toc: true
    pdf_document:
        number_sections: true
        latex_engine: xelatex
        toc: true
---

# License

This work is licensed under a Creative Commons
Attribution-ShareAlike 4.0 International License
(http://creativecommons.org/licenses/by-sa/4.0/).

# What is This?

Chapter 1 of the textbook (Agresti) already reviews bits
of the theory of probability and statistics,
but [IMHO](https://en.wiktionary.org/wiki/IMHO) that chapter is not basic
enough for some students.  Hence this.

This course is about analysis of discrete data, mostly count data.
So most of intro stats is not relevant.  This course also has
as a prerequisite some course beyond intro stats.  So we assume
you are all more sophisticated than the typical student in a intro stats
course.  But just having had some course (a long list of options, no
particular course required) doesn't
mean you are familiar with any particular material not in intro stats.

So we start by reviewing the material in intro stats that is in
what this course is about.
Along the way we add a few alternative methods that intro stats omits,
which will be prominent in this course.

# Statistical Inference

It may be annoying that much of the content of this course can be
summarized as *it's complicated.*  In order to keep all of these complications
organized, we will make a lot of lists dividing up these complications.

*Statistical inference* is the process of saying something about the
probability distribution that describes the data.  [More on probability
distributions later](ch1.html).

Statistical inference comes in two main kinds and [several minor
kinds](https://en.wikipedia.org/wiki/Statistical_inference).
The major kinds are usually called

 * frequentist and

 * Bayesian.

The minor kinds will not be used in this course.  The two kinds above
may not have been mentioned in your previous statistics courses.  If
not, everything they did was frequentist, and they just ignored Bayesian.
Hence you are expected to know something about frequentist, but may not
know anything about Bayesian.

All of statistical inference has the problem that it can never get the
right answer: estimators are equal to the parameters they estimate with
very low probability.
Hence every statistical inference comes with an indication
of how wrong it may be.  There are three main kinds of frequentist
statistical inference

 * point estimates,

 * confidence intervals (sometimes called interval estimates), and

 * tests of statistical hypotheses (hypothesis tests, for short, also
   called significance tests).

You have seen all of these in previous statistics courses.

There is only one kind of Bayesian inference.  Bayesians say that
probability theory is the correct way to describe uncertainty.
So a probability distribution is the correct way to describe our uncertainty
about the true unknown distribution of the data or any of its parameters.
We are uncertain about the true unknown parameter value both before and after
we see the data.
But our uncertainties before and after are different.
We know more, are less uncertain, after the data have been seen.
The distribution that describes our uncertainty before is called
the *prior* distribution (or just prior),
the distribution that describes our uncertainty after is called
the *posterior* distribution (or just posterior).
The mathematical operation that connects prior, posterior, and data
is called *Bayes' rule* or *Bayes' theorem*.  It involves calculus
and so will not be explicit in this course.  You will learn all about
that when you take a theory course.
We will have a whole unit on [Bayesian computation](mcmc-bayes.html)
but very little theory.  A taste of Bayesian inference is given
in the [section on Bayesian inference about a
population proportion](#bayesian-inference) below.

# Frequentist Inference

## Point Estimates

Point estimates are often obvious, especially in intro stats.

 * The sample mean is a natural estimate of the population mean.

 * The sample median is a natural estimate of the population median.

 * The sample proportion is a natural estimate of the population proportion.

And so forth.  But often point estimates, especially in complicated models
with multiple parameters,

 * have no formula that expresses them ― they can only be obtained by
   using numerical optimization with a computer ― or

 * any formula that expresses them is difficult to derive (uses calculus)
   and is not at all intuitive or obvious.

Examples of these [start with these notes](likelihood-compute.html).

The most important thing to remember about point estimates is

 * *point estimates are not the parameters they estimate* or

 * $\hat{\theta}$ is not $\theta$.

The whole of statistical theory is about dealing with this issue.
Every statistical estimate is wrong.  And we have to quantify how wrong.

## Hypothesis Tests

Hypothesis tests compare two statistical models (families of probability
distributions).  Formally they are called the null hypothesis and the
alternative hypothesis.  Tests come in three kinds:

 * tests about the value of one parameter, which are further subdivided into
     - two-tailed tests,
     - lower-tailed tests, and
     - upper-tailed tests,

 * tests comparing nested statistical models (which includes the bullet
   point above), and

 * none of the above.

We will not deal with the last except to say that, in principle,
a hypothesis test can compare any two hypotheses (logical statements)
about the true unknown distribution of the data.  But the only tests
of practical interest are the first two kinds, and those are the only
ones covered in statistics courses, including this one.  So we will
not mention "none of the above" again except when discussing
[model selection and model averaging](glmbb.html), which is not,
strictly speaking, hypothesis testing, but rather a competitor.

## Intervals

Confidence intervals are only about one parameter.
They do generalize to the notion of *confidence regions* for parameter
vectors, but those are rarely used.  Confidence intervals also come in two
kinds:

 * two-sided and

 * one-sided,

but the latter are rarely used and not covered in most statistics courses.
But we will cover them at the end of the course in the [unit about inference
when data are on the boundary](infinity.html).

## Duality of Tests and Intervals

Hypothesis tests and confidence intervals are two ways of looking at the
same math.  Hypothesis tests can be used to make confidence intervals and
vice versa.

 * If you have a recipe for doing hypothesis tests at
   significance level $\alpha$ and for any null hypothesis $\theta$,
   then a confidence interval having coverage $1 - \alpha$ is the
   set of $\theta$ that the hypothesis test does not reject.

 * Conversely, if you have a recipe for doing a confidence interval
   with coverage $1 - \alpha$, then a hypothesis test with null hypothesis
   $\theta$ and significance level $\alpha$ accepts the null hypothesis
   if $\theta$ is in the confidence interval
   and rejects it otherwise.

In more detail, two-tailed tests are dual to two-sided intervals.  And
two-tailed tests have null and alternative hypotheses
\begin{align*}
   H_0 : \text{true unknown parameter value} = \theta
   \\
   H_1 : \text{true unknown parameter value} \ne \theta
\end{align*}
So to make a confidence interval we (conceptually) perform the test
for each $\theta$ in the parameter space, and the interval is the
set of $\theta$ that the test accepts.

Thus, in theory, we have to do an infinite number of tests,
but, in practice, we don't.  We can just manipulate some formulas,
or at worst solve a few equations numerically.
Examples in the [next section](#binomial-confint).

Going the other way is much simpler.  The hypothesis test rejects
the null hypothesis of equality to $\theta$ if and only if $\theta$ is not
in the confidence interval.

Also in more detail, one-tailed tests are dual to one-sided intervals.  A
lower-tailed test has null and alternative hypotheses
\begin{align*}
   H_0 : \text{true unknown parameter value} = \theta
   \\
   H_1 : \text{true unknown parameter value} < \theta
\end{align*}
So to make a confidence interval we (conceptually) perform the test
for each $\theta$ in the parameter space, and the interval is the
set of $\theta$ that the test accepts.

Usually, if the test accepts a particular value of $\theta$, then it will
also accept any value of $\theta$ greater than that.
Thus a lower-tailed test
is dual to a lower-bound confidence interval (the interval has the form
$(l, \infty)$ for some lower bound $l$).
Similarly, an upper-tailed test
is dual to an upper-bound confidence interval (the interval has the form
$(- \infty, u)$ for some upper bound $u$).  If the parameter space
$\Theta$ is not infinite in both directions, then

 * a lower-bound interval has the form
   $\{\, \theta \in \Theta : l < \theta \,\}$ and

 * an upper-bound interval has the form
   $\{\, \theta \in \Theta : \theta < u \,\}$.

If you are happy with what this section has said so far, you may skip
to the next section, but if the description of null and alternative hypotheses
for a one-tailed test is not what you saw in other courses, read on.

Some textbooks say the null and alternative hypotheses for a lower-tailed
test are
\begin{align*}
   H_0 : \text{true unknown parameter value} \ge \theta
   \\
   H_1 : \text{true unknown parameter value} < \theta
\end{align*}
(we have changed $=$ to $\ge$ in the null hypothesis).  Usually, these
hypotheses give rise to exactly the same test procedure as the hypotheses
given above.  If they do not, then these hypotheses are a really dumb idea
and one has no idea how to do the test.

The theoretical considerations
centering around hypotheses of this form are [discussed in these theoretical
statistics lecture notes,
pp. 181–185](https://www.stat.umn.edu/geyer/5102/slides/s2.pdf#page=181),
but there is no need to look at that unless you are really curious about
this issue.

## Multiple Tests or Intervals

The mathematics of frequentist statistics is based on the following dogma.

 * Do only one hypothesis test or confidence interval.

 * Choose the procedure to be done before the data are collected.

 * Do it, and report it.

Anything else is merely exploratory.

Doing multiple procedures and cherry picking the results you like is
scientific fraud if you do not admit it and a red flag that your results
are meaningless if you do admit it.

You can rigorously do multiple hypothesis tests or confidence intervals if you
[make adjustments to your procedures to account for
that](likelihood-compute.html#simultaneous-coverage), which you may have
seen in courses beyond intro stats and which we will cover in this
course.

# Inference About a Population Proportion {#proportion}

The binomial distribution may or may not have been explicitly described
in the courses you have had previously.  Many intro stats courses
discuss this under the description *inference about a population proportion*.
But all advanced courses, including this one, say *binomial distribution*.
We will cover [inference for the binomial distribution](binom.html)
in exhaustive detail later.  For now we just review what was said about
it in intro stats, and we add a few alternative methods that were not
covered in intro stats.

Let $\pi$ denote the population proportion of some quantity and
$\hat{\pi}_n$ denote the sample proportion for sample size $n$.
Then it is an instance of the central limit theorem that
\begin{equation} \tag{$*$}
   \frac{\hat{\pi}_n - \pi}{\sqrt{\pi (1 - \pi) / n}}
\end{equation}
has an approximate standard normal distribution (mean zero, variance one)
when $n$ is sufficiently
large.  How large it needs to be depends on what the true unknown parameter
value $\pi$ is, what question you are asking, and how accurate you want
the approximation to be.  The reason why this is an instance of the central
limit theorem is that $\hat{\pi}_n$ is a sample mean.

Here $\pi$ is the unknown parameter, not $3.1415926\ldots$.
In this we are following Agresti who follows the rule that parameters
are always denoted by Greek letters (not mostly Greek letters,
like some textbooks do).

In this course, there is no exact inference (except for [fuzzy tests
and confidence intervals](fuzzy.html))
unlike what you have seen in courses about
analyzing continuous data.  For normally distributed data, you have
t and F tests, which are exact (assuming exact homoscedastic normality).
But there is nothing like that for categorical data.  We have to make
do with approximations, and this starts here.

*If you are trying to do a t or F test in this course, then you are
making a mistake.*

## Hypothesis Tests

### The Usual Procedure

The generally accepted way to do a hypothesis test about a sample proportion
is to use ($*$) as the test statistic with $\pi$ being the hypothesized value
under the null hypothesis rather than the true unknown parameter value.

Suppose we want to test the null hypothesis $\pi = 1/2$ and the data
are as follows.
```{r binom-test}
pi <- 1 / 2
n <- 10
x <- 8
pi.hat <- x / n
z <- (pi.hat - pi) / sqrt(pi * (1 - pi) / n)
z
# upper-tailed p-value
pnorm(z, lower.tail = FALSE)
prop.test(x, n, p = pi, alternative = "greater", correct = FALSE)$p.value
# lower-tailed p-value
pnorm(z)
prop.test(x, n, p = pi, alternative = "less", correct = FALSE)$p.value
# two-tailed p-value
pnorm(- abs(z)) + pnorm(abs(z), lower.tail = FALSE)
2 * pnorm(abs(z), lower.tail = FALSE)
2 * pnorm(- abs(z))
2 * min(pnorm(z), pnorm(z, lower.tail = FALSE))
prop.test(x, n, p = pi, correct = FALSE)$p.value
```
(The last four results involving `pnorm` are equal because of symmetry of
the normal distribution.  The `correct = FALSE` tells `prop.test` not to
use its default correction for continuity, which few, if any,
textbooks or experts recommend.)

The interpretation of the test is that, if we can believe an upper-tailed
test is valid ([if we can believe that this procedure was selected as the
only hypothesis to be tested
before the data were collected](#multiple-tests-or-intervals) — or at least
could have been ― they weren't data snooping),
then $P = `r pnorm(- z) |> round(3)`$ is below the conventional 0.05
level for "statistical significance".
But you might consider this weak, moderate, or strong evidence against
the null hypothesis, depending on what the data
are about, details of the experiment, and so forth.

If the authors of the test (or our second-guessing of them) consider
the two-tailed test appropriate,
then $P = `r 2 * pnorm(- z) |> round(3)`$ is above the conventional 0.05
level for "statistical significance".
But you might still consider this weak or moderate evidence against
the null hypothesis, again depending on what the data
are about, details of the experiment, and so forth.

Later on we will see that this test is a special case of a general
procedure called the Rao test or the score test.

```{r save.rao.test, echo=FALSE}
save.rao.p.value <- 2 * pnorm(- abs(z))
save.rao.p.value <- formatC(save.rao.p.value, format = "f", digits = 3)
```

### The Likelihood Ratio Test {#binomial-lrt}

We interrupt our review of intro stats to show you something you
may not have seen before.  This procedure is very widely used, and we
will use it a lot in this course.  You may have used it yourself.  It
is what R generic function `anova` does when handed the results of
calls to R function `glm`.  An application of large sample theory
of statistics (the central limit theorem plus Taylor series approximation)
says that
\begin{equation} \tag{${*}{*}$}
   2 \left[ x \log\left(\frac{\hat{\pi}_n}{\pi}\right)
   +
   (n - x) \log\left(\frac{1 - \hat{\pi}_n}{1 - \pi}\right) \right]
\end{equation}
has an approximate chi-square distribution distribution with one degree of
freedom.

So another good way to do a hypothesis test about a sample proportion
is to use (${*}{*}$) as the test statistic with $\pi$ being the hypothesized
value under the null hypothesis rather than the true unknown parameter value.

You might ask: who ordered that?  Where does (${*}{*}$) come from?
The answer is from theory.  [More on this later](likelihood.html).

Since this test statistic is large when $\hat{\pi}_n$ is far from $\pi$
in either direction (as shown by a
[plot below](#likelihood-intervals)), this is inherently a two-tailed test.
One-tailed tests are covered in the
[next section](#binomial-signed-lrt).
```{r binom-lrt}
tstat <- 2 * (x * log(pi.hat / pi) + (n - x) * log((1 - pi.hat) / (1 - pi)))
# two-tailed hypothesis test
pchisq(tstat, df = 1, lower.tail = FALSE)
```

```{r save.wilks.test, echo=FALSE}
save.wilks.p.value <- pchisq(tstat, df = 1, lower.tail = FALSE)
save.wilks.p.value <- formatC(save.wilks.p.value, format = "f", digits = 3)
```

These tests are asymptotically equivalent, meaning that
for sufficiently large $n$ they give nearly the same results for the same data.
Here they give somewhat different results,
$P = `r 2 * pnorm(- z) |> round(3)`$ for the usual (Rao) test and
$P = `r pchisq(tstat, df = 1, lower.tail = FALSE) |> formatC(format = "f",
digits = 3)`$ for
the likelihood ratio test.

**Warning:** the likelihood ratio test statistic calculation above blows up,
giving `0 * -Inf = NaN` when either the number of successes $x$ or
the number of failures $n - x$ is zero, but
this is incorrect, just bad calculation.
Since
$$
   x \log(x) \to 0, \qquad \text{as $x \to 0$}
$$
the term in the sum having observed sucesses zero or observed failures
zero should be zero, so to be careful we calculate
```{r binom-lrt-careful}
tstat <- 0
if (x > 0) tstat <- tstat + 2 * x * log(pi.hat / pi)
if (x < n) tstat <- tstat + 2 * (n - x) * log((1 - pi.hat) / (1 - pi))
# two-tailed hypothesis test
pchisq(tstat, df = 1, lower.tail = FALSE)
```

Let's see what happens when we make the sample size larger
```{r binom-both}
n <- 100
x <- 61
pi.hat <- x / n
z <- (pi.hat - pi) / sqrt(pi * (1 - pi) / n)
tstat <- 2 * (x * log(pi.hat / pi) + (n - x) * log((1 - pi.hat) / (1 - pi)))
# two-tailed p-value for usual (Rao) test
2 * pnorm(- abs(z))
# two-tailed p-value for likelihood ratio (Wilks) test
pchisq(tstat, df = 1, lower.tail = FALSE)
```

Closer.  And the $P$-values get closer and closer as $n$ goes to infinity.

A crazy way to calculate the results above.  We do this to show how these
are special cases of the way R compares complicated models
(general tests of model comparison done by R generic function `anova`)
```{r anova}
gout0 <- glm(cbind(x, n - x) ~ 0, family = binomial)
gout1 <- glm(cbind(x, n - x) ~ 1, family = binomial)
anova(gout0, gout1, test = "Rao")[2, "Pr(>Chi)"]
anova(gout0, gout1, test = "LRT")[2, "Pr(>Chi)"]
```

These look odd, but that is just he way R function `glm` works.
For binomial data,

 * the response is either a zero-or-one-valued vector, or a two-column
   matrix of counts, first column successes, second column failures,
   so here we have a $1 \times 2$ matrix `cbind(x, n - x)`,

 * in the formula 0 means no intercept and 1 means has intercept
   (which is the default, but we need the 1 here because there are no
   other terms in the formula),

 * no intercept just happens to be the null hypothesis we want but
   we could test a different null hypothesis using the `offset`
   argument to R function `glm`.
  

### The Signed Likelihood Ratio Test {#binomial-signed-lrt}

Since the square of a standard normal random variable has the chi-squared
distribution with one degree of freedom (by definition), the square root
of a chi-square random variable with one degree of freedom has
the same distribution as the absolute value of a standard normal
random variable.  If we attach a random sign, plus and minus being
equally probable, to that square root, we get back a standard normal
random variable.  By ($*$) being approximately standard normal,
$\hat{\pi}_n - \pi$ is approximately equally likely to be plus or minus
(when $\pi$ is the true unknown parameter value).
Thus the square root of (${*}{*}$) with the sign of $\hat{\pi}_n - \pi$
attached is approximately standard normal for large $n$.

Thus we have another test statistic
\begin{equation} \tag{${*}{*}{*}$}
   \mathop{\rm sign}(\hat{\pi}_n - \pi)
   \sqrt{
   2 x \log\left(\frac{\hat{\pi}_n}{\pi}\right)
   +
   2 (n - x) \log\left(\frac{1 - \hat{\pi}_n}{1 - \pi}\right)
   }
\end{equation}
which if squared gives back the test statistic (${*}{*}$) of the
likelihood ratio test.  But now this test statistic has an approximately
normal distribution and can be used for one-tailed tests.

Because the sign of the test statistic is the sign of $\hat{\pi}_n - \pi$,
a lower-tailed test has a low $P$-value when $\hat{\pi}_n$ is a lot
less than $\pi$, and
an upper-tailed test has a low $P$-value when $\hat{\pi}_n$ is a lot
greater than $\pi$, which is what these tests are supposed to do.

We return to our original data
```{r signed-lrt}
n <- 10
x <- 8
pi.hat <- x / n
tstat <- sign(pi.hat - pi) * sqrt(2 * (x * log(pi.hat / pi) +
    (n - x) * log((1 - pi.hat) / (1 - pi))))
# upper-tailed p-value
pnorm(tstat, lower.tail = FALSE)
# lower-tailed p-value
pnorm(tstat)
# two-tailed p-value
2 * pnorm(- abs(tstat))
```

We included the two-tailed $P$-value only to show that a two-tailed
signed likelihood ratio test gives exactly the same results as a
likelihood ratio test.

### The Wald Test

This section has no example code because this procedure is not
recommended for hypothesis tests about a population proportion.
We only have this section to mark that there is yet another general
test procedure to learn.

The Wald test, if we did it (which we won't), would use the test
statistic ($\dagger$) in the following section but would otherwise
calculate $P$-values just like the Rao test and the signed likelihood
ratio test.

But that is not recommended by textbooks.  The justification for using
($*$) instead of ($\dagger$) is that the $P$-value is supposed to be
calculated using the assumption that $\pi$ is the true unknown parameter
value, that is, the true unknown parameter value is the value hypothesized
by the null hypothesis.  Under that assumption the standard deviation of
$\hat{\pi}_n - \pi$ is $\sqrt{\pi (1 - \pi) / n}$.  Hence that should
be used to standardize it, giving ($*$) rather than ($\dagger$).

Later on we will see that this test is a special case of a general
procedure called the Wald test.
When $n$ is large there is nothing wrong with Wald tests.
If you have ever looked at the output of R function `summary`
applied to the output of R function `lm` or `glm`, then you have
looked at the $P$-values in the rightmost column.  All come from Wald tests.

So if there were anything really wrong with Wald tests, there would
be a big problem.  But there isn't (when $n$ is large).

Here, of course, $n = `r n`$ is not "large".

## Confidence Intervals {#binomial-confint}

### Usual Intervals

We now return to reviewing intro stats.
Another application of large sample theory of statistics
(the central limit theorem plus Slutsky's theorem, also called the
plug-in principle), says that
\begin{equation} \tag{$\dagger$}
   \frac{\hat{\pi}_n - \pi}{\sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}}
\end{equation}
has an approximate standard normal distribution
when $n$ is sufficiently large
(the only difference between ($*$) and ($\dagger$) is that the
former has $\pi$ in the denominator where the latter has $\hat{\pi}_n$).

The concern that the Wald test does not use the null hypothesis in calculating
standard deviation, which leads textbooks to (sometimes) avoid recommending it,
does not apply to confidence intervals, because they have no null hypothesis.

We find the confidence interval dual to this test by inverting the test.
A two-tailed test rejects the null hypothesis at level $\alpha$ when
the absolute value of ($\dagger$) is large.  So we accept $\pi$ when
$$
   \left\lvert
   \frac{\hat{\pi}_n - \pi}{\sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}}
   \right\rvert < z_{\alpha / 2}
$$
where $z_{\alpha / 2}$ denotes the upper $\alpha / 2$ quantile of the
standard normal distribution.  The set of all accepted $\pi$ is the
confidence interval with (approximate, large $n$) coverage $1 - \alpha$.

Solving
$$
   - z_{\alpha / 2} <
   \frac{\hat{\pi}_n - \pi}{\sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}}
   < z_{\alpha / 2}
$$
for $\pi$ gives
$$
   \hat{\pi}_n - z_{\alpha / 2} \sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}
   < \pi <
   \hat{\pi}_n + z_{\alpha / 2} \sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}
$$
so we get a confidence interval with coverage $1 - \alpha$ having endpoints
$$
   \hat{\pi}_n \pm z_{\alpha / 2} \sqrt{\hat{\pi}_n (1 - \hat{\pi}_n) / n}
$$

Try that.
```{r binom-wald-interval}
conf.level <- 0.95
crit <- qnorm((1 + conf.level) / 2)
crit
se.pi.hat <- sqrt(pi.hat * (1 - pi.hat) / n)
pi.hat + c(-1, 1) * crit * se.pi.hat
```
(the `c(-1, 1)` is R for $\pm$).

This interval is called a Wald interval, like the test it is dual to
being called a Wald test.

Note that the upper endpoint of the confidence interval is outside
the parameter space.  But it means the same as if we had reported the
confidence interval as
(`r (pi.hat - crit * se.pi.hat) |> round(3)`, 1)

In addition to this minor blemish, this Wald interval
for a population proportion
has a very big problem.  When $\hat{\pi}_n$ is equal to zero or one
(when $x$ is equal to 0 or $n$) the interval has width zero, which
is useless and stupid.
But this happens with negligible probability when $n$ is large.

```{r save.wald.interval, echo=FALSE}
save.wald.interval <- pi.hat + c(-1, 1) * crit * se.pi.hat
save.wald.interval <- formatC(save.wald.interval, format = "f", digits = 3)
```

We continue this analysis with
a [section on fixing up Wald intervals](#fixup-wald) below.

### Score Intervals

Thus some textbooks recommend an interval that is harder to calculate.
This inverts the test that uses ($*$) as the test statistic.  This is
harder.  The interval is the set of $\pi$ that satisfy
$$
   \left\lvert
   \frac{\hat{\pi}_n - \pi}{\sqrt{\pi (1 - \pi) / n}}
   \right\rvert < z_{\alpha / 2}
$$
It is not obvious how to solve this, but square both sides to get
rid of the absolute value giving
$$
   \frac{(\hat{\pi}_n - \pi)^2}{\pi (1 - \pi) / n} < z_{\alpha / 2}^2
$$
and multiply both sides by the left-hand side denominator giving
$$
   (\hat{\pi}_n - \pi)^2 < z_{\alpha / 2}^2 \pi (1 - \pi) / n
$$
Since both sides are quadratic in $\pi$, this gives us a quadratic
equation to solve.  We won't bother with the details (they are in
[these theory notes,
slides 113—116](https://www.stat.umn.edu/geyer/5102/slides/s2.pdf#page=113)).
The endpoints of this interval are
$$
   \frac{\hat{\pi}_n + \frac{z_{\alpha / 2}^2}{2 n}
   \pm z_{\alpha / 2}
   \sqrt{ \frac{\hat{\pi}_n (1 - \hat{\pi}_n)}{n} +
   \frac{z_{\alpha / 2}^2}{4 n^2}}}
   {\left(1 + \frac{z_{\alpha / 2}^2}{n}\right)}
$$
What a mess!  But it does have some advantages over the Wald interval.

 * It never has zero width.

 * It is always contained in the parameter space.

 * It is asymptotically equivalent to (no better, no worse, when $n$ is large)
   the Wald interval.

Try it.
```{r binom-score}
prop.test(x, n, conf.level = conf.level, correct = FALSE)$conf.int
(pi.hat + crit^2 / (2 * n) + c(-1, 1) * crit *
    sqrt(pi.hat * (1 - pi.hat) / n + crit^2 / (4 * n^2))) / (1 + crit^2 / n)
```
We used the complicated formula just to show that R function `prop.test`
with argument `correct = FALSE` does indeed calculate this interval.

Like the hypothesis test this is dual to, this interval is called
the score (or Rao) interval.  It is sometimes also called the Wilson
interval after the first person to propose this special case of the
general concept.

Despite the complexity of this interval, some intro stats textbooks
recommend it because of the bad performance (in some cases when $n$ is
small) of the Wald interval.

```{r save.rao.interval, echo=FALSE}
save.rao.interval <-
    prop.test(x, n, conf.level = conf.level, correct = FALSE)$conf.int
save.rao.interval <- formatC(save.rao.interval, format = "f", digits = 3)
```

### Likelihood Intervals

No intro stats books that I know of recommend these intervals,
but they are widely used and we will use them in this course.
They are obtained by inverting the test having test statistic (${*}{*}$).
The interval is the set of $\pi$ that satisfy
$$
   2 \left\lvert
   x \log\left(\frac{\hat{\pi}_n}{\pi}\right)
   +
   (n - x) \log\left(\frac{1 - \hat{\pi}_n}{1 - \pi}\right)
   \right\rvert < z_{\alpha / 2}^2
$$
(we could have written $\chi^2_\alpha$ instead of $z_{\alpha / 2}^2$
where the former means the upper $\alpha$ quantile of the chi-squared
distribution with one degree of freedom, because squaring a standard
normal random variable to get a chi-squared random variable maps two
tails onto one tail).  But there is no way to solve this in closed form.
We have to do it by numerical solution of equations ([more on this
later](binom.html#likelihood-interval)).

```{r binom-likelihood-interval}
invlogit <- function(theta) 1 / (1 + exp(- theta))
library(MASS)
gout <- glm(cbind(x, n - x) ~ 1, family = binomial)
suppressMessages(confint(gout, level = conf.level)) |> invlogit()
```

```{r save.wilks.interval, echo=FALSE}
save.wilks.interval <-
    suppressMessages(confint(gout, level = conf.level)) |> invlogit()
save.wilks.interval <- formatC(save.wilks.interval, format = "f", digits = 3)
```

This code is rather bizarre.  The R function that does this,
method `glm` of R generic function `confint` in R package `MASS`
only does likelihood intervals for regression coefficients of
generalized linear model fits so we have to fake it out by

 * treating this problem as a logistic regression with no covariates
   (only an intercept),

 * letting it make a likelihood interval for that intercept parameter,

 * mapping this interval from the logit scale to the probability scale
   ([more on logit in the next section](#outside)).

Later we will see how to do this by a more direct but still complicated
method
([here](binom.html#likelihood-interval) and
 [here](likelihood.html#likelihood-interval) and
 [here](likelihood-compute.html#likelihood-intervals) and
 [here](likelihood-compute.html#simultaneous-coverage)).

What we said about score intervals, also applies to likelihood intervals.

 * It never has zero width.

 * It is always contained in the parameter space.

 * It is asymptotically equivalent to (no better, no worse, when $n$ is large)
   the Wald interval or the score interval.

So all three intervals are equally good for sufficiently large $n$
and the score and likelihood intervals are equally good for all $n$
(although they only have close to their nominal coverage for large $n$).

Since the computing here was completely mysterious, we show that it worked
by plotting (${*}{*}$) as a function of $\pi$ and showing that we have
correctly computed the interval.
```{r label="show-lrt",fig.align="center",fig.cap="Likelihood Ratio Test Statistic Function and Likelihood Interval.  Solid curve is (${*}{*}$) plotted as function of $\\pi$.  Horizontal dashed line is at chi-squared critical value.  Vertical dashed lines are at the endpoints of the calculated likelihood interval."}
foo <- function(pi) 2 * (x * log(pi.hat / pi) +
    (n - x) * log((1 - pi.hat) / (1 - pi)))
crit <- qchisq(conf.level, df = 1)
curve(foo, from = 0.4, to = 1, n = 1001, ylim = c(0, 1.5 * crit),
    ylab = expression(LRT(pi)), xlab = expression(pi))
abline(h = crit, lty = "dashed")
foo <- suppressMessages(confint(gout, level = conf.level)) |> invlogit()
abline(v = foo, lty = "dashed")
```

### Fixing Up the Wald Interval {#fixup-wald}

#### Fixing Up Interval Goes Outside Parameter Space {#outside}

The trick used in the preceding section, where we made a confidence interval
for one parameter and then mapped that to a confidence interval for another
parameter can be used to make sure Wald intervals are contained in the
parameter space.

The new parameter
$$
   \theta = \mathop{\rm logit}(\pi) = \log\left(\frac{\pi}{1 - \pi}\right)
$$
is called the *canonical parameter* of the binomial distribution.
This is a term from the
[theory of exponential families of distributions](expfam.html).
This function is so important that it has a special name,
pronounced with a soft g: low-jit.
It is also the default link function for `family = binomial` for R function
`glm`, which is why it was used in
the [preceding example](#likelihood-intervals).

This function is invertible, and its inverse function is,
as was seen in the preceding example,
$$
   \pi = \mathop{\rm logit}\nolimits^{-1}(\theta)
   = \frac{\exp(\theta)}{1 + \exp(\theta)}
   = \frac{1}{1 + \exp(- \theta)}
$$

As $\pi$ goes from zero to one, $\theta$ goes from $- \infty$ to $+ \infty$.
So the interval for $\theta$ can never be outside the parameter space
for $\theta$, and must map into the parameter space for $\pi$.

So we make a confidence interval for $\theta$, which R function `glm`
knows how to do, and then transform that to being a confidence interval
for $\pi$.
```{r wald-fixup}
crit <- qnorm((1 + conf.level) / 2)
pout <- predict(gout, se.fit = TRUE)
(pout$fit + c(-1, 1) * crit * pout$se.fit) |> invlogit()
```

This trick of putting the endpoints of an interval through a function
to give another interval only works when the function is monotone.
But every scalar-to-scalar invertible function is monotone.
So we do not need to check that.

#### Fixing Zero-Width Wald Intervals

As mentioned above, Wald intervals have zero width when the data are on
the boundary of the sample space.
If $x = 0$, then $\hat{\pi}_n = 0$.
If $x = n$, then $\hat{\pi}_n = 1$.
In either case, $\hat{\pi}_n (1 - \hat{\pi}_n) = 0$
so the Wald interval has zero width.

Geyer (2009, [DOI: 10.1214/08-EJS349](https://doi.org/10.1214/08-EJS349))
proposed a fixup for this.  It [works for all of the statistical models
discussed in this course](infinity.html#theory-confidence).  But for
binomial data it simplifies to the following.

 * If $x = 0$, then a $1 - \alpha$ confidence interval is
   $0 < \pi < 1 - \alpha^{1 / n}$.

 * If $x = n$, then a $1 - \alpha$ confidence interval is
   $\alpha^{1 / n} < \pi < 1$.

## Summary of Frequentist Inference for Proportions

```{r wald.test, echo=FALSE}
z <- (pi.hat - pi) / sqrt(pi.hat * (1 - pi.hat) / n)
save.wald.p.value <- 2 * pnorm(- abs(z))
save.wald.p.value <- formatC(save.wald.p.value, format = "f", digits = 3)

save.wald.interval.theta <-
    (pout$fit + c(-1, 1) * crit * pout$se.fit) |> invlogit()
save.wald.interval.theta <-
    formatC(save.wald.interval.theta, format = "f", digits = 3)
```

There is no one right way to do frequentist statistics.  Just for
two-tailed tests and two-sided intervals we had

 * the usual (or Rao or score) test, $P = `r save.rao.p.value`$,

 * the likelihood ratio (or Wilks) test, $P = `r save.wilks.p.value`$,

 * the Wald test, $P = `r save.wald.p.value`$ (calculations not shown),

 * the usual (or Wald) confidence interval,
   (`r save.wald.interval[1]`, `r save.wald.interval[2]`),

 * the Wald confidence interval for $\theta$ mapped to an interval for $\pi$,
   (`r save.wald.interval.theta[1]`, `r save.wald.interval.theta[2]`),

 * the score (or Rao) confidence interval,
   (`r save.rao.interval[1]`, `r save.rao.interval[2]`), and

 * the likelihood (or Wilks) confidence interval,
   (`r save.wilks.interval[1]`, `r save.wilks.interval[2]`).

But the [dogma of only one procedure](#multiple-tests-or-intervals)
says you only do one of these, and you choose which one to do before
seeing the data (unless there is a correction for doing more than one).

## Bayesian Inference

Thomas Bayes ([Wikipedia article](https://en.wikipedia.org/wiki/Thomas_Bayes))
died in 1761 by which time he had written an unpublished note about the
binomial distribution and what would now be called Bayesian inference for
it using a flat prior.  The note was found by a friend and read to the Royal
Society of London in 1763 and published in its *Philosophical Transactions*
in 1764 thus becoming widely known.

### Prior

We have learned a little about Bayesian inference since then, and
we now would rather allow any prior distribution that is a beta
distribution
([Wikipedia article](https://en.wikipedia.org/wiki/Beta_distribution)).

This is mathematically convenient because it then follows that the posterior
distribution is also a beta distribution.  A family that has this property
(if the prior is in the family, then so is the posterior) is called
the *conjugate family of distributions* to the distribution of the data.
So here the beta family of distributions (for the parameter)
is conjugate to the binomial family of distributions (for the data).

The beta distribution is a
distribution on the interval $(0, 1)$ that has two
strictly positive parameters $\alpha_1$ and $\alpha_2$ that control
the mean and variance
of the distribution, but not in any obvious way.  Moreover,
if $f_{\alpha_1, \alpha_2}$ denotes the probability density function, then
$$
   \lim_{x \to 0}
   f_{\alpha_1, \alpha_2}(x) = \begin{cases}
   \infty, & 0 < \alpha_1 < 1
   \\
   \alpha_2, & \alpha_1 = 1
   \\
   0, & 1 < \alpha_1 < \infty
   \end{cases} 
$$
and, similarly,
$$
   \lim_{x \to 1}
   f_{\alpha_1, \alpha_2}(x) = \begin{cases}
   \infty, & 0 < \alpha_2 < 1
   \\
   \alpha_1, & \alpha_2 = 1
   \\
   0, & 1 < \alpha_2 < \infty
   \end{cases} 
$$
so the behavior is qualitatively different near the edges of the
sample space for different values of $\alpha_1$ and $\alpha_2$.
Nevertheless, $f_{\alpha_1, \alpha_2}(x)$ is continuous in $\alpha_1$,
$\alpha_2$, and $x$ for $0 < \alpha_1 < \infty$ and $0 < \alpha_2 < \infty$
and $0 < x < 1$.  It is only the limits above that are discontinuous.

For $\alpha_1 = \alpha_2 = 1$, the beta distribution is uniform (flat)
on the interval $(0, 1)$.  Thus our allowing any beta distribution
includes the original calculation made by Bayes.

Recall that in Bayesian inference, it is the parameters that are
considered random (and the data not random after they are seen).
Thus these beta distributions are (in this Bayesian application)
considered distributions of the parameter $\pi$.  That is why
we want distributions whose sample space is the parameter space
for the binomial distribution (the interval $0 < \pi < 1$).

It is confusing when the prior and posterior distribution,
which are distributions
of the parameters of the given problem (here $\pi$), themselves have
parameters (here $\alpha_1$, and $\alpha_2$).  To avoid confusion, the
latter are called *hyperparameters*.  So

  * $\pi$ is the parameter and

  * $\alpha_1$ and $\alpha_2$ are the hyperparameters.

Bayesians treat parameters and hyperparameters very differently.

  * Parameters are considered random.  Their distributions (prior
    and posterior) describe our uncertainty about their true values.

  * Hyperparameters are not considered random.  They just say which
    prior and posterior we are talking about.

### Posterior

When applied to this problem Bayes' rule says,

 * if the prior is the beta distribution with hyperparameters
   $\alpha_1$ and $\alpha_2$,

 * then the posterior is the beta distribution
   with hyperparameters $x + \alpha_1$ and $n - x + \alpha_2$,

where (as above) the binomial data are $x$ successes in $n$ trials.

It takes a large amount of theory to derive this, which is beyond
the scope of this course.  It is done in theory courses
(for example [these notes,
slides 6—12](https://www.stat.umn.edu/geyer/5102/slides/s4.pdf#page=6)).

Now our Bayesian inference is a whole distribution (the posterior).
We can show it by plotting its probability density function.
We start with the $\alpha_1 = \alpha_2 = 1$ case that Bayes used
```{r binom-bayes-caption,echo=FALSE}
# need to do this way, because cannot substitute values of R objects
# into character string inside fig.cap in next code chunk.
# don't really understand but does not work.
foo <- paste0("Binomial data with x = ", x,
    " and n = ", n, ".  Solid curve posterior, dashed curve prior.")
```
```{r binom-bayes,fig.align="center",fig.cap=foo}
alpha1 <- 1
alpha2 <- 1
curve(dbeta(pi, x + alpha1, n - x + alpha2), from = 0, to = 1,
    xname = "pi", xlab = expression(pi), ylab = "probability density")
curve(dbeta(pi, alpha1, alpha2), from = 0, to = 1,
    xname = "pi", add = TRUE, lty = "dashed")
```

### Concentration

Hopefully, you can just see from the picture that the posterior is
more concentrated than the prior (the prior is more spread out).

If not, the variance of the posterior distribution is
$$
   \frac{(x + \alpha_1) (n - x + \alpha_2)}
   {(n + \alpha_1 + \alpha_2)^2 (n + 1 + \alpha_1 + \alpha_2)}
   \approx
   \frac{\hat{\pi}_n (1 - \hat{\pi}_n)}{n}
$$
for large $n$.  Thus we see that when $n$ is large

 * the prior does not matter and

 * the Bayesian (approximate) posterior variance is the same
   as the frequentist (approximate) sampling variance.

These two facts are a special case of a general phenomenon called
the Bernstein—von Mises theorem ([Wikipedia
page](https://en.wikipedia.org/wiki/Bernstein%E2%80%93von_Mises_theorem)).

We also see that the Bayesian agrees with the frequentist about the
square root law: statistical precision varies as the square root of
the sample size.

### Bayesian Learning Paradigm

Prior and posterior distributions are the same sort of thing.
In fact a prior can be a posterior.

A distribution can be the posterior distribution incorporating some
data we have seen but also the prior distribution for analyzing some data
we have not yet seen.

The posterior for the previous analysis serves as the prior for our
next analysis.  Whatever you call it, it describes our current
uncertainty about the true unknown parameter value.
Suppose we observe some more data.
```{r binom-bayes-newdata}
# update hyperparameters incorporating old data
alpha1 <- x + alpha1
alpha2 <- n - x + alpha2
# new data
x <- 3
n <- 7
```
```{r binom-bayes-caption-too,echo=FALSE}
foo <- paste0("Posterior Distribution. Binomial data with x = ", x,
    " and n = ", n, ".  Solid curve posterior, dashed curve prior.")
```
```{r binom-bayes-too,fig.align="center",fig.cap=foo}
curve(dbeta(pi, x + alpha1, n - x + alpha2), from = 0, to = 1,
    xname = "pi", xlab = expression(pi), ylab = "probability density")
curve(dbeta(pi, alpha1, alpha2), from = 0, to = 1,
    xname = "pi", add = TRUE, lty = "dashed")
```

Note that we can write the hyperparameters more symmetrically as

 * $x + \alpha_1 = \text{number of successes} + \alpha_1$ and

 * $n - x + \alpha_2 = \text{number of failures} + \alpha_2$.

Suppose we observe a sequence of data obtained in batches with
$x_i$ successes and $y_i$ failures in batch $i$.  Then

 * for the first batch.
     - prior has hyperparameters $\alpha_1$ and $\alpha_2$ and
     - posterior has hyperparameters $x_1 + \alpha_1$ and $y_1 + \alpha_2$
 * for the second batch.
     - prior has hyperparameters $x_1 + \alpha_1$ and $y_1 + \alpha_2$
       (same as previous posterior)
     - posterior has hyperparameters $x_2 + x_1 + \alpha_1$
       and $y_2 + y_1 + \alpha_2$
 * for the third batch.
     - prior has hyperparameters $x_2 + x_1 + \alpha_1$
       and $y_2 + y_1 + \alpha_2$
       (same as previous posterior)
     - posterior has hyperparameters $x_3 + x_2 + x_1 + \alpha_1$
       and $y_3 + y_2 + y_1 + \alpha_2$
 * and so forth.

At each stage the posterior for old data is the prior for new data.
And ― *this is very important* ― we get the same results for the same
data no matter how we subdivide the data into batches.

All Bayesian inference works this way, not just this particular example.
If Bayesian inference did not have this property, it would not be so
philosophically satisfying.  Since it does, it can serve as a model
for all learning.  Data arrives.  Uncertainty is reduced.
And Bayes' rule says exactly how.

### Point Estimates

Bayesians have posterior distributions, but sometimes one does not
want a whole picture.  Also, when there is more than one parameter,
pictures get harder to draw.

So one can just compute some numbers, for example
```{r bayes-point}
# current posterior hyperparameters
alpha1 <- x + alpha1
alpha2 <- n - x + alpha2
# posterior mean
alpha1 / (alpha1 + alpha2)
# posterior median
qbeta(0.5, alpha1, alpha2)
# posterior mode
((alpha1 - 1) / (alpha1 + alpha2 - 2)) |> max(0) |> min(1)
# posterior standard deviation
sqrt(alpha1 * alpha2 / (alpha1 + alpha2)^2 / (alpha1 + alpha2 + 1))
# posterior quartiles
qbeta(c(1, 3) / 4, alpha1, alpha2)
# posterior interquartile range
qbeta(c(1, 3) / 4, alpha1, alpha2) |> diff()
```
And this could keep on going forever.  Hence the modern emphasis
on posterior distributions rather than point estimates.

# Inference About a Two-Dimensional Contingency Table {#two-d-table}

## Contingency Tables

A *contingency table* is just a table of counts.  They come in any dimension,
1, 2, 3, 4, $\ldots$.
One-dimensional is covered in
[these notes](multi-simple.html#one-way-contingency-table).
Three-and-higher-dimensional is covered in
[these notes](ch4.html) and in
[these notes](ch9.html).

## Data

Here are some data [http://www.stat.umn.edu/geyer/5421/mydata/multi-simple-3.txt](http://www.stat.umn.edu/geyer/5421/mydata/multi-simple-3.txt).

We read them into R as follows.
The data set read in by the R function `read.table` below
(if this does not work on your computer [see this announcement on the
course home page](http://www.stat.umn.edu/geyer/5421/#download)).
```{r two-data-frame}
foo <- read.table(
    url("http://www.stat.umn.edu/geyer/5421/mydata/multi-simple-3.txt"),
    header = TRUE)
class(foo)
foo
```

These are data in the form that R function `glm` likes, a data frame.
If you have heard of the [tidyverse](https://www.tidyverse.org/), it
is based on the fundamental principle that all data should be put in this
form (otherwise it is not "tidy").

But humans and some other R functions and the textbook (Agresti) like
contingency tables.  We turn it into that as follows.
```{r two-data-table}
bar  <- xtabs(y ~ color + opinion, data = foo)
class(bar)
bar
```

In data frame form, we have three variables

 * the count variable `y`,

 * the categorical variable `color`, and

 * the categorical variable `opinion`.

(in a data frame, every column is a variable).

In the table form

 * the data are now in a two-dimensional array, and the count
   variable has no name, but corresponds to `y` in the data frame.

 * the categorical variables in the data frame are also gone, they
   correspond to the row and column labels in the table.

```{r two-data-table-labels}
dimnames(bar)
```

 * variable names in the data frame correspond to dimension names in the table
   and

 * variable values in the data frame correspond to row and column labels
   in the table.

Since we have a contingency table with
$`r nrow(bar)` \times `r ncol(bar)` = `r prod(dim(bar))`$ cells,
we also have that many possible parameters.

Since no one parameter is inherently more interesting than any other,
we will not be interested in confidence intervals, only hypothesis tests.

Since there are a lot of parameters, it is not clear what hypotheses to test.
The theory is complicated.  The actual procedures are a lot simpler, so
we will do practice before theory.  This means we will not know how to
interpret the tests until we have done the theory.  You cannot interpret
a hypothesis test until you are very clear about what the null hypothesis is.

## R Function `chisq.test`

```{r two-doit}
baz <- chisq.test(bar)
baz
```
Whatever the null hypothesis is (we don't know that yet),
we can reject it at level = 0.05.

R function `chisq.test` likes the table form of the data.

To go further, we need the table of estimated cell mean values,
also called expected values, or just *expected* for short.
```{r two-expected}
moo <- baz$expected
moo
```
These are the best estimates of the expected values *under the null hypothesis*
(and we don't even know what that is yet).

For short, we call the observed data *observed*.  The form of the
Pearson chi-squared test statistic is
\begin{equation} \tag{${\ddagger}$}
   \sum_\text{all cells}
   \frac{(\text{observed} - \text{expected})^2}{\text{expected}}
\end{equation}
We didn't actually need to know this because R function `chisq.test` knows
that.
```{r two-pearson-statistic}
tstat <- sum((bar - moo)^2 / moo)
all.equal(tstat, baz$statistic, check.attributes = FALSE)
```

The Pearson chi-squared test is a special case of the Rao (or score) test.
There is also a likelihood ratio (or Wilks) test.
It has a test statistic of this form.
\begin{equation} \tag{${\ddagger}{\ddagger}$}
   \sum_\text{all cells}
   2 \cdot \text{observed} \cdot
   \log \left( \frac{\text{observed}}{\text{expected}} \right)
\end{equation}
As always, these two tests are *asymptotically equivalent* (give nearly
the same results for sufficiently large $n$, and here the sample size
is large).
```{r two-n}
sum(bar)
```
So we use the same null distribution (distribution of the test statistic
under the null hypothesis), which is the chi-squared distribution with
```{r two-df}
baz$parameter
```
degrees of freedom (not completely sure why it calls degrees of freedom
"parameter", apparently so it can use the same print function for both
tests about one parameter and tests of model comparison).

So the LRT (likelihood ratio test) $P$-value is
```{r two-lrt}
tstat <- sum(2 * bar * log(bar / moo))
pchisq(tstat, df = baz$parameter, lower.tail = FALSE)
```

```{r two-both, echo=FALSE}
save.pearson.p.value <- baz$p.value
save.pearson.p.value <- formatC(save.pearson.p.value, format = "f", digits = 3)
save.lrt.p.value <- pchisq(tstat, df = baz$parameter, lower.tail = FALSE)
save.lrt.p.value <- formatC(save.lrt.p.value, format = "f", digits = 3)
```
So not much difference, $P = `r save.pearson.p.value`$ for the Pearson
chi-squared test (special case of Rao test) and
$P = `r save.lrt.p.value`$ for the likelihood ratio test.

**Warning:** the likelihood ratio test statistic calculation above blows up,
giving `0 * -Inf = NaN` when any of the observed values are zero, but
this is incorrect, just bad calculation.
Since
$$
   x \log(x) \to 0, \qquad \text{as $x \to 0$}
$$
the term in the sum for each cell having observed value zero should be zero,
so to be careful we calculate
```{r two-lrt-careful}
tstat <- 2 * bar * log(bar / moo) # no sum yet
tstat[is.nan(tstat)] <- 0
tstat <- sum(tstat) # now OK to sum
```

## R Function `glm` with `family = poisson`

So that was that.  We are now done with how to do two tests (Pearson and LRT)
for these data.  But that is not how we are going to do higher dimensional
tables.  For that we are going to use R function `glm`.  So we show how to
do that too.
```{r two-glm}
gout0 <- glm(y ~ color + opinion, family = poisson, data = foo)
gout1 <- glm(y ~ color * opinion, family = poisson, data = foo)
sout.lrt <- anova(gout0, gout1, test = "LRT")
sout.rao <- anova(gout0, gout1, test = "Rao")
```
Note that R function `glm` wants the data in data frame form.
```{r two-glm-too}
sout.lrt
sout.rao
all.equal(sout.rao[2, "Pr(>Chi)"], baz$p.value)
all.equal(sout.lrt[2, "Pr(>Chi)"],
    pchisq(tstat, df = baz$parameter, lower.tail = FALSE))
all.equal(sout.rao[2, "Pr(>Chi)"], baz$p.value,
    tolerance = 1e-5)
```
One or the other of these R functions is a bit sloppy about calculating
the Rao test.  No idea which.  Some sloppiness is inherent in the fact
that computer arithmetic is not exact (no more than 16 decimal place
precision, but errors can add up in complicated calculations; also
convergence can be sloppy in optimization algorithms ― R function
`chisq.test` does not do optimization but R function `glm` does).

## Theory

### Sampling Distributions for Hypothesis Tests {#sampling-dist}

So now we get to theory.  What are these hypothesis tests testing?
And that depends on

 * the family of distributions assumed for the data and

 * the null hypothesis (which is in the assumed family).

And here is where it gets complicated *even in intro stats*.
There are three different families of distributions that can be used

 * Poisson,

 * multinomial, and

 * product multinomial.

But intro stats only distinguishes two

 * tests of independence (of the categorical variables in the data frame,
   which are the row and column labels in the table) and

 * tests of homogeneity of proportions (that the rows (or columns but not
   both) all have the same proportions in the same positions).

The correspondence between intro stats and sophisticated courses like this
is

 * tests of independence correspond to Poisson *and* multinomial and

 * tests of homogeneity of proportions correspond to product multinomial.

So, first, a little bit about sampling models (families of distributions).
More on that [here](ch1.html#the-relationship-of-the-various-sampling-schemes)
and [here](sampling.html).

 * Poisson sampling is for when you go out and collect data for a fixed
   amount of time and the time at which you collect each data point and
   which cell of
   the table it goes in are unrelated to any other data point.

 * multinomial sampling is like Poisson except the total sample size
   (sum of all counts) is fixed in advance and you continue sampling
   until you get that many and then stop.

 * product multinomial sampling is like multinomial except the row sums
   or column sums but not both of the table are fixed in advance and
   you continue sampling until you get that many in each row (respectively,
   column) and then stop (in that row (respectively, column); you keep
   sampling the others until they get the required number).

So the difference is what is random and what is fixed.

   * In Poisson, all counts are random.

   * In multinomial, all counts are random subject to the sum being fixed.

   * In product multinomial, all counts are random subject to the row
     (respectively column) sums being fixed.

And now comes the theoretical shocker.  If we could do it (but we cannot)
*exact* inference (based on *exact* sampling distributions of test statistics)
would differ for the three sampling schemes. But the asymptotic approximations
are *the same!*  That is, whether we use Poisson, multinomial, or product
multinomial for the sampling scheme (family of distributions) and whether
we use the Wilks (also called likelihood ratio), Rao (also called score)
or Wald test, we have (for large $n$)

 * nearly the same test statistic,

 * the same large sample approximate null distribution of the test
   statistic (chi-squared with some degrees of freedom, which is the same
   for all), and hence

 * nearly the same (approximate, large $n$) $P$-value.

That is why all of our computing examples above *ignored the sampling scheme*.
You get the same answer (using large sample approximation) for all three
sampling schemes.

As with tests about binomial proportions, we did not illustrate Wald tests.
This is because textbooks do not recommend them, even though they are
asymptotically equivalent (work no better and *no worse* for large $n$)
than the other two.

Also R function `anova` does not offer the Wald test as an option.

### The Null Hypothesis.

#### All Sampling Schemes

We are going to discuss the null hypothesis for *data in table form*.
The notation would be different ([discussed later in the course](ch4.html))
for data frame form.
But, of course, the data are the data, no matter how printed or stored
in the computer.  So data frame versus table make no difference to the
statistical analysis.

Mathematically, the table is a matrix.  The observed values are $y_{i j}$
where $i$ is the row index and $j$ is the column index, that is,
$y_{i j}$ is the count in the cell in row $i$ and column $j$.
The corresponding
parameter is $\mu_{i j}$.  These parameters are the mean values (or
expected values) for the cells.  The families of distributions are

 * Poisson: the data in the cells of the table are independent Poisson
   random variables.

 * Multinomial: the data in the cells of the table are components of
   a multinomial random vector.

 * Product multinomial: the rows (or columns) of the table are independent
   multinomial random vectors.

More on what [independent means in probability and statistics](ch1.html#iid)
here.  You should have some idea about that from previous statistics courses.

All three sampling distributions are determined by (the collection of)
means $\mu_{i j}$.  Thus those are parameters for all the models.

The null hypothesis uses the multiplicative assumption
$$
   \mu_{i j} = \alpha_i \beta_j, \qquad \text{for all $i$ and $j$},
$$
where the $\alpha_i$ and $\beta_j$ are the parameters of the family of
distributions that is the null hypothesis.

The reason the tests discussed above are called tests of independence is
[whenever means are products, that comes from independence](ch1.html#iid).

There is something a bit funny about this parameterization in that it
is not *identifiable*.  (You may have seen this topic discussed in other
statistics courses using the term *collinearity* or
[multicollinearity](https://en.wikipedia.org/wiki/Multicollinearity).)
If we divide all of the $\alpha_i$ by a constant
$c$ and multiply all the $\beta_j$ by the same constant we get the same
means
\begin{equation} \tag{$\S$}
   \mu_{i j} = \alpha_i \beta_j = \frac{\alpha_i}{c} \cdot c \beta_j
\end{equation}
So the $\alpha_i$ and $\beta_j$ uniquely determine the $\mu_{i j}$
but *not vice versa*.
Thus we cannot estimate all of the alphas and all of the betas:
they are not *identifiable*.  But we can fix any one of the alphas or
betas at an arbitrary positive value, and then estimate the rest
(this just amounts to choosing a particular constant $c$).

#### Multinomial Sampling

When we are assuming multinomial sampling, we often use a different
parameterization
$$
   \mu_{i j} = n \pi_{i j}
$$
so
$$
   \pi_{i j} = \frac{\mu_{i j}}{n}
$$
where (as usual) $n$ is the sample size, the total number of counts in the
whole table, the total number of individuals surveyed.

Recall that in multinomial sampling $n$ is not a parameter.  It is a number
fixed in advance of any data being collected.  It is a known constant,
unlike the parameters which are *unknown* constants.

Also probabilities sum to one, so we must have
\begin{align}
   \textstyle \sum_{i j} \pi_{i j} & = 1
   \\
   \textstyle \sum_{i j} \mu_{i j} & = n
\end{align}

Of course since
the $\pi_{i j}$ determine the $\mu_{i j}$ and vice versa, they are equivalent;
they parameterize the same model.
And if the $\mu_{i j}$ have the multiplicative form,
then so do the $\pi_{i j}$ and vice versa.  But for these
we sometimes redefine the alphas and betas, writing
$$
   \pi_{i j} = \alpha_i \beta_j
$$
so
$$
   \mu_{i j} = n \alpha_i \beta_j
$$

Now we have the restriction that probabilities must sum to one, so
$$
   \textstyle \sum_{i j} \pi_{i j} = 1 =
   \left(\sum_i \alpha_i\right)
   \left(\sum_j \beta_j\right)
$$
and we see that one way to obtain this is to have the $\alpha_i$ to sum
to one and the $\beta_j$ to sum to one.  We don't have to do this.
We only need the $\pi_{i j}$ to sum to one.
But this is convenient and makes $\alpha$ and $\beta$ probability vectors
(have components that are nonnegative and sum to one).

Introduce some convenient notation
\begin{align*}
  y_{i +} & = \textstyle \sum_j y_{i j}
  \\
  y_{+ j} & = \textstyle \sum_i y_{i j}
  \\
  y_{+ +} & = \textstyle \sum_{i j} y_{i j}
\end{align*}
then the
[best estimates of the parameters of the null hypothesis](likelihood.html) 
are
\begin{align}
  \hat{\alpha}_i & = \frac{y_{i +}}{y_{+ +}}
  \\
  \hat{\beta}_j & = \frac{y_{+ j}}{y_{+ +}}
  \\
  \hat{\pi}_{i j} & = \hat{\alpha}_i \hat{\beta}_j
  \\
  \hat{\mu}_{i j} & = n \hat{\alpha}_i \hat{\beta}_j
\end{align}
We have hats on everything to remind us that these are just parameter estimates
(not the true unknown parameter values).

#### Poisson Sampling

The same numbers are the parameter estimates for Poisson sampling.
If we want to not mention the $\pi_{i j}$ because they only make
sense for multinomial sampling, then we can plug earlier equations into
the last equation to get
\begin{equation} \tag{${\S}{\S}$}
  \hat{\mu}_{i j} = \frac{y_{i +} y_{+ j}}{y_{+ +}}
\end{equation}

In either case our formulas agree with the output of R functions `chisq.test`
and `glm`.
```{r two-check-expected}
mu.hat <- outer(rowSums(bar), colSums(bar)) / sum(bar)
all.equal(mu.hat, baz$expected, check.attributes = FALSE)
qux <- transform(foo, mu.hat = gout0$fitted.values)
moo <- xtabs(mu.hat ~ color + opinion, data = qux)
moo <- unclass(moo)
all.equal(mu.hat, moo, check.attributes = FALSE)
```

Sometimes it can be a real pain to get R to check that two different
calculations give the same results (except for inaccuracy of computer
arithmetic).  Here the means produced by R function `glm` are a vector.
So we had to use R function `xtabs` to put the means in the right cells
of the contingency table.  And we had to use R function `unclass` to
get rid of the class attribute for R object `moo` because otherwise
R function `all.equal` says one is a `"matrix"` and the other `"xtabs"`
and won't do the comparison.

#### Product Multinomial Sampling

Now we still have multiplicative form of the null hypothesis and
(${\S}{\S}$) still gives the best estimates of the parameters under the
null hypothesis.  But now we write the null hypothesis
\begin{equation} \tag{${\S}{\S}{\S}$}
   \mu_{i j} = n_i \beta_j
\end{equation}
if we are assuming that row sums of the data ($n_i$) were fixed in
advance and
$$
   \mu_{i j} = n_j \alpha_i
$$
if we are assuming that column sums of the data ($n_j$) were fixed in
advance.  Now our only parameters are $\beta_j$ in the first and
$\alpha_i$ in the second.

And we say the null hypothesis specifies *homogeneity of proportions*,
that is, the same vector $\beta$ with components $\beta_j$ applies to each
row of the table in the first case,
and, the same vector $\alpha$ with components $\alpha_i$ applies to each
column of the table in the first case,

All of this seems complicated, but for all the sampling schemes
the null hypothesis has the form
$$
   \text{expected}_{i j} = (\text{something})_i (\text{something else})_j
$$
regardless of what notation we use, and regardless of what parts are
unknown parameters and what parts are known sample sizes.

Also for all three sampling schemes

 * the likelihood ratio test statistic is the same,

 * the (approximate, large $n$) null distribution of the test statistic
   is the same, and hence

 * the (approximate, large $n$) $P$-value is the same

Of course this is only for the different sampling schemes.
For the different test statistics (Wilks, Rao, Wald),
repeating what was said above, for large $n$,

 * the test statistics are nearly the same,

 * the approximate null distribution of the test statistics
   are exactly the same, and hence

 * the approximate $P$-values are nearly the same.

### The Alternative Hypothesis

For all of these procedures, the alternative hypothesis is "anything goes",
that is the parameters can be anything other than in the null hypothesis.

### Degrees of Freedom

#### General Rule

The general rule for degrees of freedom for Wilks, Rao, and Wald tests is:
degrees of freedom is the number of identifiable parameters in
the alternative hypothesis minus the number of identifiable
parameters in null hypothesis.

Here *identifiable* means different identifiable parameter vectors
correspond to different distributions.  There
are no constraints on their values.  In multinomial and product multinomial
we do have some constraints.  And those constraints change the numbers
of identifiable parameters.

#### Poisson

Suppose our table has $r$ rows and $c$ columns.  Then the $\mu_{i j}$
can be any nonnegative numbers.  So

 * the alternative hypothesis has $r c$ identifiable parameters, and

 * the null hypothesis has $r + c - 1$ identifiable parameters.

The latter is because
the adjustable constant $c$ in ($\S$) can be chosen arbitrarily,
and this makes one of the alphas (or one of the betas but not both)
a known constant (rather than an unknown parameter to be estimated).
For example, if we choose $c = \alpha_1$, then the new $\alpha$ vector
(the old $\alpha$ vector divided by $c$) has first component $\alpha_1 = 1$.

Hence the degrees of freedom is
$$
   (r c) - (r + c - 1) = (r - 1) (c - 1)
$$

#### Multinomial

Now

 * the alternative hypothesis has $r c - 1$ identifiable parameters, and

 * the null hypothesis has $r + c - 2$ identifiable parameters.

The reason for the first is the constraint that the $\pi_{i j}$ must
sum to one, so any single $\pi_{i j}$ is determined by the rest.

The reason for the second is the constraints that the $\alpha_i$ must
sum to one, and so must the $\beta_j$.  Hence there are $r - 1$
identifiable alphas and $c - 1$ identifiable betas.

Hence the degrees of freedom is
$$
   (r c - 1) - (r + c - 2) = (r - 1) (c - 1)
$$

#### Product Multinomial

We do the case where row sums are fixed in advance so the null hypothesis is
(${\S}{\S}{\S}$).  The other case is similar.  Now

 * the alternative hypothesis has $r (c - 1)$ identifiable parameters, and

 * the null hypothesis has $c - 1$ identifiable parameters.

The latter is because the $\beta_j$ must
sum to one, so any single $\beta_j$ is determined by the rest.

The alternative hypothesis has different $\beta$ vectors for different
rows, and each has $c - 1$ identifiable parameters.  So $r (c - 1)$.

Hence the degrees of freedom is
$$
   r (c - 1) - (c - 1) = (r - 1) (c - 1)
$$

We waded through all of this algebra to give a concrete example.
But we needn't have.
[It is a theorem](https://www.stat.umn.edu/geyer/8112/notes/tests.pdf)
that Wilks, Rao, and Wald tests comparing the same hypotheses always
have the same asymptotic distribution, which is chi-squared with the
same degrees of freedom.
[And it is another theorem](sampling.html#thm:df) that the degrees
of freedom do not depend on the sampling scheme; it is the same
for Poisson, multinomial, or product multinomial sampling.

So once we have done one (correct) degrees of freedom calculation,
we know all the others have to agree *by theory*.

Also, R functions can usually figure out the degrees of freedom for us.
We only need to figure it out by theory when there is no available R function
to do the analysis.

# Interpretation

A standard thing to ask students on homeworks or tests is to "interpret"
your results, that is, not just give some numbers or pictures but rather
to explain what they mean.

## Bayesian

Bayesian inference almost needs no interpretation.  A posterior distribution
or a prior distribution is a probability distribution, and if someone doesn't
know what that is, then they don't have enough education to understand
Bayesian inference ([here](ch1.html#pmf-pdf) are our notes on that).

One point is that you must make clear what parameterization is being used.
We have [already seen](#outside) that there can be more than one parameter
used to describe a statistical model.  So you have to be clear about which
one.

## Frequentist Confidence Intervals

This is even more essential for confidence intervals.  A confidence interval
is an interval estimate of a *parameter*.  If you don't say which parameter,
then your interval is useless.  In our
[example of fixing up a Wald interval](#outside) we had an interval for $\pi$
```{r interval-for-pi}
(pout$fit + c(-1, 1) * crit * pout$se.fit) |> invlogit()
```
and this corresponds to an interval for $\theta$
```{r interval-for-theta}
(pout$fit + c(-1, 1) * crit * pout$se.fit)
```
Clearly, if you don't say which interval is for which parameter, what you
are saying is meaningless.

Other points to make about interpreting confidence intervals are just
facts about frequentist inference

 * the confidence interval may not cover the true unknown parameter value
   (it's called a 95% confidence interval because it misses 5% of the time),

 * unlike Bayesians, we are not treating the parameter as random, what
   is random is the endpoints of the confidence interval, and

 * you should not refer to the parameter as some Greek letter but rather
   something your clients or readers can understand: say "the proportion
   of Minnesotans whose favorite ice cream flavor is mint chocolate chip"
   (or whatever the parameter is).

## Frequentist Hypothesis Tests

A $P$-value interprets itself.  That is why the concept was invented.
A $P$-value is not made more understandable by declaring that it means
"statistically significant" results or attaching other adjectives, like
"borderline statistically significant" or some such.

People *want* clear cut answers.  Your theory is either proved or not.
But statistics *does not give clear cut answers*.
Trying to make it seem like it does is always wrong.

$P$-values are numbers between zero and one.  Low $P$-values (near zero)
are evidence against the null hypothesis, the lower the $P$-value, the
stronger the inference.  High $P$-values (near one) are still evidence
against the null hypothesis, but the higher the $P$-value, the weaker
the inference.  For $P$-values above $0.1$ the evidence against the null
hypothesis is so weak as to be unimpressive.  Even more so for higher
$P$-values.

Some intro stats books outlaw the phrase "accept the null hypothesis"
(appropriate when the $P$-value is large) and insist on the phrase
"fail to reject the null hypothesis".

The point is, that the $P$-value is only relevant to the data used in the
analysis.  A large $P$-value only says that *these data* do not show any
"statistically significant" lack of fit of the null hypothesis.  More data
(if you were to get some) might show lack of fit (or might not).

The traditional dividing line between $P$-values being "low" or "high" is
0.05.  But this tradition has nothing to recommend it.  It is just laziness
on the part of those using it.  True, 0.05 is a nice "round" number,
but it is only a round number if we use the decimal system, which we
consider natural because people have ten fingers.  So using 0.05 as a cutoff
for $P$-values is trying to decide important scientific (business, sports,
whatever) questions by counting on your fingers.

I once wrote

> Anyone who thinks there is an important difference between
> $P = 0.049$ and $P = 0.051$ understands neither science nor statistics.

But my coauthor cut it — not because it was wrong but because it might offend.

Don't you ever seem to run afoul of this dictum.  Never say anything that
can be interpreted as saying there is an important difference between
$P = 0.049$ and $P = 0.051$.

$P$-values fall on a continuum, and no point on that continuum is special.

As to what $P$-values are, the answer is tricky.

When we have a point null hypothesis (that completely specifies a probability
model), then the $P$-value is a probability.  It is the probability *assuming
the null hypothesis is true* of seeing a result that gives at least
as much *apparent* evidence *against* the null hypothesis.  In short, it
is the probability of a false positive, of *falsely* claiming to have
statistically significant results.

When we are doing approximate rather than exact inference ― using large
sample approximation, as we are always doing in this course (except
for [fuzzy](fuzzy.html)) ― then a $P$-value is the same except for
approximately.  It is *approximately* the probability of a false positive.

So the lower the $P$-value, the lower the probability of *falsely* claiming
to have discovered something.

But that probability of getting the wrong answer is never zero.

When we have a composite null hypothesis, for example,
\begin{align*}
   H_0 : \text{true unknown parameter value} \ge \theta
   \\
   H_1 : \text{true unknown parameter value} < \theta
\end{align*}
for a one-tailed test, a $P$-value is not a probability but rather an
upper bound on probabilities, but we won't worry about that.
This does not apply to approximate (large sample) tests.
For those a $P$-value is still an approximate probability.

One last thing about hypothesis tests.  If you don't mention what the
null hypothesis is mathematically and what (if anything) this means
in scientific (business, sports, whatever) terms, then you haven't explained
anything about the hypothesis test.  It is all about the null hypothesis.
So if you don't understand the null hypothesis, then you don't understand
anything about the hypothesis test.

For example, for the [null hypothesis of independence](#the-null-hypothesis)
discussed
above, one should say which random variables, *color* and *opinion*, have
been found to be dependent ($P = `r sout.lrt[2, "Pr(>Chi)"] |> round(4)`$).

One should also be very careful to avoid overclaiming what the test says.

 * We all know *correlation is not causation* and *stochastic dependence
   isn't either*.  So don't say that!

 * Some dependence does not mean that any particular theory (causal or not)
   about the particular form of the dependence is true.  Rejecting the null
   hypothesis means there is *some* dependence of *some sort or other*,
   not necessarily any particular form.
