---
title: "Stat 5421 Lecture Notes: Completion of Exponential Families"
author: "Charles J. Geyer"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  bookdown::html_document2:
    number_sections: true
    md_extensions: -tex_math_single_backslash
    mathjax: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML
  bookdown::pdf_document2:
    number_sections: true
    md_extensions: -tex_math_single_backslash
linkcolor: blue
urlcolor: blue
bibliography: expfam.bib
csl: journal-of-the-royal-statistical-society.csl
link-citations: true
---

# License

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

# R

```{r libraries}
# library(alabama)
library(rcdd)
```

 * The version of R used to make this document is `r getRversion()`.

 * The version of the `alabama` package used to make this document is
   `r packageVersion("alabama")`.

 * The version of the `numDeriv` package used to make this document is
   `r packageVersion("numDeriv")`.

 * The version of the `rcdd` package used to make this document is
   `r packageVersion("rcdd")`.

 * The version of the `rmarkdown` package used to make this document is
   `r packageVersion("rmarkdown")`.

# Introduction

The purpose of this handout is to present two toy problems from @agresti
which are two-dimensional exponential families (hence graphable) and illustrate
situations where the maximum likelihood estimate (MLE) for the canonical
parameter does not exist, although it does exist as a limit of distributions
in the family.  The MLE canonical parameter values can be thought of being "at
infinity."  The MLE mean value parameter values are on the relative boundary
of convex support of the family (the smallest closed convex set containing
the canonical statistic with probability one).

@agresti calls these models examples of "complete separation" and
"quasi-complete separation."  These terms do not generalize beyond logistic
regression (or classification with two classes).  The general terminology
that applies to all exponential families that describes these toy models
is

 * MLE distribution concentrated at a single point (instead of complete
   separation) and

 * MLE distribution not concentrated at a single point but
   concentrated on a proper subset of the support of the original model
   (instead of quasi-complete separation).

We use the methods of @geyer-gdor to analyze these toy models.

# Example I {#one}

## Data

Section 6.5.1 of @agresti introduces the notion of complete separation
with the following example.
```{r example.i}
x <- seq(10, 90, 10)
x <- x[x != 50]
x
y <- as.numeric(x > 50)
y
```

The following figure shows these data.
```{r label="one", echo=FALSE, fig.align = "center", fig.cap = "Logistic Regression Data for Example I"}
plot(x, y)
```

## Attempt to Fit Model

Suppose we want to do "simple" logistic regression (one predictor $x$ plus
intercept, so the model is two-dimensional).  Let's try to do it naively.
```{r naive, error = TRUE}
gout <- glm(y ~ x, family = binomial)
summary(gout)
```
R function `glm` does give a warning.  But what are you supposed
to do about it?  If the output of the R function `summary` is to
be taken seriously, you cannot tell whether either regression coefficient
is nonzero.  As we shall see, that is complete nonsense.
Both parameters have to go to infinity (one to plus infinity, the other
to minus infinity) in order to maximize the likelihood.

In fact, these data are not analyzable by the R function `glm` and
its associated functions (generic functions having methods for
class `"glm"`).  So we will use the theory of Barndorff-Nielsen
completions of exponential families from @geyer-gdor.

## Submodel Canonical Statistic

That theory tells us that we
must look at the set of all possible values of the canonical statistic
$M^T y$ where $M$ is the model matrix and $y$ is the response vector.
For the model, $M$ has two columns:
the first column is all ones (the "intercept" column) and the second column
is $x$.  So let's find that set.  There are $2^n$ possible values where $n$
is the dimension of the response vector because each component of $y$ can be
either zero or one.  The following code makes all of those vectors.
```{r example.i.all.y}
yy <- matrix(0:1, nrow = 2, ncol = length(x))
colnames(yy) <- paste0("y", x)
yy <- expand.grid(as.data.frame(yy))

head(yy)
dim(yy)
```

But there are not so many distinct values of the submodel canonical statistic.
```{r example.i.all.m.transpose.y}
m <- cbind(1, x)
mtyy <- t(m) %*% t(yy)
t1 <- mtyy[1, ]
t2 <- mtyy[2, ]
t1.obs <- sum(y)
t2.obs <- sum(x * y)
```
```{r label="two",echo=FALSE,fig.align="center",fig.cap="Possible values of the submodel canonical statistic vector $M^T y$ for the data shown in Figure \\@ref(fig:one).  Solid dot is the observed value of the submodel canonical statistic vector."}
t.key <- paste(t1, t2, sep = ":")
t.idx <- match(unique(t.key), t.key)
t1.uniq <- t1[t.idx]
t2.uniq <- t2[t.idx]

plot(t1.uniq, t2.uniq, xlab = expression(t[1] == sum(y)),
    ylab = expression(t[2] == sum(x * y)))
points(t1.obs, t2.obs, pch = 19)
```

Figure \@ref(fig:two) shows these possible values of the submodel canonical
statistic.

And now we are stuck.  Figure \@ref(fig:two) seems to show that the observed
data vector is an extreme value, but we cannot easily figure out the direction
of recession.

## Computing the Limiting Conditional Model

Following @geyer-gdor we compute the support of the limiting conditional
model (LCM) and a generic direction of recession (GDOR) as follows
```{r linearity.and.gdor}
tanv <- m
tanv[y == 1, ] <- (- tanv[y == 1, ])
vrep <- cbind(0, 0, tanv)
lout <- linearity(vrep, rep = "V")
p <- ncol(tanv)
hrep <- cbind(-vrep, -1)
hrep <- rbind(hrep, c(0, 1, rep(0, p), -1))
objv <- c(rep(0, p), 1)
pout <- lpcdd(hrep, objv, minimize = FALSE)
gdor <- pout$primal.solution[1:p]
```

This computes two vectors `lout` and `gdor`.  The first is an index vector
giving the components of the response vector that are random in the LCM
(not conditioned on being equal to their observed values).
```{r show.linearity}
lout
```
Here we have found that no components of the response vector are random
in the LCM.  The LCM is the completely degenerate distribution that says
the response vector could not have had any value other than its observed
value.
The second is a GDOR.
```{r show.gdor}
gdor
```
This shows the direction the parameters ("coefficients") go to infinity
to obtain the LCM.  The first coefficient (the intercept) goes to minus
infinity `r abs(gdor[1] / gdor[2])` times faster than the second coefficient
(the slope) goes to plus infinity.

This computation of the GDOR is not general.  It only works in this
complete separation case.

## Confidence Regions

### Theory {#theory}

@geyer-gdor proposes a method of making confidence regions when
the MLE for the canonical parameter does not exist (when the data are on
the relative boundary of the convex support).  Suppose $\delta$ is a GDOR,
$Y$ is the response vector, $y$ its observed value, $M$ is the model matrix,
and $\eta = M \delta$.  If we use $\langle Y, \eta \rangle$ as a test statistic
for an upper-tailed test, then we know
$\langle Y, \eta \rangle \le \langle y, \eta \rangle$ almost surely by
definition of direction of recession.
We denote by $H$ the hyperplane that is the support of the LCM, the event
$\langle Y, \eta \rangle = \langle y, \eta \rangle$.

The $P$-value for the hypothesis test is $\Pr(Y \in H)$ since its observed
value is its largest possible value.
We can make a $100 (1 - \alpha)\%$ confidence region by inverting the
level $\alpha$ test.
So the confidence region is the set of parameter
values that put probability at least $\alpha$ on the support of the LCM.

### For Submodel Canonical Parameters

Let $\theta$ denote the saturated model canonical parameter vector
and $\beta$ the submodel canonical parameter vector, which are related
by $\theta = M \beta$.

Define a function $f$ by
$$
   f(\beta) = \log \Pr\nolimits_\beta(Y \in H)
$$
We want to plot the curve which is the locus of points $\beta$ such
that $f(\beta) = \log \alpha$.  Of course, we have
\begin{equation} \label{eq:effing}
   f(\beta)
   =
   \sum_{i = 1}^n \bigl[ y_i \log p_i + (1 - y_i) \log (1 - p_i) \bigr]
\end{equation}
As in the [notes on Agresti Chapter 4](ch4.html) we are careful computing
these probabilities.

Hence we can write this probability
calculation as an R function as follows.
```{r prob.calc.beta}}
alpha <- 0.05
fred <- function(beta) {
    theta <- m %*% beta
    logp <- ifelse(theta < 0, theta - log1p(exp(theta)),
        - log1p(exp(- theta)))
    logq <- ifelse(theta < 0, - log1p(exp(theta)),
        - theta - log1p(exp(-theta)))
    sum(logp[y == 1]) + sum(logq[y == 0]) - log(alpha)
}
```

Now we find one point inside the confidence region.
We can do that by starting at any point and moving in the GDOR.
```{r find.point.in.curve}
beta <- c(0, 0)
while (fred(beta) < 0) beta <- beta + gdor
beta
fred(beta)
```

Then we find points on the curve in various directions by solving
a univariate equation.
```{r find.points.on.curve,cache=TRUE}
beta.start <- beta
theta <- seq(0, 2 * pi, length = 10001)
ss <- double(length(theta))
for (i in seq(along = theta)) {
    sally <- function(s) {
        beta <- beta.start + s * c(sin(theta[i]), cos(theta[i]))
        fred(beta)
    }
    uout <- try(uniroot(sally, c(0, 1), extendInt = "downX"), silent = TRUE)
    ss[i] <- if (inherits(uout, "try-error")) NaN else uout$root
}
beta1 <- beta.start[1] + ss * sin(theta)
beta2 <- beta.start[2] + ss * cos(theta)
```
```{r betaConfidenceRegion,echo=FALSE,fig.align="center",fig.cap="Confidence Region for Coefficients"}
idx <- abs(beta1) <= 5 & abs(beta2) <= 5
plot(beta1[idx], beta2[idx], xlab = expression(beta[1]),
    ylab = expression(beta[2]), type = "n")
polygon(beta1, beta2, density = 3)
```

### For Submodel Mean Value Parameters

Let us map the confidence region shown in Figure \@ref(fig:betaConfidenceRegion)
to the submodel mean value parameter scale.  We know from exponential
family theory that if $Y$ is the response vector, and $\mu = E(Y)$ is
the saturated model mean value parameter vector, and $M$ is the model
matrix, then the submodel mean value parameter vector is $\tau = M^T y$.
```{R beta.to.tau}
betas <- cbind(beta1, beta2)
thetas <- m %*% t(betas)
thetas <- t(thetas)
# as with betas, rows of thetas are parameter vectors
mus <- 1 / (1 + exp(- thetas))
taus <- t(m) %*% t(mus)
taus <- t(taus)
# as with betas, rows of taus are parameter vectors
# in this case submodel mean value parameter vectors
dim(taus)
```

And we plot this confidence region along with the support
of the submodel canonical statistic vector
(Figure \@ref(fig:tauConfidenceRegion)).
```{r tauConfidenceRegion,echo=FALSE,fig.align="center",fig.cap="Confidence Region for Submodel Mean Value Parameter (gray).  Dots as in Figure \\@ref(fig:two)."}
# replace NaN values by MLE value = observed value
idx <- which(apply(taus, 1, anyNA))
taus <- rbind(taus[seq(1, min(idx) - 1), ], c(t1.obs, t2.obs),
    taus[seq(max(idx) + 1, nrow(taus)), ])
plot(t1.uniq, t2.uniq, type = "n",
    xlab = expression(tau[1]), ylab = expression(tau[2]))
polygon(taus[ , 1], taus[ , 2], border = NA, col = "gray")
points(t1.uniq, t2.uniq)
points(t1.obs, t2.obs, pch = 19)
```

### For Saturated Model Mean Value Parameters {#mean-one}

After all of the above, the analogous confidence region for the
saturated model mean value parameter vector $\mu = E(y)$, where
$y$ is the response vector, goes easier (the curve has already
been calculated for $x$ values in the observed data
(the R object `mus`), but we would like to have predictions
for all possible $x$ values, not just those in the observed data).
Thus we are also doing "prediction intervals" for unobserved predictor
values $x$.  These have the form
\begin{equation}
   \mu(x) = \mathop{\rm logit}\nolimits^{-1}(\beta_1 + \beta_2 x),
   (\#eq:mux)
\end{equation}
where $\mathop{\rm logit}^{- 1}$ denotes the inverse of the logit function.

In order for the calculations to not get out of hand, we do calculate
on a grid of values.  Although \@ref(eq:mux) makes sense for any real
$x$, we only cover $x$ values near the observed predictor values
(interpolating not extrapolating).
```{r moo.noo,error=TRUE}
x.new <- seq(-20, 120, length = 1401)
# treat x = 40 and x = 60 specially
x.new <- c(x.new, c(40, 60) + 1e-3, c(40, 60) - 1e3)
x.new <- sort(x.new)
m.new <- cbind(1, x.new)
thetas.new <- m.new %*% t(betas)
thetas.new <- t(thetas.new)
mus.new <- 1 / (1 + exp(- thetas.new))
dim(mus.new)
```
Somewhere along these curves of $\mu$ vectors, maximum and minimum values
are achieved.  As with the $\tau$ vectors, we take an (assumed) limit by
rounding and hoping it is correct, without doing the required real analysis.
```{r moo.noo.max.min}
mu.new.low <- apply(mus.new, 2, min, na.rm = TRUE)
mu.new.hig <- apply(mus.new, 2, max, na.rm = TRUE)
```

So here is the plot (Figure \@ref(fig:muConfidenceRegion)).
```{r muConfidenceRegion,echo=FALSE,fig.align="center",fig.cap="95% Confidence Region (gray) for Saturated Model Mean Value Parameter \\@ref(eq:mux) for Example I.  Dots are as in Figure \\@ref(fig:one)"}
plot(x, y, type = "n",
    xlab = expression(x), ylab = expression(mu(x)),
    xlim = c(-10, 110))
xxx <- c(x.new, rev(x.new))
yyy <- c(mu.new.low, rev(mu.new.hig))
polygon(xxx, yyy, border = NA, col = "gray")
points(x, y)
```

This figure may seem more interpretable than
Figure \@ref(fig:betaConfidenceRegion) or
Figure \@ref(fig:tauConfidenceRegion), but it actually carries
a lot less information.
One-dimensional confidence
intervals do not contain all of the information provided by the
(two-dimensional) confidence region.  Here, even though we have an infinite
number of confidence intervals, one for $\mu(x)$ for each $x$, all of these
intervals together do not tell us what either two-dimensional confidence
region does (Figure \@ref(fig:betaConfidenceRegion) or
Figure \@ref(fig:tauConfidenceRegion)).
Confidence intervals are just less informative than confidence regions
(everywhere in statistics, not just in this context).

Because these confidence intervals come from a
confidence region, they have simultaneous 95\% coverage,
even though there are an infinite number of them.

Another way to say what is wrong with Figure \@ref(fig:muConfidenceRegion)
is the following.  We think the true unknown $\beta$ is in the
confidence region shown in Figure \@ref(fig:betaConfidenceRegion)
(and, of course, are wrong about this 5% of the time).
Any such $\beta$ corresponds to a smooth curve $x \mapsto \mu(x)$
given by \@ref(eq:mux).  So how Figure \@ref(fig:muConfidenceRegion)
should be interpreted is the collection of such smooth curves that
fit inside the gray region are the confidence region (of curves).

# Example II

## Data

@agresti introduces the notion of quasi-complete separation
with the following example, which adds two data points to the data for
Example I.
```{r example.ii}
x <- c(x, 50, 50)
y <- c(y, 0, 1)
```

Figure \@ref(fig:three) shows these data.
```{r label=three,echo=FALSE,fig.align="center",fig.cap="Logistic Regression Data for Example II"}
plot(x, y)
```

## Attempt to Fit the Model

Again we want to do "simple" logistic regression (one predictor $x$ plus
intercept, so the model is two-dimensional).  Again, if we try to do it
naively, the R function `glm` complains.
```{r naive.two}
gout <- glm(y ~ x, family = binomial)
summary(gout)
```

But, again, the warning is useless because R does not give us any help in
doing anything valid with these data.

## Submodel Canonical Statistic

As in Section \@ref(one) we find the support of the submodel canonical
statistic.
```{r example.ii.all.y}
yy <- matrix(0:1, nrow = 2, ncol = length(x))
colnames(yy) <- paste0("y", x)
yy <- expand.grid(as.data.frame(yy))

head(yy)
dim(yy)

m <- cbind(1, x)
mtyy <- t(m) %*% t(yy)
t1 <- mtyy[1, ]
t2 <- mtyy[2, ]
t1.obs <- sum(y)
t2.obs <- sum(x * y)
```
```{r label="four",echo=FALSE,fig.align="center",fig.cap="Possible values of the submodel canonical statistic vector $M^T y$ for the data shown in Figure \\@ref(fig:three).  Solid dot is the observed value of the submodel canonical statistic vector."}
t.key <- paste(t1, t2, sep = ":")
t.idx <- match(unique(t.key), t.key)
t1.uniq <- t1[t.idx]
t2.uniq <- t2[t.idx]

plot(t1.uniq, t2.uniq, xlab = expression(t[1] == sum(y)),
    ylab = expression(t[2] == sum(x * y)))
points(t1.obs, t2.obs, pch = 19)
```

Figure \@ref(fig:four) shows these possible values of the submodel canonical
statistic.

## Computing the Limiting Conditional Model

Again following @geyer-gdor we compute the support of the limiting conditional
model (LCM) and a generic direction of recession (GDOR) as follows
```{r linearity.and.gdor.too}
tanv <- m
tanv[y == 1, ] <- (- tanv[y == 1, ])
vrep <- cbind(0, 0, tanv)
lout <- linearity(vrep, rep = "V")
p <- ncol(tanv)
hrep <- cbind(0, 0, -tanv, -1)
hrep[lout, 1] <- 1
hrep[lout, p + 3] <- 0
hrep <- rbind(hrep, c(0, 1, rep(0, p), -1))
objv <- c(rep(0, p), 1)
pout <- lpcdd(d2q(hrep), d2q(objv), minimize = FALSE)
gdor <- q2d(pout$primal.solution[1:p])
```

As before,
this computes two vectors `lout` and `gdor`.  The first is an index vector
giving the components of the response vector that are random in the LCM
(not conditioned on being equal to their observed values).
```{r show.linearity.too}
lout
```
Here we have found that two components of the response vector (those for
predictor $x = 50$) are random in the LCM.  The LCM is a partially but not
completely degenerate distribution.
The second is a GDOR.
```{r show.gdor.too}
gdor
```
This happens to be the same GDOR as before because the data that are
conditioned on being equal to their observed values in the LCM are the
same in both models.

## Calculating the Maximum Likelihood Estimate

Having found the MLE (for the canonical parameter vector) does not exist,
we find the MLE in the Barndorff-Nielsen completion by finding the MLE
in the LCM.  Following Section 3.14.2 in @geyer-gdor, this is simple.
```{r mle.lcm}
gout.lcm <- glm(y ~ x, family = binomial, subset = lout)
summary(gout.lcm)
```

This says that for the two data points that are treated as random in
the LCM (those for $x = 50$) the MLE has zero intercept and slope,
which corresponds to success probabilty one-half, which is no surprise
because there was one success and one failure for these two data points.

## Confidence Regions

### For Submodel Canonical Parameters

Now we have two kinds of confidence regions.

 * One is our new kind of confidence region, the theory of which
   is expounded in Section \@ref(theory) above.

 * The other is the same old confidence interval calculated from the LCM.

Of course there are many recipies for old-style confidence intervals.
We will use the Wald interval for the beta in the LCM.

```{r wald}
foo <- summary(gout.lcm)$coefficients
old.ci.low <- foo[1, "Estimate"] - qnorm(1 - alpha / 2) * foo[1, "Std. Error"]
old.ci.hig <- foo[1, "Estimate"] + qnorm(1 - alpha / 2) * foo[1, "Std. Error"]
c(old.ci.low, old.ci.hig)
```

Those are the endpoints of our "old" 95% confidence interval for $\beta_1$.
Recall that $\beta_2$ is constrained to be zero in the LCM (it is not unknown).
But we do have a direction of constancy of the LCM.  It is the GDOR of the OM.
(When only one parameter is "not defined because of singularities" as
R function `summary` applied to an object produced by R function `glm` says,
the GDOR is the only direction of constancy of
the LCM.  When there are more such parameters, there are more directions
of constancy of the LCM.)  Hence our confidence region in the original
parameter space consists of all points of the form
$$
   (\beta_1, 0) + s (\eta_1, \eta_2),
$$
where $\eta$ is the GDOR, $s$ is any real number, and $\beta_1$ is in our
"old" confidence interval having endpoints calculated above.
```{r betaConfidenceRegionToo,echo=FALSE,fig.align="center",fig.cap="Confidence Regions for Coefficients for Example II.  New region curved boundary.  Old region straight boundaries."}
idx <- abs(beta1) <= 5 & abs(beta2) <= 5
plot(beta1[idx], beta2[idx], xlab = expression(beta[1]),
    ylab = expression(beta[2]), type = "n")
polygon(beta1, beta2, density = 3)

abline(- old.ci.low * gdor[2] / gdor[1], gdor[2] / gdor[1])
abline(- old.ci.hig * gdor[2] / gdor[1], gdor[2] / gdor[1])
u <- par("usr")
xxx <- c(u[1], old.ci.low + gdor[1] / gdor[2] * u[3], u[2], u[2], u[1])
yyy <- c((u[1] - old.ci.low) * gdor[2] / gdor[1], u[3], u[3],
    (u[2] - old.ci.hig) * gdor[2] / gdor[1],
    (u[1] - old.ci.hig) * gdor[2] / gdor[1])
polygon(xxx, yyy, density = 3, angle = 135)
```

Now that we have two confidence regions (the "new" one that says how close to
to infinity we are in the OM, and the "old" one that restricts the parameters
of the LCM), what do we do with them?

@geyer-gdor suggests combining the two regions using Bonferroni correction.
For example, the intersection of the two regions in the figure would be
a 90% confidence interval.

In order to combine these two confidence regions we need to understand
them a bit better.  The upper boundary of the "old" region is the straight line
$$
   (\beta_1, \beta_2) =
   (u, 0) + s (\eta_1, \eta_2), \qquad s \in \mathbb{R},
$$
where $u$ is the upper bound of the Wald interval (`r old.ci.hig`).  Or 
$$
   (\beta_1 - u, \beta_2) = s (\eta_1, \eta_2), \qquad s \in \mathbb{R},
$$
or
$$
   \frac{\beta_1 - u}{\beta_2} = \frac{\eta_1}{\eta_2}
$$
Since this is an upper bound, we have
$$
   \beta_2 \le \frac{\eta_2}{\eta_1} (\beta_1 - u)
$$
and similarly for the lower bound, so
$$
   \frac{\eta_2}{\eta_1} (\beta_1 - l) \le
   \beta_2 \le \frac{\eta_2}{\eta_1} (\beta_1 - u)
$$
where $l$ is the lower bound ($`r old.ci.low`$)
So now we can calculate the combined region.
```{r combine.regions}
beta2 <- pmin(gdor[2] / gdor[1] * (beta1 - old.ci.hig), beta2)
beta2 <- pmax(gdor[2] / gdor[1] * (beta1 - old.ci.low), beta2)
```

```{r betaConfidenceRegionTooToo,echo=FALSE,fig.align="center",fig.cap="90% Confidence Region for Coefficients for Example II."}
idx <- abs(beta1) <= 5 & abs(beta2) <= 5
plot(beta1[idx], beta2[idx], xlab = expression(beta[1]),
    ylab = expression(beta[2]), type = "n")
polygon(beta1, beta2, density = 3)
```

### For Submodel Mean Value Parameters

We skip this for these data because it won't be that much different
than for the previous data.

### For Saturated Model Mean Value Parameters

Now we follow Section \@ref(mean-one) above except for using the
confidence region shown in Figure \@ref(fig:betaConfidenceRegionTooToo)
rather than the
confidence region shown in Figure \@ref(fig:betaConfidenceRegion).
```{r moo.noo.too,error=TRUE}
# treat x = 50 specially
x.new <- c(x.new, 50 + 1e-3, 50 - 1e3)
x.new <- sort(x.new)
m.new <- cbind(1, x.new)
betas <- cbind(beta1, beta2)
thetas.new <- m.new %*% t(betas)
thetas.new <- t(thetas.new)
mus.new <- 1 / (1 + exp(- thetas.new))
dim(mus.new)
```
This gives all the mus that come from the confidence region, and we find
the ranges as before.
```{r moo.noo.max.min.too}
mu.new.low <- apply(mus.new, 2, min, na.rm = TRUE)
mu.new.hig <- apply(mus.new, 2, max, na.rm = TRUE)
```
And then plot as before.
```{r muConfidenceRegionToo,echo=FALSE,fig.align="center",fig.cap="90% Confidence Region (gray) for Saturated Model Mean Value Parameter \\@ref(eq:mux) for Example II.  Dots are as in Figure \\@ref(fig:three)"}
plot(x, y, type = "n",
    xlab = expression(x), ylab = expression(mu(x)),
    xlim = c(-10, 110))
xxx <- c(x.new, rev(x.new))
yyy <- c(mu.new.low, rev(mu.new.hig))
polygon(xxx, yyy, border = NA, col = "gray")
points(x, y)
```

This is just one of many confidence regions we could have drawn.
There is no reason why we should combine "new" and "old" confidence regions
with equal coverage probabilities.  If they have
coverage probabilities $1 - \alpha_\text{old}$ and $1 - \alpha_\text{new}$,
then Bonferroni says their intersection will have coverage probability
$1 - \alpha_\text{old} - \alpha_\text{new}$.
And we are free to choose $\alpha_\text{old}$ and $\alpha_\text{new}$
as we please.







# Bibliography


