of $\hat{\theta}_{M}.$ This solution may be seen by considering the first
order optimality condition (i.e. setting the gradient of
$\log\tilde{p}_{M}(\theta)$ to zero).

###  S7.2 Asymptotic approximation quality

We here show that, in the large sample limit, $\tilde{\theta}^{*}$ provides a
very close approximation of the MAP estimate, $\theta^{*}.$

######  Proposition S7.1 (Asymptotic approximation quality).

Consider Bayesian logistic regression with a Gaussian prior
$\theta\sim\mathcal{N}(0,I_{N})$. Let $x_{1},x_{2},\dots$ be a sequence of
random i.i.d. covariates satisfying $\mathbb{E}[x_{m}x_{m}^{\top}]\succ 0$ and
with bounded third moment, and let $y_{1},y_{2},\dots$ be responses
distributed as in Equation 18. Denote by $X_{M}\mathrel{\mathop{\mathchar
58\relax}}=[x_{1},x_{2},\dots,x_{M}]^{\top}$ and
$Y_{M}\mathrel{\mathop{\mathchar 58\relax}}=[y_{1},y_{2},\dots,y_{M}]^{\top}$
the covariates and labels of the first $M$ data points. Consider the MAP
estimate of $\theta$ after observing $M$ data points,

$\displaystyle\theta^{*}_{M}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname*{arg\,max}_{\theta}p(\theta|Y_{M},X_{M})\text{ and
the approximation }\tilde{\theta}_{M}^{*}\mathrel{\mathop{\mathchar
58\relax}}=\left[I_{N}+\tilde{\Sigma}_{M}\right]^{-1}\hat{\theta}_{M},$ (S33)

where $\hat{\theta}_{M}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname*{arg\,max}_{\theta}p(Y_{M}|X_{M};\theta)$ and
$\tilde{\Sigma}_{M}\mathrel{\mathop{\mathchar
58\relax}}=\left[-\nabla_{\theta}^{2}\log
p(Y_{M}|X_{M};\theta)\big{|}_{\theta=\hat{\theta}_{M}}\right]^{-1}.$ Then
$\|\tilde{\theta}_{M}^{*}-\theta_{M}^{*}\|\in O_{p}(M^{-2}),$ where $O_{p}$
denotes stochastic convergence in probability.

The $O_{p}(M^{-2})$ convergence rate established in S7.1 is very fast in
comparison to the $O_{p}(M^{-\frac{1}{2}})$ convergence rate of the MLE, as
well as to the $O_{p}(M^{-1})$ rate of convergence of the MAP to the posterior
mean. Notably, this asymptotic rate is consistent with rates observed in
simulation (Figure 5(a)).

###### Proof.

We here show that $\|\theta^{*}_{M}-\tilde{\theta}^{*}_{M}\|$ is
$O_{p}(M^{-2})$. Our route to proving this relies on Lemma S7.1 (Trippe et al.
2019, Lemma E.1), which will provide a sequence of bounds on
$\|\theta^{*}_{M}-\tilde{\theta}^{*}_{M}\|$ that depend on the norms of the
gradients of $\log p_{M}(\cdot)$ at $\tilde{\theta}^{*}_{M},$
$c_{M}\mathrel{\mathop{\mathchar 58\relax}}=\|\nabla_{\theta}\log
p_{M}(\tilde{\theta}_{M}^{*})\|,$ and a sequence of strong log-concavity
constants $\alpha_{M}$ for $\log p_{M}(\cdot)$ which hold on the interval
$\\{t\theta^{*}_{M}+(1-t)\tilde{\theta}^{*}_{M}|t\in[0,1]\\}.$ In particular,
Lemma S7.1 provides that
$\|\theta^{*}_{M}-\tilde{\theta}^{*}_{M}\|\leq\frac{c_{M}}{\alpha_{M}}$ and we
obtain the result by showing that $\alpha_{M}$ grows as $\Omega_{p}(M)$ and
$c_{M}$ drops as $O_{p}(M^{-1})$.

We first use Lemma S7.3 to show that the strong log-concavity constants of
$\log p_{M}$ in a neighborhood of radius $\epsilon$ of
$\theta,B_{\epsilon}(\theta)$ grow as $\Omega_{p}(M).$ This allows us to
establish that $\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|$ is
$O_{p}(M^{-1})$ (Lemma S7.4). Since both $\hat{\theta}_{M}$ and
$\theta^{*}_{M}$ converge strongly to $\theta$ under these conditions (see
e.g. Van der Vaart (2000, Theorem 10.10)), the interval
$\\{t\theta^{*}_{M}+(1-t)\tilde{\theta}^{*}_{M}|t\in[0,1]\\}$ is then
contained within $B_{\epsilon}(\theta)$ with probability approaching $1.$
Consequently, the constants of strong log concavity of $\log p_{M}$ on this
interval, which we take as $\alpha_{1},\alpha_{2},\dots,$ must grow as
$\Omega_{p}(M)$ as well.

Now all that remains is to show that $c_{M}$ drops as $O_{p}(M^{-1}).$ Recall
from above that $\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|$ is $O(M^{-1}).$
This fact and the boundedness of the higher derivatives of $\nabla\log p_{M}$
will allow us to use Taylor’s theorem to obtain the desired rate.

However, before proceeding to a more detailed derivation of this rate, we
introduce some additional notation. Let $\phi(y,a)$ denote the GLM mapping
function, such that

$\displaystyle\phi(y,a=x^{\top}\theta)$ $\displaystyle=\log p(y|x,\theta)$
$\displaystyle=-\log(1+\text{exp}\\{-yx^{\top}\theta\\})$

and note that all higher derivatives with respect to $a$ are bounded. In
particular, third derivative satisfies

$\phi^{\prime\prime\prime}(a)\mathrel{\mathop{\mathchar
58\relax}}=\frac{d^{3}}{da^{3}}\phi(y,a)\leq\frac{1}{6\sqrt{3}},$

where we have dropped $y$ as an argument, because these higher derivatives do
not depend on $y$.

We now proceed to derive a stochastic rate of convergence of
$\|\nabla_{\theta}\log p_{M}(\tilde{\theta}_{M}^{*})\|$. We obtain this
through a long derivation involving a series of upper bounds.

$\displaystyle\|\nabla_{\theta}\log p_{M}(\tilde{\theta}_{M}^{*})\|$
$\displaystyle=\|\nabla_{\theta}(\log
p_{M}-\log\tilde{p}_{M})(\tilde{\theta}^{*}_{M})\|$
$\displaystyle=\|\nabla_{\theta}(\log
p_{M}-\log\tilde{p}_{M})(\hat{\theta}_{M})+(\tilde{\theta}_{M}^{*}-\hat{\theta}_{M})^{\top}\nabla_{\theta}^{2}(\log
p_{M}-\log\tilde{p}_{M})(\theta_{M}^{\prime})\|$ $\displaystyle\text{ // By
Taylor's theorem, for some
}\theta_{M}^{\prime}\in\\{t\hat{\theta}_{M}+(1-t)\tilde{\theta}^{*}_{M}|t\in[0,1]\\}$
$\displaystyle=\|(\tilde{\theta}_{M}^{*}-\hat{\theta}_{M})^{\top}\nabla_{\theta}^{2}(\log
p_{M}(\theta_{M}^{\prime})-\log\tilde{p}_{M}(\theta_{M}^{\prime}))\|$
$\displaystyle\text{ // Since
}\nabla_{\theta}\log\tilde{p}_{M}(\hat{\theta})=\nabla_{\theta}\log
p_{M}(\hat{\theta})$
$\displaystyle=\|(\tilde{\theta}_{M}^{*}-\hat{\theta}_{M})^{\top}\left[\nabla_{\theta}^{2}\log
p(Y_{M}|X_{M},\theta_{M}^{\prime})-\nabla_{\theta}^{2}\log
p(Y_{M}|X_{M},\hat{\theta}_{M})\right]\|$ $\displaystyle\text{ // Since
}\log\tilde{p}_{M}\text{ is a second degree approximation defined at
}\hat{\theta}_{M}$
$\displaystyle\leq\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|\left[\sum_{m=1}^{M}\|\nabla_{\theta}^{2}\log
p(y_{m}|x_{m},\theta^{\prime}_{M})-\nabla_{\theta}^{2}\log
p(y_{m}|x_{m},\hat{\theta}_{M})\|_{\mathrm{OP}}\right]$
$\displaystyle=\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|\left[\sum_{m=1}^{M}\|\theta_{M}^{\prime}-\hat{\theta}_{M}\|\cdot\|\int_{t=0}^{1}\frac{\partial}{\partial
t}\nabla_{\theta}^{2}\log
p(y_{m}|x_{m},\theta)\big{|}_{\theta=t\hat{\theta}_{M}+(1-t)\theta_{M}^{\prime}}\|_{\mathrm{OP}}\right]$
// By the fundamental theorem of calculus
$\displaystyle\leq\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|^{2}\left[\sum_{m=1}^{M}\|\int_{t=0}^{1}\frac{\partial}{\partial
t}\nabla_{\theta}^{2}\log
p(y_{m}|x_{m},\theta)\big{|}_{\theta=t\hat{\theta}_{M}+(1-t)\theta_{M}^{\prime}}\|_{\mathrm{OP}}\right]$
$\displaystyle\leq\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|\left[\sum_{m=1}^{M}\|x_{m}\|^{3}(\text{max}_{a}\phi^{\prime\prime\prime}(a))\right]$
$\displaystyle=\frac{1}{6\sqrt{3}}\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|^{2}\left[\sum_{m=1}^{M}\|x_{m}\|^{3}\right]$
$\displaystyle\leq O_{p}(\frac{1}{M^{2}})O_{p}(M)=O_{p}(\frac{1}{M}),$

where the final line requires that the covariates have bounded third moment. ∎

#### Supporting Lemmas

######  Lemma S7.1 (Trippe et al., 2019, Lemma E.1).

Let $f,g$ be twice differentiable functions mapping
$\mathbb{R}^{N}\rightarrow\mathbb{R}$ and attaining minima at
$\theta_{f}=\operatorname*{arg\,min}_{\theta}f(\theta)$ and
$\theta_{g}=\operatorname*{arg\,min}_{\theta}g(\theta)$, respectively.
Additionally, assume that $f$ is $\alpha$–strongly convex for some $\alpha>0$
on the set $\\{t\theta_{f}+(1-t)\theta_{g}|t\in[0,1]\\}$ and that
$\|\nabla_{\theta}f(\theta_{g})-\nabla_{\theta}g(\theta_{g})\|_{2}=\|\nabla_{\theta}f(\theta_{g})\|_{2}\leq
c$. Then

$\displaystyle\|\theta_{f}-\theta_{g}\|_{2}\leq\frac{c}{\alpha}.$ (S34)

######  Lemma S7.2 (uniform law of large numbers).

Let $H_{M}(\theta)$ be as defined in Equation S32 and define
$H(\theta)\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}[\nabla_{\theta}^{2}\log p(y_{1}|x_{1};\theta)]$, where
the expectation is taken under the true $\theta$. If
$\mathbb{E}[x_{1}x_{1}^{\top}]$ exists and is positive definite then

$\sup_{\theta^{\prime}\in
B_{\epsilon}(\theta)}\|\frac{1}{M}H_{M}(\theta^{\prime})-H(\theta^{\prime})\|_{2}\overset{a.s.}{\rightarrow}0.$

according to $p$, where $B_{\epsilon}(\theta)$ is a closed neighborhood of
$\theta$ of radius $\epsilon,$ for any $\epsilon>0.$

###### Proof.

Since the each of the $M$ data points $\\{(x_{m},y_{m})\\}_{m=1}^{\infty}$ are
i.i.d. by assumption, $M^{-1}H_{M}$ converges point-wise by the law of large
numbers. However, we are additionally interested in uniform convergence; a
number of different uniform laws of large numbers suffice for this. Because
$H$ is continuously differentiable in $\theta$ (recall that for any $x_{m}$,
$\frac{d^{3}}{d\theta^{3}}\log p(y_{m}|x_{m},\theta)$ is bounded) it is
therefore Lipschitz continuous on the bounded set $B_{\epsilon}(\theta)$. As
such one can construct a bounded envelope for $H$ on this set, which amounts
to a sufficient condition for uniform convergence on $B_{\epsilon}$, see Van
der Vaart (2000, Theorem 19.4 - Glivenko-Cantelli). We refer the reader to Van
der Vaart (2000, Chapter 19) for technical background, and in particular to
Van der Vaart (2000, Example 19.8) which walks through an example closely
related to the present case. ∎

######  Lemma S7.3.

Consider logistic regression with random covariates, $x_{1},x_{2},\dots.$ Let
$B_{\epsilon}(\theta)$ be a closed neighborhood of radius $\epsilon>0$ around
$\theta$ and for each $M$ define

$\alpha_{M}\mathrel{\mathop{\mathchar 58\relax}}=\inf_{\theta^{\prime}\in
B_{\epsilon}(\theta)}\lambda_{min}\left[\nabla_{\theta}^{2}\log
p_{M}(\theta^{\prime})\right]$

to be the constant of strong log-concavity constant of $\log p_{M}(\cdot)$ on
$B_{\epsilon}(\theta),$ where $\lambda_{min}(\cdot)$ denotes the smallest
eigenvalue of its matrix argument. If the covariates are i.i.d. and satisfy
$\mathbb{E}[x_{1}x_{1}^{\top}]\succ 0,$ then $\alpha_{M}$ is $\Omega_{p}(M).$

###### Proof.

Consider the scaled Hessians of $\log p_{M}(\cdot)$, $M^{-1}H_{M}(\cdot).$ By
Lemma S7.2, $M^{-1}H_{M}(\cdot)$ converges uniformly to its expectation,
$H(\theta)\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}[\nabla_{\theta}^{2}\log p(y_{1}|x_{1},\theta)]$ on
$B_{\epsilon}(\theta).$ Since $H(\theta)\succ 0$ on $B_{\epsilon}(\theta),$ we
have that

$\inf_{\theta^{\prime}\in
B_{\epsilon}(\theta)}\lambda_{min}(\frac{1}{M}H_{M}(\theta))\overset{a.s.}{\rightarrow}\inf_{\theta^{\prime}\in
B_{\epsilon}(\theta)}\lambda_{min}(H_{M}(\theta))>0.$

Therefore $\alpha_{M}\mathrel{\mathop{\mathchar
58\relax}}=\inf_{\theta^{\prime}\in
B_{\epsilon}(\theta)}\lambda_{min}\left[\nabla_{\theta}^{2}\log
p_{M}(\theta^{\prime})\right]$ is $\Omega_{p}(M).$ ∎

######  Lemma S7.4.

Let $\hat{\theta}$ and $\tilde{\theta}^{*}$ be the MLE and the approximation
to the MAP defined in Equation 19, respectively. If the covariates,
$x_{1},x_{2},\dots$ are i.i.d. and satisfy $\mathbb{E}[x_{1}x_{1}^{\top}]\succ
0,$ then $\|\hat{\theta}_{M}-\tilde{\theta}^{*}_{M}\|$ is $O_{p}(M^{-1}).$

###### Proof.

Recall that

$\tilde{\theta}^{*}_{M}=\left[I_{N}+\tilde{\Sigma}_{M}\right]^{-1}\hat{\theta}_{M},$

where $\tilde{\Sigma}_{M}\mathrel{\mathop{\mathchar
58\relax}}=H_{M}(\hat{\theta}_{M})^{-1}.$ Lemma S7.3 provides that the
constants of strong log-concavity for $\log p_{M}$ grow as $\Omega_{p}(M)$ in
a neighborhood of $\theta.$ Therefore, since $\hat{\theta}_{M}$ converges
strongly to $\theta,$ we can see that $\lambda_{min}(H_{M}(\hat{\theta}_{M}))$
is $\Omega_{p}(M).$ Next, we rewrite

$\displaystyle\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|$
$\displaystyle=\|\left[I_{N}+\tilde{\Sigma}_{M}\right]^{-1}\hat{\theta}_{M}-\hat{\theta}_{M}\|$
$\displaystyle=\|\left[I_{N}+H_{M}(\hat{\theta}_{M})\right]^{-1}\hat{\theta}_{M}\|$
$\displaystyle\leq\|\left[I_{N}+H_{M}(\hat{\theta}_{M})\right]^{-1}\|_{\mathrm{OP}}\|\hat{\theta}_{M}\|$
$\displaystyle\leq\frac{\|\hat{\theta}_{M}\|}{\lambda_{min}\left(H_{M}(\hat{\theta}_{M})\right)}.$

which one can see is $O_{p}(M^{-1})$ since $\|\hat{\theta}_{M}\|$ is bounded
in probability. ∎

###  S7.3 Proof of Theorem 5.3

Before proving the theorem we begin by explicitly writing out the win and our
proposed bound defined in Section 5.2. For clarity, we introduce a subscript
$M$ to index the size of the dataset on which these quantities are computed.
Specifically, recalling that in this case we have $A=I_{N}$ and
$C=(I_{N}+\tilde{\Sigma}_{M})^{-1},$ and noting that therefore
$A-C=(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1},$ we have

$\displaystyle b_{M}(\alpha)$
$\displaystyle=2\mathrm{tr}[(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}]+{}$
$\displaystyle
2z_{\frac{1-\alpha}{2}}\sqrt{U_{M}(\|G_{M}(\hat{\theta}_{M})\|_{\tilde{\Sigma}_{M}}^{2},\frac{1-\alpha}{2})+2\|\tilde{\Sigma}_{M}^{\frac{1}{2}}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}^{\frac{1}{2}}\|_{F}^{2}}-\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|^{2}$

where $G_{M}(\hat{\theta}_{M})\mathrel{\mathop{\mathchar
58\relax}}=(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\hat{\theta}_{M}$ and

$\displaystyle U_{M}$
$\displaystyle(\|G_{M}(\hat{\theta}_{M})\|_{\tilde{\Sigma}_{M}}^{2},1-\alpha)\mathrel{\mathop{\mathchar
58\relax}}=\inf_{\delta>0}\bigg{\\{}\delta\,\bigg{|}\,\|G_{M}(\hat{\theta}_{M})\|_{\tilde{\Sigma}_{M}}^{2}\leq(\delta+\|\tilde{\Sigma}_{M}^{\frac{1}{2}}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}^{\frac{1}{2}}\|_{F}^{2})+{}$
(S35) $\displaystyle
z_{1-\alpha}\sqrt{2\|\tilde{\Sigma}_{M}^{\frac{1}{2}}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}^{\frac{1}{2}}\|_{F}^{2}+4\|\tilde{\Sigma}_{M}^{\frac{1}{2}}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}^{\frac{1}{2}}\|_{\mathrm{OP}}^{2}\delta}\bigg{\\}}$
(S36)

is an approximate high-confidence upper bound on
$\|G_{M}(\hat{\theta}_{M})\|_{\tilde{\Sigma}_{M}}^{2}.$ For convenience, we
abbreviate
$U_{M}(\|G_{M}(\hat{\theta}_{M})\|_{\tilde{\Sigma}_{M}}^{2},1-\alpha)$ by
$U_{M}.$

Next, we recall that we may decompose the win in squared error loss for using
$\theta^{*}_{M}$ in place of $\hat{\theta}_{M}$ as

$W_{M}(\theta)=2\epsilon_{M}^{\top}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\hat{\theta}-\|\theta^{*}_{M}-\hat{\theta}_{M}\|^{2},$

where $\epsilon_{M}\mathrel{\mathop{\mathchar
58\relax}}=\hat{\theta}_{M}-\theta.$

###### Proof.

Proving the theorem amounts to showing that for any $\theta$ and
$\alpha\in(0,1)$,

$\lim_{M\rightarrow\infty}\mathbb{P}_{\theta}\left[W_{M}(\theta)\geq
b_{M}(\alpha)\right]\geq\alpha.$

Lemma S7.6 provides that $M^{1.5}(W_{M}(\theta)-b_{M}(\alpha))$ converges in
distribution to $2\sqrt{\theta^{\top}H(\theta)^{-3}\theta}(\delta-
z_{\frac{1-\alpha}{2}}),$ for $\delta\sim\mathcal{N}(0,1).$ Thus for any
$\theta,$
$\mathbb{P}_{\theta}\left[W_{M}(\theta)-b_{M}(\alpha)>0\right]\rightarrow(1-\Phi(z_{\frac{1-\alpha}{1}}))=1-\frac{1-\alpha}{2}>\alpha.$
This establishes that $b_{M}(\cdot)$ has above nominal coverage
asymptotically, as desired. ∎

######  Lemma S7.5.

$|U_{M}-\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}|$ is
$O_{p}(M^{-3.5}).$

###### Proof.

Recall that we can rearrange Equation S35 to see that $U_{M}$ satisfies

$\displaystyle\|(I+\tilde{\Sigma}_{M}^{-1})^{-1}\hat{\theta}_{M}\|_{\tilde{\Sigma}_{M}}^{2}$
$\displaystyle=U_{M}+2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{F}^{2}+{}$
$\displaystyle\sqrt{\|\tilde{\Sigma}_{M}^{4}(I_{N}+\tilde{\Sigma}_{M})^{2}\|_{F}^{2}+4\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}^{2}U_{M}}$

where we have simplified
$\tilde{\Sigma}_{M}^{\frac{1}{2}}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}^{\frac{1}{2}}$
to $\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}.$

We next further simplify the condition above by replacing two quantities with
simplifying approximations plus lower order terms. First note that we may
write

$\displaystyle\|(I+\tilde{\Sigma}_{M}^{-1})^{-1}\hat{\theta}_{M}\|_{\tilde{\Sigma}_{M}}^{2}$
$\displaystyle=\|\tilde{\Sigma}_{M}\hat{\theta}_{M}-\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\hat{\theta}_{M}\|_{\tilde{\Sigma}_{M}}^{2}$
$\displaystyle=\|\tilde{\Sigma}_{M}\hat{\theta}_{M}\|_{\tilde{\Sigma}_{M}}^{2}+\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\hat{\theta}_{M}\|_{\tilde{\Sigma}_{M}}^{2}-2\hat{\theta}_{M}^{\top}\tilde{\Sigma}_{M}^{4}(I_{N}+\tilde{\Sigma}_{M})^{-1}\hat{\theta}_{M}$
$\displaystyle=\|\tilde{\Sigma}_{M}(\theta+\epsilon_{M})\|_{\tilde{\Sigma}_{M}}^{2}+O_{p}(M^{-4})$
$\displaystyle=\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}+\|\tilde{\Sigma}_{M}\epsilon_{M}\|_{\tilde{\Sigma}_{M}}^{2}+2\epsilon_{M}^{\top}\tilde{\Sigma}_{M}^{3}\theta+O_{p}(M^{-4})$
$\displaystyle=\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}+O_{p}(M^{-3.5}).$

Second, we write

$\displaystyle\sqrt{\|\tilde{\Sigma}_{M}^{4}(I_{N}+\tilde{\Sigma}_{M})^{2}\|_{F}^{2}+4\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}^{2}U_{M}}$
$\displaystyle=\sqrt{O_{p}(M^{-8})+4\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}^{2}U_{M}}$
$\displaystyle=2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}\sqrt{U_{M}}+O_{p}(M^{-4}).$

As such, we may see that $U_{M}$ satisfies

$\displaystyle\begin{split}\|\tilde{\Sigma}_{M}\theta_{M}\|_{\tilde{\Sigma}_{M}}^{2}-U_{M}&=2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{F}^{2}+2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}\sqrt{U_{M}}+O_{p}(M^{-3.5})\\\
&=2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}\sqrt{U_{M}}+O_{p}(M^{-3.5})\end{split}$
(S37)

where we have dropped
$2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{F}^{2}$ since it
is $O_{p}(M^{-4}).$

We next observe that $U_{M}$ must be $O_{p}(M^{-3}).$ Otherwise, the event
that $\|\tilde{\Sigma}_{M}\theta_{M}\|^{2}_{\tilde{\Sigma}_{M}}-U_{M}<0$ must
occur infinitely often (since
$\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}$ is $O_{p}(M^{-3})$);
in turn, this condition would imply that
$\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{\mathrm{OP}}\sqrt{U_{M}}<0$
occurs infinitely often, which provides a contradiction.

Finally, in tangent with Equation S37, that $U_{M}$ is $O_{p}(M^{-3})$ allows
us to see that $\left|U_{M}-\|\Sigma\theta\|^{2}_{\tilde{\Sigma}_{M}}\right|$
is $O_{p}(M^{-3.5})$, as desired. ∎

######  Lemma S7.6.

Let $\alpha\in(0,1)$ and $\theta\in\mathbb{R}^{N}.$ Consider the sequence of
wins, $W_{M}(\theta)$, and bounds, $b_{M}(\alpha),$ computed for logistic
regression. Then

$M^{1.5}(W_{M}(\theta)-b_{M}(\alpha))\overset{d}{\rightarrow}2\sqrt{\theta^{\top}H(\theta)^{-3}\theta}(\delta-
z_{\frac{1-\alpha}{2}}),$

where $\delta\sim\mathcal{N}(0,1).$

###### Proof.

We prove the lemma by first writing $W_{M}$ and $b_{M}$ using simplifying
approximations and lower order terms. The result is obtained by manipulating a
scaling of the difference between the two expressions and considering the
limit in $M.$

Note first that we may write

$\displaystyle W_{M}(\theta)\mathrel{\mathop{\mathchar 58\relax}}$
$\displaystyle=2\epsilon^{\top}(\theta^{*}_{M}-\hat{\theta}_{M})-\|\theta^{*}_{M}-\hat{\theta}_{M}\|^{2}$
$\displaystyle=2\epsilon^{\top}(\tilde{\theta}^{*}_{M}-\hat{\theta}_{M})-\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2})$
$\displaystyle=2\epsilon^{\top}(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\hat{\theta}_{M}-\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2})$
$\displaystyle=2\epsilon^{\top}\tilde{\Sigma}_{M}\hat{\theta}_{M}-\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2})$
$\displaystyle=2\epsilon^{\top}\tilde{\Sigma}_{M}\theta-\|\tilde{\theta}^{*}_{M}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2}).$

Next we write

$\displaystyle b_{M}(\alpha)$
$\displaystyle=2\mathrm{tr}\left[(I_{N}+\tilde{\Sigma}_{M}^{-1})^{-1}\tilde{\Sigma}_{M}\right]+2z_{\frac{1-\alpha}{2}}\sqrt{U_{M}+2\|\tilde{\Sigma}_{M}^{2}(I_{N}+\tilde{\Sigma}_{M})^{-1}\|_{F}^{2}}-\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|^{2}$
$\displaystyle=2z_{\frac{1-\alpha}{2}}\sqrt{\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}+O_{p}(M^{-3.5})}-\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2})$
$\displaystyle=2z_{\frac{1-\alpha}{2}}\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}-\|\tilde{\theta}_{M}^{*}-\hat{\theta}_{M}\|^{2}+O_{p}(M^{-2}).$

where the second line uses Lemma S7.5.

By considering a scaled difference between these two terms we find,

$\displaystyle M^{1.5}(W_{M}(\theta)-b_{M}(\alpha))$
$\displaystyle=2M^{1.5}\epsilon^{\top}\tilde{\Sigma}_{M}\theta-2M^{1.5}z_{\frac{1-\alpha}{2}}\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}+O_{p}(M^{-\frac{1}{2}})$
$\displaystyle\overset{d}{\rightarrow}2M^{1.5}\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}(\delta-
z_{\frac{1-\alpha}{2}})$

for $\delta\sim\mathcal{N}(0,1),$ by recognizing that $\epsilon_{M}$ is
asymptotically normal with mean zero and covariance $\Sigma_{M},$ and
therefore that $2\epsilon^{\top}\tilde{\Sigma}_{M}\theta$ is asymptotically
normal with variance $\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}^{2}.$

Finally, the result obtains by noting that Lemma S7.2 implies that

$M^{1.5}\|\tilde{\Sigma}_{M}\theta\|_{\tilde{\Sigma}_{M}}=\sqrt{\theta^{\top}(H_{M}(\hat{\theta})/M)^{-3}\theta}\overset{a.s.}{\rightarrow}\sqrt{\theta^{\top}H(\theta)^{-3}\theta}.$

∎

###  S7.4 Empirical validation of logistic regression bound in simulation

(a)

(b)

(c)

Figure S5: c-values for logistic regression in two sets of simulations. With
$N=2,$ (a) empirical rates of convergence of distances amongst various
estimates and the true parameter. With $N=25$ and $M=1000$ (b) c-values are
able to detect improvements, sometimes with high confidence, and (c) the
approximate bound has greater than nominal coverage. See Section S7.4 for
details.

We here demonstrate the fast convergence of our approximation to the MAP in
logistic regression on simulated data. We also include supplementary results
illustrating the favorable performance of c-values in this setting, which is
made possible by this fast convergence. Figure 5(a) shows the distance between
various estimates and the true parameter for a range of sample sizes in
simulation. Due to the log-log scale, the slopes of the series in this plot
reflect the polynomial rates of convergence. Notably we see the fast
$O_{p}(M^{-2})$ rate of convergence of our approximation to the MAP estimate,
$\tilde{\theta}^{*}_{M},$ to the exact MAP estimate, $\theta^{*}_{M}.$

Figure 5(b) demonstrates that our approach is able to detect improvements
(i.e. we can obtain high c-values). Furthermore, our proposed bound has
similar coverage properties as in the Gaussian case (Figure 5(c)). In the
experiments for Figures 5(b) and 5(c), we simulated the parameter as
$\theta\sim\mathcal{N}(0,\frac{1}{2}I_{N})$ and, in each replicate, simulated
the covariates for each data point, indexed by $m,$ as
$x_{m}\overset{i.i.d.}{\sim}\mathcal{N}(0,N^{-2}I_{N}).$

Two of the series in Figure 5(a) are distances between the posterior mean of
$\theta$ and other estimates, $\mathbb{E}[\theta|X,Y]=\int p(\theta|X,Y)\theta
d\theta.$ Because this model is non-conjugate, the estimate does not have an
analytic form. As such we approximated these quantities with Gauss-Hermite
quadrature. For each sample size $M,$ we performed $25$ replicate simulations.

In the experiments that went into Figures 5(b) and 5(c), we used $N=25$ and
$M=1000.$ See logistic_regression_approximations.ipynb and
logistic_regression_c_values_and_operating_characteristics.ipynb for details.

##  Appendix S8 Additional details on applications

In this section, we provide additional details associated with the
applications in Section 6.

###  S8.1 Estimation from educational testing data

#### Conservatism of c-values with the empirical Bayes step.

The application in Section 6.1 diverges from the scenarios covered by our
theory in Sections 3 and 4 in its use of the empirical Bayes step to estimate
$\beta,\tau,$ and $\sigma.$ As a result, our theory does provide that $c(y)$
satisfies the guarantee of Theorem 2.2. However, given the favorable
asymptotic and empirical properties of the empirical Bayes procedure
established in Section 5.1, we conjectured that the looseness in the lower
bound $b(y,\alpha)$ would be sufficiently large to compensate for any error
introduced by these departures from the assumptions of our theory.

Figure S6:  Calibration of the lower bounds $b(y,\alpha)$ in small area
inference with an empirical Bayes step ($5000$ replicates). The coverage on
the y-axis is a Monte Carlo estimate of
$\mathbb{P}_{\theta}\left[W(\theta,y)\geq b(y,\alpha)\right].$ Each series
corresponds to a set of simulations within which we excluded a different
subset of schools based on a minimum number of students tested.

To investigate this, we performed a simulation study in which we used this
empirical Bayes step and confirmed that the c-values retained at least nominal
coverage (Figure S6). To ensure that the simulated data had similar
characteristics to the real data, we simulated $5000$ datasets by drawing
hypothetical school level means according the assumed generative model with
the parameters ($\beta,\tau$ and $\sigma$) fit on the real dataset. In each
simulation, we re-estimated the fixed effects and variances (again using
lme4), and computed the associated MLE, Bayes estimates, and bounds across a
range of confidence levels. We then computed the empirical coverage of these
bounds and found them to be conservative across all tested levels.

#### Additional preprocessing and calibration details.

Hoff (2021) considered only schools at which $2$ or more students took the
reading test. We excluded an additional $8$ schools with fewer than $5$
students tested because we expected that the high variance in these
observations could introduce too much slack into our bound as result of the
poor conditioning of $\Sigma^{\frac{1}{2}}(A-C)\Sigma^{\frac{1}{2}}$ (recall
the operator norm bound in Equation 11, derived in Equation S26). Consistent
with this hypothesis we computed a c-value of $0.88$ when we included these
additional schools, and when we further restricted to the $657$ schools with
at least $10$ students tested we computed a c-value $0.999992.$ To further
validate this hypothesis of increased conservatism we simulated additional
datasets with these different thresholds on school size and evaluated the
calibration of computed bounds (Figure S6). We observed the coverage for the
simulations with smallest threshold was noticeably higher at large $\alpha$,
in agreement with this hypothesis.

###  S8.2 Estimation of violent crime rates in Philadelphia

#### Dependence on the order in which estimates are compared.

In Section 6.2 we chose to report one among three estimates as described in
Remark 6.1. We note however that this paradigm is sensitive to the order in
which the different estimates are considered. For this set of three models, if
we had first compared $\theta^{\circ}(y)$ as the alternative to
$\hat{\theta}(y)$ as the default we would have rejected $\hat{\theta}(y)$
(with $c=0.99942$), and then again sided against updating our estimate a
second time with a low c-value ($c=0.0$) for comparing $\theta^{*}(y)$ as the
alternative against $\theta^{\circ}(y)$ as the default. The potential cost of
ending up with a worse estimate as a result of considering these estimates in
sequence may be understood as a cost of looking at the data an additional
time.

#### Selection of prior parameters from historical data.

The parameters $\sigma_{\delta}^{2},\sigma_{z}^{2},\sigma_{y}^{2}$ were
selected based on historical data. Specifically, we estimated $\sigma_{y}^{2}$
and $\sigma_{z}^{2}$ as the averages of the sample variances of the violent
and non-violent report rates, respectively, computed within each census block
in the preceding years. For the first model described in Section 6.2, we then
estimated $\sigma_{\delta}^{2}$ using these same historical data to reflect
the prior belief that half of the variability across the unknown rates is
common across the two response types.

For the second model considered, we selected the signal variance and length
scale of this covariance function by drawing hypothetical datasets of crime
levels from the prior predictive distributions and selecting those which
produced the most reasonable looking patterns. In particular, we chose the
length scale to be one sixth of the maximum distance between the centroids of
census blocks, and the signal variance to reflect the prior belief that one
third of the variability in the unknown rates was explained by the spatial
component. In addition, we choose a smaller value for $\sigma_{\delta}^{2}$ in
this second model, so that the total implied variance would be the same. See
supplementary code in Philly_reported_crime_estimation.ipynb for additional
details.

#### Derivation of $\theta^{*}$ (posterior mean in the first model).

As mentioned in the main text, since the prior and likelihoods for this model
are independent across each census block we can compute the posterior mean for
each block independently.

Let $\pi(\cdot)$ denote the joint density of all variables. Then, since
$z_{n}\raisebox{0.50003pt}{\rotatebox[origin={c}]{90.0}{$\models$}}y_{n}\big{|}\theta_{n}$,
we have that

$\displaystyle\pi(\theta_{n}|y_{n},z_{n})$
$\displaystyle\propto\pi(\theta_{n}|z_{n})\pi(y_{n}|\theta_{n},z_{n})$
$\displaystyle=\pi(\theta_{n}|z_{n})\pi(y_{n}|\theta_{n}).$

Next observe that by construction,
$z_{n}-\theta_{n}=\epsilon^{z}_{n}+\delta_{n}^{z}-\delta_{n}^{y}\sim\mathcal{N}(0,2\sigma_{\delta}^{2}+\sigma_{z}^{2})$
and so
$\theta_{n}|z_{n}\sim\mathcal{N}(z_{n},2\sigma_{\delta}^{2}+\sigma_{z}^{2})$.
Since again by construction we have that
$y_{n}|\theta_{n}\sim\mathcal{N}(\theta_{n},\sigma_{y}^{2}),$ Gaussian
conjugacy provides that

$\theta_{n}|y_{n},z_{n}\sim\mathcal{N}(\mathbb{E}[\theta_{n}|y_{n},z_{n}],\mathrm{Var}[\theta_{n}|y_{n},z_{n}]),$

where

$\displaystyle\mathrm{Var}[\theta_{n}|y_{n},z_{n}]$
$\displaystyle=\frac{1}{\sigma_{y}^{-2}+(2\sigma_{\delta}^{2}+\sigma_{z}^{2})^{-1}}$
$\displaystyle=\frac{\sigma_{y}^{2}(2\sigma_{\delta}^{2}+\sigma_{z}^{2})}{\sigma_{y}^{2}+2\sigma_{\delta}^{2}+\sigma_{z}^{2}}$

and

$\displaystyle\mathbb{E}[\theta_{n}|y_{n},z_{n}]$
$\displaystyle=\mathrm{Var}[\theta_{n}|y_{n},z_{n}](\mathrm{Var}[\theta_{n}|z_{n}]^{-1}\mathbb{E}[\theta_{n}|z_{n}]+\mathrm{Var}[y_{n}|\theta_{n}]^{-1}y_{n})$
$\displaystyle=\frac{\sigma_{y}^{2}(2\sigma_{\delta}^{2}+\sigma_{z}^{2})}{\sigma_{y}^{2}+2\sigma_{\delta}^{2}+\sigma_{z}^{2}}\left[(2\sigma_{\delta}^{2}+\sigma_{z}^{2})^{-1}z_{n}+\sigma^{-2}y_{n}\right]$
$\displaystyle=\frac{2\sigma_{\delta}^{2}+\sigma_{z}^{2}}{2\sigma_{\delta}^{2}+\sigma_{y}^{2}+\sigma_{z}^{2}}y_{n}+\frac{\sigma_{y}^{2}}{2\sigma_{\delta}^{2}+\sigma_{y}^{2}+\sigma_{z}^{2}}z_{n}$

as desired.

Analogously, for the second model considered in Section 6.2 we find the
posterior mean as

$\theta^{\circ}(y)=\left[I_{N}+\sigma_{y}^{2}(2K+2\sigma_{\delta}^{2}I_{N}+\sigma_{z}^{2}I_{N})^{-1}\right]^{-1}y+\left[I_{N}+\sigma_{y}^{-2}(2K+2\sigma_{\delta}^{2}I_{N}+\sigma_{z}^{2}I_{N})\right]^{-1}z.$

#### Additional dataset details.

The data considered in this application are counts of police responses
categorized as associated with violent crimes and violent crimes in October
2018. These were obtained from opendataphilly.org. The observed data we model
are the inverse hyperbolic sine transform of the number of recorded police
responses per square mile. For all practical purposes, these values can be
interpreted as log densities (see, e.g., Burbidge et al. (1988)).

###  S8.3 Gaussian process kernel selection for estimation of ocean currents

We here provide additional details of the Gaussian process covariance
functions used in Section 6.3. The first covariance function described, which
incorporated covariation at two scale is defined, for both the longitudinal
and latitudinal components ($i$ in $\\{1,2\\}$) and for each pair of buoys $n$
and $n^{\prime}$, as

$\displaystyle k(\theta^{(i)}_{n},\theta^{(i)}_{n^{\prime}})=$
$\displaystyle\sigma_{1}^{2}\exp\left\\{-\frac{1}{2}\left[\frac{(\text{lat}_{n}-\text{lat}_{n^{\prime}})^{2}}{r_{1,\text{lat}}^{2}}+\frac{(\text{lon}_{n}-\text{lon}_{n^{\prime}})^{2}}{r_{1,\text{lon}}^{2}}+\frac{(t_{n}-t_{n^{\prime}})^{2}}{r_{1,t}^{2}}\right]\right\\}$
$\displaystyle+$
$\displaystyle\sigma_{2}^{2}\exp\left\\{-\frac{1}{2}\left[\frac{(\text{lat}_{n}-\text{lat}_{n^{\prime}})^{2}}{r_{2,\text{lat}}^{2}}+\frac{(\text{lon}_{n}-\text{lon}_{n^{\prime}})^{2}}{r_{2,\text{lon}}^{2}}+\frac{(t_{n}-t_{n^{\prime}})^{2}}{r_{2,t}^{2}}\right]\right\\},$

where $\sigma_{1}^{2},r_{1,\text{lat}},r_{1,\text{lon}}$ and $r_{1,t}$
parameterize the mesoscale variation in currents whereas
$\sigma_{2}^{2},r_{2,\text{lat}},r_{2,\text{lon}}$ and $r_{2,t}$ parameterize
the submesoscale variation. As in Lodise et al. (2020), the latitudinal and
longitudinal components of $F$ are modeled as a priori independent. We choose
these parameters by maximal marginal likelihood (Rasmussen & Williams 2006,
Chapter 5) on an independent subset of the GLAD dataset. Estimates of the
underlying currents are obtained as the posterior mean of $F$ under this
model, which we take as the alternative, $\theta^{*}(y)$.

The second covariance function captures covariation among observations only at
the mesoscale. In this case, the Gaussian process prior has covariance
function

$k(\theta^{(i)}_{n},\theta^{(i)}_{n^{\prime}})=\sigma_{1}^{2}\exp\left\\{-\frac{1}{2}\left[\frac{(\text{lat}_{n}-\text{lat}_{n^{\prime}})^{2}}{r_{1,\text{lat}}^{2}}+\frac{(\text{lon}_{n}-\text{lon}_{n^{\prime}})^{2}}{r_{1,\text{lon}}^{2}}+\frac{(t_{n}-t_{n^{\prime}})^{2}}{r_{1,t}^{2}}\right]\right\\}+\sigma_{2}^{2}\mathbbm{1}[n=n^{\prime}],$

which maintains the same marginal variance but excludes submesoscale
covariances. We take the posterior mean under this model as the default
estimate $\hat{\theta}(y)$. See submesoscale_GP_c_value.ipynb for further
implementation details.