---
title: "Empirical Bayes Deconvolution"
author: "Balasubramanian Narasimhan and Bradley Efron"
date: '`r Sys.Date()`'
output:
  html_document:
  fig_caption: yes
  theme: cerulean
  toc: yes
  toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Empirical Bayes Deconvolution}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r echo=FALSE}
knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    error = FALSE,
    tidy = FALSE,
    cache = FALSE
)
```

## A simulation example

We start with a simulated Poisson example where the $\Theta_i$ are
drawn from a chi-squared density with 10 degrees of freedom and the
$X_i|\Theta_i$ are Poisson with expectation $\Theta_i:$

$$
\Theta_i \sim \chi^2_{10} \mbox{ and } X_i|\Theta_i \sim \mbox{Poisson}(\Theta_i)
$$

The $\Theta_i$ for this setting, with `N = 1000` observations can be
generated as follows.


```{r}
set.seed(238923) ## for reproducibility
N <- 1000
Theta <- rchisq(N,  df = 10)
```

Next, the $X_i|\Theta_i$, for each of `nSIM = 1000` simulations can
be generated as below.

```{r}
nSIM <- 1000
data <- sapply(seq_len(nSIM), function(x) rpois(n = N, lambda = Theta))
```

We take the discrete set $\mathcal{T}=(1, 2, \ldots, 32)$ as the
$\Theta$-space and apply the `deconv` function in the package
`deconvolveR` to estimate $g(\theta).$

```{r}
library(deconvolveR)
tau <- seq(1, 32)
results <- apply(data, 2,
                 function(x) deconv(tau = tau, X = x, ignoreZero = FALSE,
                                    c0 = 1))
```

The default setting for `deconv` uses the `Poisson` family and a
natural cubic spline basis of degree 5 as $Q.$ The regularization
parameter for this example (`c0`) is set to 1.

Some warnings are emitted by the `nlm` routine used for optimization,
but they are mostly inconsequential.

Since `deconv` works on a sample at a time, the result above is a list
of lists from which various statistics can be extracted.  Below, we
construct a table of values for various values of $\Theta$.

```{r}
g <- sapply(results, function(x) x$stats[, "g"])
mean <- apply(g, 1, mean)
SE.g <- sapply(results, function(x) x$stats[, "SE.g"])
sd <- apply(SE.g, 1, mean)
Bias.g <- sapply(results, function(x) x$stats[, "Bias.g"])
bias <- apply(Bias.g, 1, mean)
gTheta <- pchisq(tau, df = 10) - pchisq(c(0, tau[-length(tau)]), df = 10)
gTheta <- gTheta / sum(gTheta)
simData <- data.frame(theta = tau, gTheta = gTheta,
                      Mean = mean, StdDev = sd, Bias = bias,
                      CoefVar = sd / mean)
table1 <- transform(simData,
                    gTheta = 100 * gTheta,
                    Mean = 100 * Mean,
                    StdDev = 100 * StdDev,
                    Bias = 100 * Bias)
```

The table below summarizes the results for some chosen values of
$\theta .$

```{r}
knitr::kable(table1[c(5, 10, 15, 20, 25), ], row.names=FALSE)
```

Although, the coefficient of variation of $\hat{g}(\theta)$ is still
large, the $g(\theta)$ estimates are reasonable.

We can compare the empirical standard deviations and biases of
$g(\hat{\alpha})$ with the approximation given by the formulas in the
paper.

```{r}
library(ggplot2)
library(cowplot)
theme_set(theme_get() +
          theme(panel.grid.major = element_line(colour = "gray90",
                                                size = 0.2),
                panel.grid.minor = element_line(colour = "gray98",
                                                size = 0.5)))
p1 <- ggplot(data = as.data.frame(results[[1]]$stats)) +
    geom_line(mapping = aes(x = theta, y = SE.g), color = "black", linetype = "solid") +
    geom_line(mapping = aes(x = simData$theta, y = simData$StdDev), color = "red", linetype = "dashed") +
    labs(x = expression(theta), y = "Std. Dev")

p2 <- ggplot(data = as.data.frame(results[[1]]$stats)) +
    geom_line(mapping = aes(x = theta, y = Bias.g), color = "black", linetype = "solid") +
    geom_line(mapping = aes(x = simData$theta, y = simData$Bias), color = "red", linetype = "dashed") +
    labs(x = expression(theta), y = "Std. Dev")
plot_grid(plotlist = list(p1, p2), ncol = 2)
```

The approximation is quite good for the standard deviations, but a
little too small for the biases.

## The Shakespeare data

Here we are given the word counts for the entire Shakespeare canon in
the data set `bardWordCount`.  We assume the $i$th distinct word
appeared $X_i \sim Poisson(\Theta_i)$ times in the canon.

```{r}
data(bardWordCount)
str(bardWordCount)
```

We take the support set $\mathcal{T}$ for $\Theta$ to be equally
spaced on the log-scale and the sample space for $\mathcal{X}$ to be
$(1,2,\ldots,100).$

```{r}
lambda <- seq(-4, 4.5, .025)
tau <- exp(lambda)
```

Using a regularization parameter of `c0=2` we can deconvolve the data
to get $\hat{g}.$

```{r}
result <- deconv(tau = tau, y = bardWordCount, n = 100, c0=2)
stats <- result$stats
```

The plot below shows the Empirical Bayes deconvoluation estimates for
the Shakerspeare word counts.

```{r}
ggplot() +
    geom_line(mapping = aes(x = lambda, y = stats[, "g"])) +
    labs(x = expression(log(theta)), y = expression(g(theta)))
```

The quantity $R(\alpha)$ in the paper (Efron, Biometrika 2015) can be
extracted from the `stats` list; in this case for a regularization
parameter of `c0=2` we can print its value:
```{r}
print(result$R)
```
The `stats` list contains other estimates quantities as well.

As noted in the paper citing this package, about
`r 100 * round(stats[161, "G"], 2)` percent of the total mass of
$\hat{g}$ lies below $\Theta = 1$, which is an underestimate. This can
be corrected for by defining
$$
\tilde{g} = c_1\hat{g} / (1 - e^{-\theta_j}),
$$
where $c_1$ is the constant that normalizes $\tilde{g}$.

```{r}
gt <- stats[, "g"] / (1 - exp(-tau))
gt <- gt / sum(gt)
```

```{r}
d <- data.frame(lambda = lambda, g = stats[, "g"], SE.g = stats[, "SE.g"])
indices <- seq(1, length(lambda), 5)
ggplot(data = d) +
    geom_line(mapping = aes(x = lambda, y = g)) +
    geom_errorbar(data = d[indices, ],
                  mapping = aes(x = lambda, ymin = g - SE.g, ymax = g + SE.g),
                  width = .01, color = "blue") +
    labs(x = expression(log(theta)), y = expression(g(theta))) +
    ylim(0, 0.006) +
    geom_line(mapping = aes(x = lambda, y = gt), linetype = "dashed", color = "red")
```

We can now plot the posterior estimates.

```{r, fig.keep='all', fig.width=7.5, fig.height=10}
gPost <- sapply(seq_len(100), function(i) local({gt <- gt * result$P[i, ]; gt / sum(gt)}))
plots <- lapply(c(1, 2, 4, 8), function(i) {
    ggplot() +
        geom_line(mapping = aes(x = tau, y = gPost[, i])) +
        labs(x = expression(theta), y = expression(g(theta)),
             title = sprintf("x = %d", i))
})
plots <- Map(f = function(p, xlim) p + xlim(0, xlim), plots, list(6, 8, 14, 20))
plot_grid(plotlist = plots, ncol = 2)
```

### Bootstrap Comparison

As a check, one can perform a parametric bootstrap using
$$
y^{*} \sim Mult_n(N, {\mathbf f})
$$
with the MLE $\hat{\mathbf f} = {\mathbf f}(\hat{\alpha})$. We use
200 bootstrap replcates to compute the standard errors for $\hat{g}$
and compare them to the theoretical values in the plot below. The
agreement is pretty good.

```{r}
set.seed(1783)
B <- 200
fHat <- as.numeric(result$P %*% d$g)
fHat <- fHat / sum(fHat)
yStar <- rmultinom(n = B, size = sum(bardWordCount), prob = fHat)
gBoot <- apply(yStar, 2,
               function(y) deconv(tau = tau, y = y, n = 100, c0 = 2)$stats[, "g"])
seG <- apply(gBoot, 1, sd)
ggplot(data = d) +
    geom_line(mapping = aes(x = lambda, y = SE.g,
                            color = "Theoretical", linetype = "Theoretical")) +
    geom_line(mapping = aes(x = lambda, y = seG,
                            color = "Bootstrap", linetype = "Bootstrap")) +
    scale_color_manual(name = "Legend",
                       values = c("Bootstrap" = "black", "Theoretical" = "red")) +
    scale_linetype_manual(name = "Legend",
                          values = c("Bootstrap" = "solid", "Theoretical" = "dashed")) +
    theme(legend.title = element_blank()) +
    labs(x = expression(log(theta)), y = expression(sigma(hat(g))))

```

### Predict ratio of new distinct words

Suppose then that a previously unknown Shakespearean corpus of length
$t\times C$ were found, $C\times 900,000$ the length of the known
canon. Assuming a Poisson process model with intensity $\Theta_i$ for
word $i$, the probability that word $i$ did not appear in the canon
but does appear in the new corpus is
$$
e^{-\Theta_i}\left(1-e^{-\Theta_it}\right);
$$
yielding after some work, an estimate for $R(t)$, the expected number of distinct new words found,
divided by $N$, the observed number of distinct words in the canon:
$$
R(t)=\sum_{j=1}^m\hat{g}_jr_j(t),
$$
with
$$
r_j=\frac{e^{-\theta_{(j)}}}{1-e^{-\theta_{(j)}}}\left(1-e^{-\theta_{(j)}t}\right).
$$

We can compute the $R(t)$ and plot it as follows.

```{r}
gHat <- stats[, "g"]
Rfn <- function(t) {
    sum( gHat * (1 - exp(-tau * t)) / (exp(tau) - 1) )
}
r <- sapply(0:10, Rfn)
ggplot() +
    geom_line(mapping = aes(x = 0:10, y = r)) +
    labs(x = "time multiple t", y = expression(R(t)))
```

And the (speculative) doubling time for Shakespeare's vocabulary is
easy to compute too.

```{r}
print(uniroot(f = function(x) Rfn(x) - 1, interval = c(2, 4))$root)
```

## Normal Example

Consider a data set that is generated via the following mechanism.

$$
z_i \sim N(\mu_i, 1), \mbox{ $i = 1,2,\ldots, N = 10,000$}
$$

with

$$
\mu_i=
	\begin{cases}
		0,        & \mbox{with probability .9}\\
	    N(-3, 1), & \mbox{with probability .1}\\
	\end{cases}
$$

```{r}
set.seed(129023)
N <- 10000
pi0 <- .90
data <- local({
    nullCase <- (runif(N) <= pi0)
    muAndZ <- t(sapply(nullCase, function(isNull) {
        if (isNull) {
            mu <- 0
            c(mu, rnorm(1))
        } else {
            mu <- rnorm(1, mean = -3)
            c(mu, rnorm(1, mean = mu))
        }
    }))
    data.frame(nullCase = nullCase, mu = muAndZ[, 1], z = muAndZ[, 2])
})
```

Below is a histogram of the data ($z$ values) and that of the $\Theta$s.

```{r}
p1 <- ggplot(mapping = aes(x = data$z)) +
    geom_histogram(mapping = aes(y  = ..count.. / sum(..count..) ),
                   color = "brown", bins = 60, alpha = 0.5) +
    labs(x = "z", y = "Density")
p2 <- ggplot(mapping = aes(x = data$mu)) +
    geom_histogram(mapping = aes(y  = ..count.. / sum(..count..) ),
                   color = "brown", bins = 60, alpha = 0.5) +
    labs(x = expression(theta), y = "Density")
plot_grid(plotlist = list(p1, p2), ncol = 2)
```

Now we deconvolve this using $\mathcal{T} = (-6, -5.75,\ldots, 3)$ and a
spike at zero and a fifth-degree polynomial.


```{r}
tau <- seq(from = -6, to = 3, by = 0.25)
atomIndex <- which(tau == 0)
result <- deconv(tau = tau, X = data$z, deltaAt = 0, family = "Normal", pDegree = 5)
```

The estimates and the standard errors of the penalized MLE $\hat{g}$
with `c0 = 1` are shown below.

```{r}
knitr::kable(result$stats)
```

Per the above table, the estimated probability of $\mu = 0$ is `r round(result$stats[atomIndex, "g"], 3)`
$\pm$ `r round(result$stats[atomIndex, "SE.g"] , 3)`  with a bias of
about `r round(result$stats[atomIndex, "Bias.g"] , 3)`.

We can now plot the $g$-estimate removing the atom at 0.

```{r}
gData <- as.data.frame(result$stats[-atomIndex, c("theta", "g")])
gData$g <- gData$g / sum(gData$g)
ggplot(data = gData) +
    geom_line(mapping = aes(x = theta, y = g)) +
    geom_line(mapping = aes(x = theta, y = dnorm(theta, mean = -3)),
                            color = "red") +
    labs(x = expression(theta), y = expression(g(theta)))
```

The density approximation is not accurate at all, however, the
posterior estimates for the g's are similar to what one obtains by the
Benjamini-Yekutieli procedure as shown below.

1. Sort the $p$-values.
```{r}
p <- pnorm(data$z)
orderP <- order(p)
p <- p[orderP]
```

2. Compute $R$, the number of discoveries and count the number of
false discoveries
```{r}
## FCR
q <- 0.05
R <- max(which(p <= seq_len(N) * q / N))
discIdx <- orderP[1:R]
disc <- data[discIdx, ]
cat("BY_q procedure discoveries", R, "cases,", sum(disc$nullCase),
    "actual nulls among them.\n")
```

3. Construct Benjamini-Yekutieli and Bayes confidence intervals.

```{r}
alphaR <- 1 - R * q / N
zAlpha <- qnorm(alphaR, lower.tail = FALSE)
zMarker <- max(disc$z)
xlim <- c(-7.6, 0.0)
ylim <- c(-10, 0.0)
BY.lo <- c(xlim[1] - zAlpha, xlim[2] - zAlpha)
BY.up <- c(xlim[1] + zAlpha, xlim[2] + zAlpha)
Bayes.lo <- c(0.5 * (xlim[1] - 3) - 1.96 / sqrt(2), 0.5 * (xlim[2] - 3) - 1.96 / sqrt(2))
Bayes.up <- c(0.5 * (xlim[1] - 3) + 1.96 / sqrt(2), 0.5 * (xlim[2] - 3) + 1.96 / sqrt(2))
```

4. Compute the estimated posterior density for of $\mu$ given $z$ and
   construct the 95% credible intervals for $\mu$ given $z$.

```{r}
d <- data[order(data$mu), ]
muVals <- unique(d$mu)
s <- as.data.frame(result$stats)
indices <- findInterval(muVals, s$theta) + 1
gMu <- s$g[indices]
st <- seq(min(data$z), -2.0, length.out = 40)
gMuPhi <- sapply(st, function(z) gMu * dnorm(z - muVals))
g2 <- apply(gMuPhi, 2, function(x) cumsum(x)/sum(x))
pct <- apply(g2, 2, function(dist) approx(y = muVals, x = dist, xout = c(0.025, 0.975)))
qVals <- sapply(pct, function(item) item$y)
```

5. Plot it

```{r}
ggplot() +
    geom_line(mapping = aes(x = xlim, y = BY.lo), color = "blue") +
    geom_line(mapping = aes(x = xlim, y = BY.up), color = "blue") +
    geom_line(mapping = aes(x = xlim, y = Bayes.lo), color = "magenta",
              linetype = "dashed") +
    geom_line(mapping = aes(x = xlim, y = Bayes.up), color = "magenta",
              linetype = "dashed") +
    geom_point(mapping = aes(x = disc$z, y = disc$mu), color = "red") +
    geom_point(mapping = aes(x = disc$z[disc$nullCase], y = disc$mu[disc$nullCase]),
                             color = "orange") +
    geom_line(mapping = aes(x = rep(zMarker, 2), y = c(-10, 1))) +
    geom_line(mapping = aes(x = st, y = qVals[1, ]), color = "brown") +
    geom_line(mapping = aes(x = st, y = qVals[2, ]), color = "brown") +
    labs(x = "Observed z", y = expression(mu)) +
    annotate("text", x = -1, y = -4.25, label = "BY.lo") +
    annotate("text", x = -1, y = 1.25, label = "BY.up") +
    annotate("text", x = -7.5, y = -6.1, label = "Bayes.lo") +
    annotate("text", x = -7.5, y = -3.4, label = "Bayes.up") +
    annotate("text", x = -2.0, y = -1.75, label = "EB.lo") +
    annotate("text", x = -2.0, y = -3.9, label = "EB.up") +
    annotate("text", x = zMarker, y = 1.25, label = as.character(round(zMarker, 2)))
```

## Another Normal example

In the first normal example, the $\theta$ distribution had a
significant atom at 0 and the rest of the density was smeared around
-3. We now investigate what happens when the $\theta$ distribution is
clearly bimodal. Below is a histogram of the the $\theta$ and
alongside a histogram of the data, generated using $X_i \sim
N(\theta_i, 1)$.

```{r}
p1 <- ggplot(mapping = aes(x = disjointTheta)) +
    geom_histogram(mapping = aes(y  = ..count.. / sum(..count..) ),
                   color = "brown", bins = 60, alpha = 0.5) +
    labs(x = expression(theta), y = "Density")
set.seed (2332)
z <- rnorm(n = length(disjointTheta), mean = disjointTheta)
p2 <- ggplot(mapping = aes(x = z)) +
    geom_histogram(mapping = aes(y  = ..count.. / sum(..count..) ),
                   color = "brown", bins = 60, alpha = 0.5) +
    labs(x = "z", y = "Density")
plot_grid(plotlist = list(p1, p2), ncol = 2)
```

We deconvolve the data, using various values for the spline degrees of
freedom.


```{r, fig.keep='all', fig.width=7.5, fig.height=10}
tau <- seq(from = -4, to = 6, by = 0.2)
plots1 <- lapply(2:7,
                 function(p) {
                     result <- deconv(tau = tau, X = z, family = "Normal", pDegree = p)
                     g <- result$stats[, "g"]
                     ggplot(mapping = aes(x = z)) +
                         geom_histogram(mapping = aes(y = ..count.. / sum(..count..)),
                                        color = "brown", bins = 60, alpha = 0.2) +
                         geom_line(mapping = aes(x = tau, y = g), color = "blue") +
                         labs(x = "z", y = "Density", title = sprintf("DF = %d", p))
                 })
plots2 <- lapply(2:7,
                 function(p) {
                     result <- deconv(tau = tau, X = z, family = "Normal", pDegree = p)
                     g <- result$stats[, "g"]
                     ggplot(mapping = aes(x = disjointTheta)) +
                         geom_histogram(mapping = aes(y = ..count.. / sum(..count..)),
                                        color = "brown", bins = 60, alpha = 0.2) +
                         geom_line(mapping = aes(x = tau, y = g), color = "blue") +
                         labs(x = expression(theta), y = "Density",
                              title = sprintf("DF = %d", p))
                })
plots <- mapply(function(x, y) list(x, y), plots1, plots2)
plot_grid(plotlist = plots, ncol=2)
```

Choosing the degrees of freedom to be 6, we now examine the effect of
regularization parameter $c0$.


```{r, fig.keep='all', fig.width=7.5, fig.height=10}
plots1 <- lapply(c(0.5, 1, 2, 4, 16, 32),
                 function(c0) {
                     result <- deconv(tau = tau, X = z, family = "Normal", pDegree = 6,
                                      c0 = c0)
                     g <- result$stats[, "g"]
                     ggplot(mapping = aes(x = z)) +
                         geom_histogram(mapping = aes(y = ..count.. / sum(..count..)),
                                        color = "brown", bins = 60, alpha = 0.2) +
                         geom_line(mapping = aes(x = tau, y = g), color = "blue") +
                         labs(x = "z", y = "Density", title = sprintf("C0 = %.1f", c0))
                 })
plots2 <- lapply(c(0.5, 1, 2, 4, 16, 32),
                 function(c0) {
                     result <- deconv(tau = tau, X = z, family = "Normal", pDegree = 6,
                                      c0 = c0)
                     g <- result$stats[, "g"]
                     ggplot(mapping = aes(x = disjointTheta)) +
                         geom_histogram(mapping = aes(y = ..count.. / sum(..count..)),
                                        color = "brown", bins = 60, alpha = 0.2) +
                         geom_line(mapping = aes(x = tau, y = g), color = "blue") +
                         labs(x = expression(theta), y = "Density",
                              title = sprintf("C0 = %.1f", c0))
                })
plots <- mapply(function(x, y) list(x, y), plots1, plots2)
plot_grid(plotlist = plots, ncol = 2)
```


## Binomial Example

The dataset `surg` contains data on intestinal surgery on 844 cancer
patients. In the study, surgeons removed _satellite_ nodes for later
testing. The data consists of pairs $(n_i, X_i)$ where $n_i$ is the
number of satellites removed and $X_i$ is the number found to be
malignant among them.

We assume a binomial model with $X_i \sim Binomial(n_i, \theta_i)$
with $\theta_i$ being the probability of any one satellite site being
malignant for the $i$th patient.

We take $\mathcal{T} = (0.01, 0.02,\ldots, 0.09)$, so $m = 99.$ We
take $Q$ to be the default 5-degree natural spline with columns
standardized to mean 0 and sum of squares equal to 1. The penalization
parameter is set to 1. The figure below shows the estimated prior
density of $g(\theta)$.

```{r}
tau <- seq(from = 0.01, to = 0.99, by = 0.01)
result <- deconv(tau = tau, X = surg, family = "Binomial", c0 = 1)
d <- data.frame(result$stats)
indices <- seq(5, 99, 5)
errorX <- tau[indices]
ggplot() +
    geom_line(data = d, mapping = aes(x = tau, y = g)) +
    geom_errorbar(data = d[indices, ],
                  mapping = aes(x = theta, ymin = g - SE.g, ymax = g + SE.g),
                  width = .01, color = "blue") +
    labs(x = expression(theta), y = expression(paste(g(theta), " +/- SE")))
```

The complete table of estimates and standard errors is also available.

```{r}
knitr::kable(d[indices, ], row.names = FALSE)
```

### Bootstrap Comparison

As a check on the estimates of standard error and bias provided by
`deconv`, we compare the results with what we obtain using a
parametric boostrap.

The boostrap is run as follows. For each of 1000 runs, 844 simulated
realizations $\hat{\Theta}^*$ are sampled from density $\hat{g}.$ Each
gave an $X_i^* \sim Binomial(n_i, \Theta_i^*)$ with $n_i$ the $i$th
sample in the original data set. Finally, $\hat{\alpha}$ was computed
using `deconv`.

```{r}
set.seed(32776)
B <- 1000
gHat <- d$g
N <- nrow(surg)

genBootSample <- function() {
    thetaStar <- sample(tau, size = N, replace = TRUE, prob = gHat)
    sStar <- sapply(seq_len(N),
                    function(i)
                        rbinom(n = 1 , size = surg$n[i], prob = thetaStar[i]))
    data.frame(n = surg$n, s = sStar)
}

bootResults <- lapply(seq_len(B),
                      function(k) {
                          surgBoot <- genBootSample()
                          mat <- deconv(tau = tau, X = surgBoot, family = "Binomial",
                                        c0 = 1)$stats
                          mat[, c("g", "Bias.g")]
                      })

gBoot <- sapply(bootResults, function(x) x[, 1])
BiasBoot <- sapply(bootResults, function(x) x[, 2])

indices <- c(seq(1, 99, 11), 99)

table2 <- data.frame(theta = tau,
                     gTheta = round(gHat * 100, 3),
                     sdFormula = round(d$SE.g * 100, 3),
                     sdSimul = round(apply(gBoot, 1, sd) * 100, 3),
                     BiasFormula = round(d$Bias.g * 100, 3),
                     BiasSimul = round(apply(BiasBoot, 1, mean) * 100, 3))[ indices, ]

```

We print out some estimated quantities for comparison.

```{r}
knitr::kable(table2, row.names = FALSE)
```



