$\displaystyle\mathbb{E}\Big{[}\frac{I(A=1)I(G=O)}{p(A=1,G=O)}\\{\mathbb{E}[Y\mid
X,A=0,G=O]-\psi_{4}\\}S(A,G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{y,x}yp(y\mid x,A=0,G=O)\partial_{t}p_{t}(x\mid A=1,G=O)$
(29)
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=1)I(G=O)}{p(A=1,G=O)}\\{\mathbb{E}[Y\mid
X,A=0,G=O]-\psi_{4}\\}S(V)\Big{]}.$

Combining (28) and (29) concludes that

$\displaystyle\partial_{t}\psi_{4_{t}}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(A=1,G=O)}\Big{\\{}\frac{I(A=0)p(A=1\mid
X,G=O)}{p(A=0\mid X,G=O)}\\{Y-\mathbb{E}[Y\mid A,X,G]\\}$
$\displaystyle+I(A=1)\\{\mathbb{E}[Y\mid
X,A=0,G=O]-\psi_{4}\\}\Big{\\}}S(V)\Big{]}.$

Therefore

$\displaystyle\frac{I(G=O)}{p(A=1,G=O)}\Big{\\{}\frac{I(A=0)p(A=1\mid
X,G=O)}{1-p(A=1\mid X,G=O)}\\{Y-\mathbb{E}[Y\mid A,X,G]\\}$
$\displaystyle+I(A=1)\\{\mathbb{E}[Y\mid X,A=0,G=O]-\psi_{4}\\}\Big{\\}}$

is the influence function of $\psi_{4}$.

For $i\in\\{1,2,3,4\\}$, denote the obtained influence functions by
$IF_{\psi_{i}}$. The influence function for $\psi_{ETT}$ can be obtained as
$IF_{\psi_{ATE}}=\sum_{i=1}^{4}IF_{\psi_{i}}$. Therefore,

$\displaystyle\frac{1}{p(A=1,G=O)}\Big{\\{}$
$\displaystyle\frac{I(G=O)I(A=0)}{1-p(A=1\mid X,G=O)}\\{M-\mathbb{E}[M\mid
X,A,G]\\}$ $\displaystyle+\frac{I(G=E)I(A=0)}{1-p(A=1\mid
X,G=E)}\\{\frac{1}{p(G=E\mid X)}-1\\}\\{M-\mathbb{E}[M\mid X,A,G]\\}$
$\displaystyle+\frac{I(G=O)I(A=0)p(A=1\mid X,G=O)}{1-p(A=1\mid
X,G=O)}\\{Y-\mathbb{E}[Y\mid A,X,G]\\}$
$\displaystyle+I(G=O)\big{\\{}\mathbb{E}[M\mid X,A=0,G=O]+\mathbb{E}[M\mid
X,A=0,G=E]$ $\displaystyle\qquad\qquad\qquad+I(A=1)\\{Y+\mathbb{E}[Y\mid
X,A=0,G=O]-\psi_{ETT}\\}\big{\\}}\Big{\\}}$

is the influence function of $\psi_{ETT}$.

∎

###### Proof of Theorem 14.

For $a\in\\{0,1\\}$, define

$\displaystyle\psi_{1}^{(a)}=\mathbb{E}\big{[}\mathbb{E}[M\mid
A=a,X,G=E]\big{|}G=O\big{]},$
$\displaystyle\psi_{2}^{(a)}=\mathbb{E}\big{[}\mathbb{E}[M\mid
A=a,X,G=O]\big{|}G=O\big{]},$
$\displaystyle\psi_{3}^{(a)}=\mathbb{E}\big{[}\mathbb{E}[Y\mid
A=a,X,G=O]\big{|}G=O\big{]}.$

We use the notation $\partial_{t}f(t)$ to denote $\frac{\partial
f(t)}{\partial t}\big{|}_{t=0}$. For parameter $\psi^{(a)}$, let
$\psi^{(a)}_{t}$ be the parameter under a regular parametric sub-model indexed
by $t$, that includes the ground-truth model at $t=0$. Let $V$ be the set of
all observed variables. In order to obtain an influence function, we need to
find a random variable $\Gamma$ with mean zero, that satisfies

$\partial_{t}\psi^{(a)}_{t}=\mathbb{E}[\Gamma S(V)],$

where $S(V)=\partial_{t}\log p_{t}(V)$.

For $\psi_{1}^{(a)}$, note that

$\displaystyle\partial_{t}\psi^{(a)}_{1_{t}}$
$\displaystyle=\partial_{t}\sum_{m,x}mp_{t}(m\mid a,x,G=E)p_{t}(x\mid G=O)$
(30) $\displaystyle=\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle\quad+\sum_{m,x}mp(m\mid a,x,G=E)\partial_{t}p_{t}(x\mid G=O).$

For the first term in (30), we have

$\displaystyle\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle=\sum_{m,x}m\partial_{t}p_{t}(m\mid
a,x,G=E)\\{\frac{1}{p(G=E\mid x)}-1\\}\frac{1}{p(G=O)}p(x,G=E)$
$\displaystyle=\sum_{m,x}mS(m\mid a,x,G=E)\\{\frac{1}{p(G=E\mid
x)}-1\\}\frac{1}{p(A=a\mid x,G=E)}\cdot\frac{1}{p(G=O)}p(m,a,x,G=E)$
$\displaystyle=\sum_{m,\tilde{a},x,g}mS(m\mid\tilde{a},x,g)\\{\frac{1}{p(G=E\mid
x)}-1\\}\frac{I(\tilde{a}=a)}{p(A=a\mid
x,G=E)}\cdot\frac{I(g=E)}{p(G=O)}p(m,\tilde{a},x,g)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}M\\{\frac{1}{p(G=E\mid X)}-1\\}S(M\mid
A,X,G)\Big{]}$ $\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}S(M\mid A,X,G)\Big{]},$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}S(A,X,G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid G=O)$ (31)
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}S(V)\Big{]}.$

For the second term in (30), we have

$\displaystyle\sum_{m,x}mp(m\mid a,x,G=E)\partial_{t}p_{t}(x\mid G=O)$
$\displaystyle=\sum_{x}\mathbb{E}[M\mid A=a,X=x,G=E]\partial_{t}p_{t}(x\mid
G=O)$ $\displaystyle=\sum_{x}\mathbb{E}[M\mid A=a,X=x,G=E]S(x\mid G=O)p(x\mid
G=O)$ $\displaystyle=\sum_{x,g}\frac{I(g=O)}{p(G=O)}\mathbb{E}[M\mid
A=a,X=x,G=E]S(x\mid g)p(x,g)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\mathbb{E}[M\mid
A=a,X,G=E]S(X\mid G)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\mathbb{E}[\mathbb{E}[M\mid A=a,X,G=E]\mid G=O]\\}S(X\mid
G)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\psi_{1}^{(a)}\\}S(X\mid G)\Big{]}.$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\psi_{1}^{(a)}\\}S(G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}mp(m\mid a,x,G=E)\partial_{t}p_{t}(x\mid
G=O)=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\psi_{1}^{(a)}\\}S(V)\Big{]}.$ (32)

Combining (31) and (32) concludes that

$\displaystyle\partial_{t}\psi^{(a)}_{1_{t}}$
$\displaystyle=\mathbb{E}\Big{[}\Big{\\{}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}$
$\displaystyle\qquad+\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\psi_{1}^{(a)}\\}\Big{\\}}S(V)\Big{]}.$

Therefore,

$\displaystyle\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}$
$\displaystyle+\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=E]-\psi_{1}^{(a)}\\}$

is the influence function of $\psi_{1}^{(a)}$.

For $\psi_{2}^{(a)}$, note that

$\displaystyle\partial_{t}\psi^{(a)}_{2_{t}}$
$\displaystyle=\partial_{t}\sum_{m,x}mp_{t}(m\mid a,x,G=O)p_{t}(x\mid G=O)$
(33) $\displaystyle=\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=O)p(x\mid G=O)$
$\displaystyle\quad+\sum_{m,x}mp(m\mid a,x,G=O)\partial_{t}p_{t}(x\mid G=O).$

For the first term in (33), we have

$\displaystyle\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=O)p(x\mid G=O)$
$\displaystyle=\sum_{m,x}mS(m\mid a,x,G=O)\frac{1}{p(A=a\mid
x,G=O)}\cdot\frac{1}{p(G=O)}p(m,a,x,G=O)$
$\displaystyle=\sum_{m,\tilde{a},x,g}mS(m\mid\tilde{a},x,g)\frac{I(\tilde{a}=a)}{p(A=a\mid
x,G=O)}\cdot\frac{I(g=O)}{p(G=O)}p(m,\tilde{a},x,g)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\cdot\frac{I(G=O)}{p(G=O)}MS(M\mid A,X,G)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\cdot\frac{I(G=O)}{p(G=O)}\\{M-\mathbb{E}[M\mid A,X,G]\\}S(M\mid
A,X,G)\Big{]},$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\cdot\frac{I(G=O)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}S(A,X,G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}m\partial_{t}p_{t}(m\mid a,x,G=O)p(x\mid G=O)$ (34)
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\cdot\frac{I(G=O)}{p(G=O)}\\{M-\mathbb{E}[M\mid A,X,G]\\}S(V)\Big{]}.$

For the second term in (33), we have

$\displaystyle\sum_{m,x}mp(m\mid a,x,G=O)\partial_{t}p_{t}(x\mid G=O)$
$\displaystyle=\sum_{x}\mathbb{E}[M\mid A=a,X=x,G=O]\partial_{t}p_{t}(x\mid
G=O)$ $\displaystyle=\sum_{x}\mathbb{E}[M\mid A=a,X=x,G=O]S(x\mid G=O)p(x\mid
G=O)$ $\displaystyle=\sum_{x,g}\frac{I(g=O)}{p(G=O)}\mathbb{E}[M\mid
A=a,X=x,G=O]S(x\mid g)p(x,g)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\mathbb{E}[M\mid
A=a,X,G=O]S(X\mid G)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=O]-\mathbb{E}[\mathbb{E}[M\mid A=a,X,G=O]\mid G=O]\\}S(X\mid
G)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=O]-\psi_{2}^{(a)}\\}S(X\mid G)\Big{]}.$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=O]-\psi_{2}^{(a)}\\}S(G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}mp(m\mid a,x,G=O)\partial_{t}p_{t}(x\mid
G=O)=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=O]-\psi_{2}^{(a)}\\}S(V)\Big{]}.$ (35)

Combining (34) and (35) concludes that

$\displaystyle\partial_{t}\psi^{(a)}_{2_{t}}$
$\displaystyle=\mathbb{E}\Big{[}\Big{\\{}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\cdot\frac{I(G=O)}{p(G=O)}\\{M-\mathbb{E}[M\mid A,X,G]\\}$
$\displaystyle\qquad+\frac{I(G=O)}{p(G=O)}\\{\mathbb{E}[M\mid
A=a,X,G=O]-\psi_{2}^{(a)}\\}\Big{\\}}S(V)\Big{]}.$

Therefore,

$\displaystyle\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\\{M-\mathbb{E}[M\mid A,X,G]\\}+\mathbb{E}[M\mid
A=a,X,G=O]-\psi_{2}^{(a)}\big{\\}}$

is the influence function of $\psi_{2}^{(a)}$.

Similarly,

$\displaystyle\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{I(A=a)}{p(A=a\mid
X,G=O)}\\{Y-\mathbb{E}[Y\mid A,X,G]\\}+\mathbb{E}[Y\mid
A=a,X,G=O]-\psi_{3}^{(a)}\big{\\}}$

is the influence function of $\psi_{3}^{(a)}$.

For $i\in\\{1,2,3\\}$ and $a\in\\{0,1\\}$, denote the obtained influence
functions by $IF_{\psi_{i}^{(a)}}$. The influence function for $\psi_{ATE}$
can be obtained as
$IF_{\psi_{ATE}}=IF_{\psi_{1}^{(1)}}-IF_{\psi_{1}^{(0)}}+IF_{\psi_{2}^{(0)}}-IF_{\psi_{2}^{(1)}}+IF_{\psi_{3}^{(1)}}-IF_{\psi_{3}^{(0)}}$.
Therefore,

$\displaystyle\frac{(-1)^{1-A}}{p(A\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\mathbb{E}[M\mid
A,X,G]\\}\\{\frac{1}{p(G=E\mid X)}-1\\}$
$\displaystyle+\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{(-1)^{1-A}}{p(A\mid
X,G=O)}\\{Y-\mathbb{E}[Y\mid A,X,G]-M+\mathbb{E}[M\mid A,X,G]\\}$
$\displaystyle+\mathbb{E}[Y\mid X,A=1,G=O]-\mathbb{E}[Y\mid X,A=0,G=O]$
$\displaystyle+\mathbb{E}[M\mid X,A=1,G=E]-\mathbb{E}[M\mid X,A=0,G=E]$
$\displaystyle+\mathbb{E}[M\mid X,A=0,G=O]-\mathbb{E}[M\mid
X,A=1,G=O]-\psi_{ATE}\big{\\}}$

is the influence function of $\psi_{ATE}$.

∎

###### Proof of Proposition 1.

First, suppose the set $\\{\mu^{G=E}_{M},\mu^{G=O}_{M},\mu^{G=O}_{Y}\\}$ is
correctly specified. We have

$\displaystyle\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}\\{M-\hat{\mu}^{G=O}_{M}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\\{M-\hat{\mu}^{G=E}_{M}(0,X)\\}$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}\\{Y-\hat{\mu}^{G=O}_{Y}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}\\{\mathbb{E}[M\mid
X,A=0,G=O]-\hat{\mu}^{G=O}_{M}(0,X)\\}$ $\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\\{\mathbb{E}[M\mid X,A=0,G=E]-\hat{\mu}^{G=E}_{M}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}\\{\mathbb{E}[Y\mid
X,A=0,G=O]-\hat{\mu}^{G=O}_{Y}(0,X)\\}$ $\displaystyle\hskip
71.13188pt+I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\theta_{ETT}.$

$\displaystyle\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\hat{\mu}^{G=E}_{M}(A,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{Y-\hat{\mu}^{G=O}_{Y}(A,X)-M+\hat{\mu}^{G=O}_{M}(A,X)\\}$
$\displaystyle\quad+\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}\\{\mathbb{E}[M\mid
X,A,G=E]-\hat{\mu}^{G=E}_{M}(A,X)\\}\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{\mathbb{E}[Y\mid
X,A,G=O]-\hat{\mu}^{G=O}_{Y}(A,X)-\mathbb{E}[M\mid
X,A,G=O]+\hat{\mu}^{G=O}_{M}(A,X)\\}$
$\displaystyle\quad+\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\theta_{ATE}.$

Second, suppose the set $\\{\pi^{G=E},\pi^{G=O},p(G=E\mid x)\\}$ is correctly
specified. We have

$\displaystyle\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}\\{M-\hat{\mu}^{G=O}_{M}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\\{M-\hat{\mu}^{G=E}_{M}(0,X)\\}$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}\\{Y-\hat{\mu}^{G=O}_{Y}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}M$
$\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}M$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}Y$
$\displaystyle\hskip 71.13188pt+\frac{I(G=O)\mathbb{E}[I(A=0)\mid
X,G=O]}{1-\hat{\pi}^{G=O}(X)}\\{-\hat{\mu}^{G=O}_{M}(0,X)\\}$
$\displaystyle\hskip 71.13188pt+\frac{I(G=E)\mathbb{E}[I(A=0)\mid
X,G=E]}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\\{-\hat{\mu}^{G=E}_{M}(0,X)\\}$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)\mathbb{E}[I(A=0)\mid
X,G=O]\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}\\{-\hat{\mu}^{G=O}_{Y}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}M$
$\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}M$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}Y$
$\displaystyle\hskip 71.13188pt+I(G=O)\\{-\hat{\mu}^{G=O}_{M}(0,X)\\}$
$\displaystyle\hskip 71.13188pt+\mathbb{E}[I(G=E)\mid
X]\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}\\{-\hat{\mu}^{G=E}_{M}(0,X)\\}$
$\displaystyle\hskip
71.13188pt+I(G=O)\mathbb{E}[I(A=1)\\{-\hat{\mu}^{G=O}_{Y}(0,X)\\}\mid X,G=O]$
$\displaystyle\hskip
71.13188pt+I(G=O)\big{\\{}\hat{\mu}^{G=O}_{M}(0,X)+\hat{\mu}^{G=E}_{M}(0,X)+I(A=1)\\{Y+\hat{\mu}^{G=O}_{Y}(0,X)\\}\big{\\}}\Big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{1}{p(A=1,G=O)}\Big{\\{}\frac{I(G=O)I(A=0)}{1-\hat{\pi}^{G=O}(X)}M$
$\displaystyle\hskip
71.13188pt+\frac{I(G=E)I(A=0)}{1-\hat{\pi}^{G=E}(X)}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}M$ $\displaystyle\hskip
71.13188pt+\frac{I(G=O)I(A=0)\hat{\pi}^{G=O}(X)}{1-\hat{\pi}^{G=O}(X)}Y$
$\displaystyle\hskip 71.13188pt+I(G=O)I(A=1)Y\Big{\\}}\Big{]}$
$\displaystyle=\theta_{ETT}.$

$\displaystyle\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}\\{M-\hat{\mu}^{G=E}_{M}(A,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{Y-\hat{\mu}^{G=O}_{Y}(A,X)-M+\hat{\mu}^{G=O}_{M}(A,X)\\}$
$\displaystyle\quad+\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}M\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\cdot\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{Y-M\\}$
$\displaystyle\quad+\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}\\{-\hat{\mu}^{G=E}_{M}(A,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\big{\\{}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{-\hat{\mu}^{G=O}_{Y}(A,X)+\hat{\mu}^{G=O}_{M}(A,X)\\}$
$\displaystyle\quad+\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}M\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\cdot\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{Y-M\\}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\\{-\hat{\mu}^{G=E}_{M}(1,X)+\hat{\mu}^{G=E}_{M}(0,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\big{\\{}\\{-\hat{\mu}^{G=O}_{Y}(1,X)+\hat{\mu}^{G=O}_{M}(1,X)+\hat{\mu}^{G=O}_{Y}(0,X)-\hat{\mu}^{G=O}_{M}(0,X)\\}$
$\displaystyle\quad+\hat{\mu}^{G=O}_{Y}(1,X)-\hat{\mu}^{G=O}_{Y}(0,X)+\hat{\mu}^{G=E}_{M}(1,X)-\hat{\mu}^{G=E}_{M}(0,X)+\hat{\mu}^{G=O}_{M}(0,X)-\hat{\mu}^{G=O}_{M}(1,X)\big{\\}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=E}(X)}\cdot\frac{I(G=E)}{p(G=O)}M\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\cdot\frac{(-1)^{1-A}}{1-A+(-1)^{1-A}\hat{\pi}^{G=O}(X)}\\{Y-M\\}\Big{]}$
$\displaystyle=\theta_{ATE}.$

Parts 3 and 4 can be proven by combining the techniques used in parts 1 and 2,
and thus we omit here.

∎

###### Proof of Theorem 15.

Recall the definition that for $a\in\\{0,1\\}$

$\psi^{(a)}=\mathbb{E}\big{[}\mathbb{E}[h(M,A,X)\mid
A=a,X,G=E]\big{|}G=O\big{]}.$

We use the notation $\partial_{t}f(t)$ to denote $\frac{\partial
f(t)}{\partial t}\big{|}_{t=0}$. For parameter $\psi^{(a)}$, let
$\psi^{(a)}_{t}$ be the parameter under a regular parametric sub-model indexed
by $t$, that includes the ground-truth model at $t=0$. Let $V$ be the set of
all observed variables. In order to obtain an influence function, we need to
find a random variable $\Gamma$ with mean zero, that satisfies

$\partial_{t}\psi^{(a)}_{t}=\mathbb{E}[\Gamma S(V)],$

where $S(V)=\partial_{t}\log p_{t}(V)$.

Note that

$\displaystyle\partial_{t}\psi^{(a)}_{t}$
$\displaystyle=\partial_{t}\sum_{m,x}h_{t}(m,a,x)p_{t}(m\mid
a,x,G=E)p_{t}(x\mid G=O)$ (36)
$\displaystyle=\sum_{m,x}\partial_{t}h_{t}(m,a,x)p(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle\quad+\sum_{m,x}h(m,a,x)\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid
G=O)$ $\displaystyle\quad+\sum_{m,x}h(m,a,x)p(m\mid
a,x,G=E)\partial_{t}p_{t}(x\mid G=O).$

For the first term in (36), we have

$\displaystyle\sum_{m,x}\partial_{t}h_{t}(m,a,x)p(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle=\sum_{m,x}\partial_{t}h_{t}(m,a,x)\frac{p(m\mid
a,x,G=E)}{p(m\mid a,x,G=O)p(a\mid x,G=O)}p(m,a,x\mid G=O)$
$\displaystyle=\sum_{z,m,x}\partial_{t}h_{t}(m,a,x)q(z,a,x)p(z,m,a,x\mid G=O)$
$\displaystyle=\sum_{z,m,\tilde{a},x}I(\tilde{a}=a)\partial_{t}h_{t}(m,\tilde{a},x)q(z,\tilde{a},x)p(z,m,\tilde{a},x\mid
G=O)$ $\displaystyle=\mathbb{E}[I(A=a)\partial_{t}h_{t}(M,A,X)q(Z,A,X)\mid
G=O]$
$\displaystyle=\mathbb{E}\big{[}I(A=a)\partial_{t}\mathbb{E}[h_{t}(M,A,X)\mid
Z,A,X,G=O]q(Z,A,X)\big{|}G=O\big{]}.$

Note that by Assumption 16 $(ii)$

$\displaystyle\mathbb{E}[Y-h(M,A,X)\mid Z,A,X,G=O]=0$
$\displaystyle\Rightarrow\partial_{t}\mathbb{E}_{t}[Y-h_{t}(M,A,X)\mid
Z,A,X,G=O]=0$
$\displaystyle\Rightarrow\mathbb{E}[\partial_{t}\\{Y-h_{t}(M,A,X)\\}\mid
Z,A,X,G=O]$ $\displaystyle\quad+\mathbb{E}[\\{Y-h(M,A,X)\\}S(Y,M\mid
Z,A,X,G=O)\mid Z,A,X,G=O]=0$
$\displaystyle\Rightarrow\mathbb{E}[\partial_{t}h_{t}(M,A,X)\mid
Z,A,X,G=O]=\mathbb{E}[\\{Y-h(M,A,X)\\}S(Y,M\mid Z,A,X,G=O)\mid Z,A,X,G=O].$

Therefore,

$\displaystyle\sum_{m,x}\partial_{t}h_{t}(m,a,x)p(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle=\mathbb{E}\big{[}I(A=a)q(Z,A,X)\mathbb{E}[\\{Y-h(M,A,X)\\}S(Y,M\mid
Z,A,X,G=O)\mid Z,A,X,G=O]\big{|}G=O\big{]}$
$\displaystyle=\mathbb{E}\big{[}I(A=a)q(Z,A,X)\\{Y-h(M,A,X)\\}S(Y,M\mid
Z,A,X,G=O)\big{|}G=O\big{]}$
$\displaystyle=\mathbb{E}\big{[}I(A=a)\frac{I(G=O)}{p(G=O)}q(Z,A,X)\\{Y-h(M,A,X)\\}S(Y,M\mid
Z,A,X,G)\big{]}.$

Also, note that

$\displaystyle\mathbb{E}\big{[}I(A=a)\frac{I(G=O)}{p(G=O)}q(Z,A,X)\\{Y-h(M,A,X)\\}S(Z,A,X,G)\big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}\partial_{t}h_{t}(m,a,x)p(m\mid a,x,G=E)p(x\mid
G=O)=\mathbb{E}\big{[}I(A=a)\frac{I(G=O)}{p(G=O)}q(Z,A,X)\\{Y-h(M,A,X)\\}S(V)\big{]}.$
(37)

For the second term in (36), we have

$\displaystyle\sum_{m,x}h(m,a,x)\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid G=O)$
$\displaystyle=\sum_{m,x}h(m,a,x)\partial_{t}p_{t}(m\mid
a,x,G=E)\\{\frac{1}{p(G=E\mid x)}-1\\}\frac{1}{p(G=O)}p(x,G=E)$
$\displaystyle=\sum_{m,x}h(m,a,x)S(m\mid a,x,G=E)\\{\frac{1}{p(G=E\mid
x)}-1\\}\frac{1}{p(A=a\mid x,G=E)}\cdot\frac{1}{p(G=O)}p(m,a,x,G=E)$
$\displaystyle=\sum_{m,\tilde{a},x,g}h(m,\tilde{a},x)S(m\mid\tilde{a},x,g)\\{\frac{1}{p(G=E\mid
x)}-1\\}\frac{I(\tilde{a}=a)}{p(A=a\mid
x,G=E)}\cdot\frac{I(g=E)}{p(G=O)}p(m,\tilde{a},x,g)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}h(M,A,X)\\{\frac{1}{p(G=E\mid X)}-1\\}S(M\mid
A,X,G)\Big{]}$ $\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{h(M,A,X)-\eta(A,X)\\}\\{\frac{1}{p(G=E\mid
X)}-1\\}S(M\mid A,X,G)\Big{]},$

where

$\displaystyle\eta(a,x)\coloneqq$ $\displaystyle\mathbb{E}[h(M,A,X)\mid
A=a,X=x,G=E]$ $\displaystyle=$ $\displaystyle\sum_{m}h(m,a,x)p(m\mid a,x,G=E)$
$\displaystyle=$
$\displaystyle\sum_{z,m,\tilde{a}}I(\tilde{a}=a)h(m,\tilde{a},x)q(z,\tilde{a},x)p(z,m,\tilde{a}\mid
x,G=E)$ $\displaystyle=$ $\displaystyle\mathbb{E}[I(A=a)h(M,A,X)q(Z,A,X)\mid
X=x,G=E].$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{h(M,A,X)-\eta(A,X)\\}\\{\frac{1}{p(G=E\mid
X)}-1\\}S(A,X,G)\Big{]}=0.$

Therefore,

$\displaystyle\sum_{m,x}h(m,a,x)\partial_{t}p_{t}(m\mid a,x,G=E)p(x\mid G=O)$
(38) $\displaystyle=\mathbb{E}\Big{[}\frac{I(A=a)}{p(A=a\mid
X,G=E)}\cdot\frac{I(G=E)}{p(G=O)}\\{h(M,A,X)-\eta(A,X)\\}\\{\frac{1}{p(G=E\mid
X)}-1\\}S(V)\Big{]}.$

For the third term in (36), we have

$\displaystyle\sum_{m,x}h(m,a,x)p(m\mid a,x,G=E)\partial_{t}p_{t}(x\mid G=O)$
$\displaystyle=\sum_{x}\sum_{m}h(m,a,x)p(m\mid a,x,G=E)S(x\mid G=O)p(x\mid
G=O)$ $\displaystyle=\sum_{x,g}\frac{I(g=O)}{p(G=O)}\eta(a,x)S(x\mid g)p(x,g)$
$\displaystyle=\mathbb{E}[\frac{I(G=O)}{p(G=O)}\eta(a,X)S(X\mid G)]$
$\displaystyle=\mathbb{E}[\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\mathbb{E}[\eta(a,X)\mid
G=O]\\}S(X\mid G)]$
$\displaystyle=\mathbb{E}[\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\psi^{(a)}\\}S(X\mid
G)].$

Note that

$\displaystyle\mathbb{E}[\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\psi^{(a)}\\}S(G)]=0.$

Therefore,

$\displaystyle\sum_{m,x}h(m,a,x)p(m\mid a,x,G=E)\partial_{t}p_{t}(x\mid
G=O)=\mathbb{E}[\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\psi^{(a)}\\}S(V)].$ (39)

Combining (36)-(39) concludes that

$\displaystyle\partial_{t}\psi^{(a)}_{t}$
$\displaystyle=\mathbb{E}\Big{[}\Big{\\{}\frac{I(G=O)}{p(G=O)}I(A=a)q(Z,A,X)\\{Y-h(M,A,X)\\}$
$\displaystyle\qquad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{p(A=a\mid
X,G=E)}\\{h(M,A,X)-\eta(A,X)\\}\\{\frac{1}{p(G=E\mid X)}-1\\}$
$\displaystyle\qquad+\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\psi^{(a)}\\}\Big{\\}}S(V)\Big{]}.$

Therefore,

$\displaystyle\frac{I(G=O)}{p(G=O)}I(A=a)q(Z,A,X)\\{Y-h(M,A,X)\\}$
$\displaystyle+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{p(A=a\mid
X,G=E)}\\{h(M,A,X)-\eta(A,X)\\}\\{\frac{1}{p(G=E\mid X)}-1\\}$
$\displaystyle+\frac{I(G=O)}{p(G=O)}\\{\eta(a,X)-\psi^{(a)}\\}$

is the influence function of $\psi^{(a)}$.

∎

###### Proof of Proposition 2.

First, suppose the pair $\\{h,p(m\mid a,x,G=E)\\}$ is correctly specified. We
have

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\\{Y-\hat{h}(M,a,X)\\}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\underbrace{\mathbb{E}[Y-\hat{h}(M,a,X)\mid
Z,A=a,G=O]}_{=0}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\underbrace{\mathbb{E}[\hat{h}(M,a,X)-\hat{\eta}(a,X)\mid
A=a,X,G=E]}_{=0}\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}$
$\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\theta^{(a)}.$

Second, suppose the set $\\{h,p(A=1\mid x,G=E),p(G=E\mid x)\\}$ is correctly
specified. We have

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\\{Y-\hat{h}(M,a,X)\\}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\underbrace{\mathbb{E}[Y-\hat{h}(M,a,X)\mid
Z,A=a,G=O]}_{=0}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\hat{h}(M,a,X)\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}$
$\displaystyle\quad+\underbrace{\mathbb{E}\big{[}\frac{I(G=O)}{p(G=O)}-\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\cdot\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\big{|}X\big{]}}_{=0}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\hat{h}(M,a,X)\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}\Big{]}$
$\displaystyle=\theta^{(a)}.$

Third, suppose the set $\\{q,p(A=1\mid x,G=E),p(G=E\mid x)\\}$ is correctly
specified. We have

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\\{Y-\hat{h}(M,a,X)\\}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)Y$
$\displaystyle\quad+\Big{\\{}\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\cdot\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}-\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\Big{\\}}\hat{h}(M,a,X)$
$\displaystyle\quad+\underbrace{\mathbb{E}\big{[}\frac{I(G=O)}{p(G=O)}-\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\cdot\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}\big{|}X\big{]}}_{=0}\hat{\eta}(a,X)\Big{]}.$

Note that

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\hat{h}(M,a,X)\Big{]}$
$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\mathbb{E}[\hat{q}(Z,a,X)\mid
M,A=a,X,G=O]\hat{h}(M,a,X)\Big{]}$
$\displaystyle=\frac{1}{p(G=O)}\sum_{m,x}\frac{p(m\mid a,x,G=E)}{p(m\mid
a,x,G=O)p(a\mid x,G=O)}\hat{h}(m,a,x)p(m,a,x,G=O)$
$\displaystyle=\frac{1}{p(G=O)}\sum_{m,x}p(m\mid
a,x,G=E)p(x,G=O)\hat{h}(m,a,x)$
$\displaystyle=\frac{1}{p(G=O)}\sum_{m,x}\frac{1}{p(a,x,G=E)}p(x,G=O)\hat{h}(m,a,x)p(m,a,x,G=E)$
$\displaystyle=\frac{1}{p(G=O)}\sum_{m,x}\frac{1}{p(a\mid
x,G=E)}\frac{p(G=O\mid x)}{p(G=E\mid x)}\hat{h}(m,a,x)p(m,a,x,G=E)$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\cdot\\{\frac{1}{\hat{p}(G=E\mid X)}-1\\}\hat{h}(M,a,X)\Big{]}.$

Therefore,

$\displaystyle\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)\\{Y-\hat{h}(M,a,X)\\}$
$\displaystyle\quad+\frac{I(G=E)}{p(G=O)}\cdot\frac{I(A=a)}{1-a+(-1)^{1-a}\hat{p}(A=1\mid
X,G=E)}\\{\hat{h}(M,a,X)-\hat{\eta}(a,X)\\}\\{\frac{1}{\hat{p}(G=E\mid
X)}-1\\}$ $\displaystyle\quad+\frac{I(G=O)}{p(G=O)}\hat{\eta}(a,X)\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O)}I(A=a)\hat{q}(Z,a,X)Y\Big{]}$
$\displaystyle=\theta^{(a)}.$

∎

###### Proof of Proposition 3.

$\displaystyle\mathbb{E}[q(Z,A,X)\mid M,A,X,G=O]=\frac{p(M\mid
A,X,G=E)}{p(M\mid A,X,G=O)p(A\mid X,G=O)}$
$\displaystyle\Rightarrow\sum_{z}q(z,a,x)p(z,m,a\mid x,G=O)=p(m\mid a,x,G=E)$
$\displaystyle\Rightarrow\sum_{z}q(z,a,x)\frac{p(z,G=O\mid m,a,x)p(m,a\mid
x)}{p(G=O\mid x)}=\frac{p(G=E\mid m,a,x)p(m,a\mid x)}{p(a,G=E\mid x)}$
$\displaystyle\Rightarrow\sum_{z,g}I(g=O)q(z,a,x)\frac{p(z,g\mid
m,a,x)}{p(G=O\mid x)}=\sum_{g}I(g=E)\frac{p(g\mid m,a,x)}{p(a,G=E\mid x)}$
$\displaystyle\Rightarrow\mathbb{E}\Big{[}\frac{I(G=O)}{p(G=O\mid
X)}q(Z,A,X)-\frac{I(G=E)}{p(A,G=E\mid X)}\Big{|}M,A,X\Big{]}=0.$

∎


# Computational methods for fast Bayesian model assessment via calibrated
posterior p-values

Sally Paganin
Department of Environmental Science, Policy, and Management,
University of California, Berkeley
Department of Biostatistics, Harvard T.H. Chan School of Public Health
and
Perry de Valpine
Department of Environmental Science, Policy, and Management,
University of California, Berkeley  The authors gratefully acknowledge please
remember to list all relevant funding sources in the unblinded version

###### Abstract

Posterior predictive p-values (ppps) have become popular tools for Bayesian
model criticism, being general-purpose and easy to use. However, their
interpretation can be difficult because their distribution is not uniform
under the hypothesis that the model did generate the data. To address this
issue, procedures to obtain calibrated ppps (cppps) have been proposed
although not used in practice, because they require repeated simulation of new
data and model estimation via MCMC. Here we give methods to balance the
computational trade-off between the number of calibration replicates and the
number of MCMC samples per replicate. Our results suggest that investing in a
large number of calibration replicates while using short MCMC chains can save
significant computation time compared to naive implementations, without
significant loss in accuracy. We propose different estimators for the variance
of the cppp that can be used to confirm quickly when the model fits the data
well. Variance estimation requires the effective sample sizes of many short
MCMC chains; we show that these can be well approximated using the single long
MCMC chain from the real-data model. The procedure for cppp is implemented in
NIMBLE, a flexible framework for hierarchical modeling that supports many
models and discrepancy measures.

Keywords: Bayesian Goodness-of-fit; Double simulation; Hierarchical models;
Model criticism; NIMBLE

##  1 Introduction

An important step in statistical modeling is to ask whether a model fits the
data reasonably. In applied Bayesian statistics using Markov chain Monte Carlo
(MCMC) algorithms, a common method has been the use of posterior predictive
checks, which compare some function of the data to a reference distribution,
such as the posterior predictive distribution (Rubin, 1984; Meng, 1994; Gelman
et al., 1996). The result of this comparison is typically summarized by a
“p-value”, which serves as an assessment of whether, for some chosen measure
of disagreement between model and data, the data are so unusual as to doubt
the veracity of the model. These posterior predictive p-values (hereafter,
ppps) are considered a measure of “surprise”, with low values indicating that
the data are incompatible with the model.

Unfortunately, it is well known that posterior predictive checks are liberal
in the sense of concluding that models are better than they really are.
Correspondingly, posterior predictive p-values lack the frequentist property
of following a uniform distribution if the model did generate the data,
clustering instead around a value of 0.5 (Gelman et al., 2013b). This makes
ppps difficult to interpret and leaves this common practice of Bayesian model
assessment on shaky ground. The ppp is a measure of “surprise” in units that
sound like probability, in the same sense as a frequentist p-value, but are
not.

For this reason, many authors have discussed the need for calibration of ppps
to set an interpretable scale (Robins et al., 2000; Hjort et al., 2006;
Steinbakk and Storvik, 2009). In particular, Hjort et al. (2006) propose
calibration of ppps via a bootstrap-like procedure: simulate many data sets
from the reference distribution and go through the steps of MCMC and
calculation of ppp to obtain a simulated null distribution of ppps. From this
null distribution, the frequentist p-value of the Bayesian ppp can be
determined, providing a more standardized (or objective) measure of model
plausibility. Since the goal is frequentist calibration of a Bayesian method,
we speak informally of a null distribution and clarify this notion in Section
2.

Calibrated ppps (cppps) are appealing because they can provide an automated
method for a large class of problems with few additional assumptions. However,
the major difficulty with cppps is the potentially prohibitive computational
burden of obtaining them. It is not uncommon for analysts to run MCMC for
hours or days, so repeating this step a large number of times can be daunting
in practice. There are some approaches (Johnson, 2007; Yuan and Johnson, 2012)
that allow performing goodness-of-fit tests without the need for a calibration
procedure. These approaches use measures of discrepancy based on pivotal
quantities, i.e. statistics of data, parameters, or both, having a known
probability distribution that is independent of unknown parameters. However,
these methods are limited in that they are restricted to particular choices of
models and discrepancies.

In this paper, we give computational methods for approximating the cppp to a
useful degree of accuracy much faster than would be achieved by a naive
implementation. We seek a model-generic method, one which makes no assumptions
of analytic tractability and can work for any discrepancy measure of interest.
The core computational trade-off is between the number of calibration
replicates (bootstrap samples) and the number of MCMC samples of each
replicate for an accurate approximation of its ppp. Guided by theoretical
examination of bias and variance in the cppp estimator, we show that only a
small MCMC sample for each calibration replicate is necessary. MCMC samples
with effective sample size (ESS)—of the indicator used in ppp
calculation—between $50$ and $200$ can be enough to make bias negligible.
Beyond that, it is optimal to spend computational effort on more calibration
replicates to reduce the variance of the cppp estimate.

Critically, in the case of a good model fit, even a fairly rough (fast)
estimate of cppp may be adequate to conclude the model is acceptable. For
example, a report that the cppp is 0.4 with 95% Monte Carlo confidence
interval from 0.3-0.5 would be rough but adequate to conclude the model is
acceptable. Therefore, it is important to estimate the Monte Carlo standard
error of the cppp estimate itself. However, this requires an estimate of the
ESS of the MCMC sample for each calibration replicate. And that creates a
conundrum because the ESS estimates from short MCMC chains will themselves be
inaccurate, hindering precise estimation of the cppp uncertainty. To address
this, we propose a “transfer” approach that uses the properties of the MCMC
chain from the real data as a basis for approximating ESS in the many MCMC
samples of the calibration replicate. The chain for the real data will be much
longer and hence more informative about mixing properties. As this step is
somewhat heuristic, one can build in a conservative buffer. This makes it
possible to estimate the cppp and its Monte Carlo uncertainty from short MCMC
runs for many calibration replicates.

Together, these two advances – using short MCMC chains for each calibration
replicate and using few calibration replicates when a rough estimate of the
cppp is sufficient – can save one or more orders of magnitude of computational
effort compared to a naive implementation to compute the cppp. By a “naive
implementation”, we mean one where the analyst runs an MCMC for each
calibration replicate that is as long as the MCMC for the real data and uses
an unnecessarily large number of calibration replicates.

The rest of the paper is organized as follows. In Section 2, we set up the ppp
and cppp problems. In Section 3, we give the ppp and cppp estimators and
develop expressions for the bias and variance of the latter. In Section 4, we
study the bias-variance trade-off of allocating computational effort between
the number of calibration replicates and the number of MCMC samples per
calibration replicate. The theory from Section 3 suggests that the entire cppp
estimation process can be understood as similar to a beta-binomial estimation
problem, which we use to give theoretical results on the bias-variance trade-
off. Section 5 discusses the estimation of the standard error of cppp
estimates, including plug-in and bootstrap estimates that use the heuristic
proposal to transfer mixing properties of the long MCMC chain to the many
short calibration MCMC chains. Section 6 gives two worked examples, one from
physics and one from ecology, with scenarios of models that should be accepted
and rejected. The procedure for cppp is implemented in the NIMBLE software (de
Valpine et al., 2022; de Valpine et al., 2017) a flexible R-based system for
hierarchical models, and code to reproduce examples in the paper is available
at https://github.com/salleuska/fastCPPP.

##  2 Model-checking, ppp and cppp

Suppose data $y$ are modeled as having been generated from $p(y|\theta)$, with
$\theta$ including unknown parameters and latent states. Typically some prior
distribution $\pi(\theta)$ is assumed, and $\theta$ is estimated in a Bayesian
framework by its posterior $p(\theta|y)$. Bayesian model checking focuses on
assessing whether or not the assumed model, defined as the combination of data
distribution $p(y|\theta)$ and prior $\pi(\theta)$, is compatible with data.
Informally, we look for evidence against the null hypothesis that the posited
model did really generate the data. We refer to Robins et al. (2000) for a
rigorous definition of the testing problem in this context, and analysis of
asymptotic properties of Bayesian p-values. For model-checking, a discrepancy
measure $D(y,\theta)$ is chosen, reflecting aspects of the data that it is
important that the model describes well. The term _discrepancy_ has been
introduced by Gelman et al. (1996), where $D(y,\theta)$ can be thought of as a
generalization of a test statistic whose distribution can depend on $\theta$.

Define $f(y^{*},\theta|y)=p(y^{*}|\theta)p(\theta|y)$ as the joint posterior
predictive distribution of $\theta$ and hypothetical data $y^{*}$ given
observed data $y$. The posterior predictive distribution of $y^{*}$ is $\int
f(y^{*},\theta|y)d\theta=f(y^{*}|y)$. The ppp is the probability that the
discrepancy measure calculated for a draw of the posterior predictive
distribution is more extreme than the discrepancy measure of the data,
averaged over the posterior $p(\theta|y)$. This is given by

$\mbox{ppp}(y)=\text{Pr}\left\\{D(Y^{*},\Theta)\geq D(y,\Theta)|y\right\\}.$
(1)

Letting $\Delta(y^{*},\theta|y)=D(y^{*},\theta)-D(y,\theta)$, we rewrite (1)
as

$\mbox{ppp}(y)=\mathbb{E}_{(Y^{*},\Theta)}\left[\mathbb{I}\left\\{\Delta(Y^{*},\Theta|y)\geq
0\right\\}\right],$ (2)

where the expectation is over $(Y^{*},\Theta)\sim f(y^{*},\theta|y)$. We
interpret the ppp as a test statistic, or simply a statistic, rather than a
“p-value”, because it does not behave like a p-value.

The cppp is designed to be the proper p-value for the ppp statistic. That is,
the cppp is defined as the probability, if the model is valid, of obtaining a
ppp at least as extreme (small) as $\mbox{ppp}(y)$. In this case, the
probability is over the space of hypothetical data, $\widetilde{Y}\sim
g(\tilde{y}|y)$, for some relevant reference density $g(\tilde{y}|y)$. The
cppp is given by

$\mbox{cppp}(y)=\text{Pr}\\{\mbox{ppp}(\widetilde{Y})\leq\mbox{ppp}(y)|y\\}=\mathbb{E}_{\widetilde{Y}}\left[\mathbb{I}\left\\{\mbox{ppp}(\widetilde{Y})\leq\mbox{ppp}(y)\right\\}\right].$
(3)

The distribution of $\mbox{ppp}(\widetilde{Y})$, with $\widetilde{Y}\sim
g(\cdot)$, is used here as a null distribution for the test statistic
$\mbox{ppp}(y)$. The choice of the reference density $g(\cdot)$ is important.

The dependence of $g(\cdot)$ on the data $y$ means that even the cppp may not
be perfectly calibrated. If one knew the true parameters, $\theta_{0}$, then
the choice $g(\tilde{y}|y)=g(\tilde{y})=p(\tilde{y}|\theta_{0})$ would
guarantee that $\mbox{cppp}(y)$ is perfectly calibrated. For example, it would
give an accurate Type I error rate, $\alpha$, if one chooses as a decision
rule to reject the model if $\mbox{cppp}(y)<\alpha$. Since $\theta_{0}$ is
unknown, in the examples below we choose the posterior predictive
distribution, $g(\tilde{y}|y)=f(\tilde{y},\theta|y)$. As long as
$\text{Pr}\\{\mbox{ppp}(\widetilde{Y})\leq\mbox{ppp}(y)|\theta\\}$ does not
vary strongly over values of $\theta$ with high posterior density, calibration
based on this choice of $g(\cdot)$ should be fairly accurate. This is similar
to the issue that a parametric bootstrap may not be perfectly accurate because
it depends on estimated parameters rather than true parameters.

Notice that different choices of the reference density $g(\cdot)$ lead to
alternative definitions of Bayesian p-values. For example, Hjort et al. (2006)
focus on the prior predictive distribution (Box, 1980), which does not depend
on data $y$, with comments on other choices such as the posterior predictive.
Other important variations and alternative Bayesian p-values are introduced in
Bayarri and Berger (2000) and further discussed by Robins et al. (2000) and
Bayarri and Castellanos (2007).

Our choice is motivated by applied problems where the prior predictive
distribution based on uninformative priors for $\theta$ would give a
distribution of $\mbox{ppp}(\widetilde{Y})$ very different from its
distribution based on $\theta$ close to some $\theta_{0}$. That is, the
distribution of $\mbox{ppp}(\widetilde{Y})$, resulting from $\widetilde{Y}\sim
p(\tilde{y}|\theta)$, may vary substantially over very large ranges of
$\theta$ such as the range supported by an uninformative prior. The
computational methods presented here can be used for any choice of $g(\cdot)$,
so we do not delve further into issues around the choice of $g(\cdot)$.

Note that the posterior predictive distribution has two roles. It is the
distribution of $(Y^{*},\Theta)$ in (2), and it is a sensible choice for
$g(\tilde{y}|y)$ in (3). These roles are distinct, as different choices of
$g(\cdot)$ are possible.

##  3 Monte Carlo estimates of ppp and cppp

Figure 1: Schematic representation of the Monte Carlo approximation for the
ppp.

In practice, we work with Monte Carlo approximations of all expected values in
(2)-(3). Figure 1 provides a summary of the quantities involved in the
estimation. Define a posterior sample of size $m$ based on the real data as
$S_{1}(y,m)=\left\\{\theta^{(i)}\right\\},i=1,\ldots,m$, resulting from an
MCMC or other algorithm such that $\theta^{(i)}\sim p(\theta|y)$ for
$i=1,\ldots,m$. The $\theta^{(i)}$s will typically be sequentially correlated
because they are sampled by MCMC. Define the collection of samples from the
posterior predictive distribution based on the real data as
$S_{2}(y,m)=\left\\{(y^{*(i)},\theta^{(i)})\right\\},i=1,\ldots,m$, with
$\theta^{(i)}$ from $S_{1}(y,m)$ and $y^{*(i)}\sim p(y^{*}|\theta^{(i)})$.
Define $k(S_{2}(y,m))=\sum_{i=1}^{m}\mathbb{I}^{(i)}$, where
$\mathbb{I}^{(i)}=\mathbb{I}\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)\geq
0\right\\}$, so that $k$ is the count of posterior predictive discrepancies
$D(y^{*(i)},\theta^{(i)})$ that are more extreme than data discrepancies
$D(y,\theta^{(i)})$. Then the standard Monte Carlo estimate of $\mbox{ppp}(y)$
is

$\widehat{\mbox{ppp}}(y)=\frac{k}{m}.$ (4)

Notice that each $\mathbb{I}^{(i)}$ draw is marginally distributed as
Bernoulli with probability $\mbox{ppp}(y)$. Hence, $k$ is a realization of a
random variable parameterized by $y$ and $m$, and we can write $K\sim
h(k|y,m)$. If the posterior draws are independent, the distribution of $K$ is
binomial, but in the general case, the $\mathbb{I}^{(i)}$s are serially
dependent. Using standard theory for MCMC output (e.g. Robert and Casella,
1999, Chapter 12), define the integrated autocorrelation time of the chain of
indicator variables $\left\\{\mathbb{I}^{(i)}\right\\},i=1,\ldots,m$, as
$\tau$, so that its effective sample size is $\mbox{ESS}=m/\tau$. Then we have
that the mean and variance of $K$ are

$\displaystyle\mathbb{E}[K|y]$ $\displaystyle=m\mbox{ppp}(y),$
$\displaystyle\mathbb{V}[K|y]$
$\displaystyle=\frac{m^{2}\mbox{ppp}(y)(1-\mbox{ppp}(y))}{\mbox{ESS}}.$

For some purposes, one may want to replace $\frac{k}{m}$ with
$\frac{k+0.5}{m+1}$, which is asymptotically equivalent but avoids zeros.

Figure 2: Schematic representation of the Monte Carlo approximation for the
cppp.

The foregoing will be used once for the real data and then repeatedly for
calibration replicates, as described next. For theoretical purposes, we assume
that the number of MCMC samples $m$ is sufficiently large for the real data so
that $\widehat{\mbox{ppp}}(y)$ is very accurate, and hence assume that
$\mbox{ppp}(y)$ computed on observed data is essentially known, or at least
has much smaller error than that of the calibration replicates.

To estimate $\mbox{cppp}(y)$, we need a second layer of Monte Carlo
approximation (see Figure 2). Define the collection of calibration replicates
as $S_{3}(y,r)=\left\\{\tilde{y}^{(j)}\right\\}$, $j=1\ldots r$ with
$\tilde{y}^{(j)}\sim g(\tilde{y}|y)$ and all $\tilde{y}^{(j)}$ mutually
independent. Given each $\tilde{y}^{(j)}$, new posterior samples
$S_{1}(\tilde{y}^{(j)},\tilde{m}^{(j)})$ and samples from the posterior
predictive $S_{2}(\tilde{y}^{(j)},\tilde{m}^{(j)})$ are drawn using a sample
size choice, $\tilde{m}^{(j)}$, possibly smaller from the $m$ used for the
real data. These new samples come from running the MCMC with the calibration
replicate data, $\tilde{y}^{(j)}$, in the role of observed data. We initialize
this MCMC with the same value of $\theta$ in $S_{1}(y,m)$ used to generate the
calibration replicate data $\tilde{y}^{(j)}$. Hence we do not need a burn-in
for the calibration replicates MCMC, since we are starting in the posterior
region of interest.

From these samples, the value of interest is
$\tilde{k}^{(j)}=k(S_{2}(\tilde{y}^{(j)},\tilde{m}^{(j)}))$. The Monte Carlo
estimate of $\mbox{cppp}(y)$ is then:

$\displaystyle\widehat{\mbox{cppp}}(y)=\frac{1}{r}\sum_{j=1}^{r}\mathbb{I}\left\\{\tilde{k}^{(j)}\leq\tilde{m}^{(j)}\mbox{ppp}(y)\right\\}=\frac{1}{r}\sum_{j=1}^{r}\mathbb{I}\left\\{\widehat{\mbox{ppp}}(\tilde{y}^{(j)})\leq\mbox{ppp}(y)\right\\}.$
(5)

It is important to notice that we are not interested in
$\tilde{k}^{(j)}/\tilde{m}^{(j)}$ as an estimate of
$\widehat{\mbox{ppp}}(\tilde{y}^{(j)})$. Rather, we care about it as an
estimate of whether $\widehat{\mbox{ppp}}(\tilde{y}^{(j)})\leq\mbox{ppp}(y)$.
From here on we assume for simplicity that $\tilde{m}^{(j)}=\tilde{m}$ for all
$j$s, i.e. all calibration replicates use the same sample size.

Since the reference density $g(\cdot)$ is chosen as the posterior predictive
distribution based on the real data, $S_{3}(y,r)$ will in practice come from a
subsample of $S_{2}(y,m)$. Thus, $S_{2}(y,m)$ will be used twice: once to
estimate the ppp of the real data, $\mbox{ppp}(y)$, and once to approximate
the distribution of $\mbox{ppp}(\widetilde{Y})$, where $\widetilde{Y}\sim
g(\tilde{y}|y)$, for the purpose of estimating $\mbox{cppp}(y)$. In the latter
role, we aim to use relatively few samples to reduce computation, $r\ll m$, so
the $S_{3}(y,r)$ samples are drawn randomly or by systematically thinning
$S_{2}(y,m)$, making it reasonable to assume the $S_{3}(y,r)$ samples are
mutually independent.

###  3.1 Bias and variance of $\widehat{\mbox{cppp}}(y)$

To obtain insights on the role of the number of calibration replicates $r$ and
MCMC samples $\tilde{m}$, we look at bias and variance of
$\widehat{\mbox{cppp}}(y)$ with respect to its theoretical definition.

The expected value of $\widehat{\mbox{cppp}}(y)$ in (5) is
$\mathbb{E}_{\widetilde{Y}}\left[\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y}\right]\right]$,
where $E_{\widetilde{Y}}$ is expectation over $\widetilde{Y}\sim
g(\tilde{y}|y)$ and $E_{\widetilde{K}|\widetilde{Y}}$ is expectation over
$\widetilde{K}\sim h(\tilde{k}|\tilde{y},m)$. The $\mbox{cppp}(y)$ definition
(3) can be written as
$\mathbb{E}_{\widetilde{Y}}\left[\mathbb{I}\left\\{\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\widetilde{K}|\widetilde{Y}\right]\leq\tilde{m}\mbox{ppp}(y)\right\\}\right]$.
The bias of $\widehat{\mbox{cppp}}(y)$ is thus

$\mathbb{E}_{\widetilde{Y}}\left[\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y}\right]-\mathbb{I}\left\\{\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\widetilde{K}|\widetilde{Y}\right]\leq\tilde{m}\mbox{ppp}(y)\right\\}\right].$
(6)

To gain some intuition about the bias, notice that the expression inside the
outer expectation $\mathbb{E}_{\widetilde{Y}}[\cdot]$ is the average signed
error of determining whether $\mbox{ppp}(\tilde{y})\leq\mbox{ppp}(y)$ using a
draw $\widetilde{K}$, i.e. using the samples $S_{2}(\tilde{y},\tilde{m})$
(which would also give the estimate $\widehat{\mbox{ppp}}(\tilde{y})$). When
$\mbox{ppp}(\tilde{y})\leq\mbox{ppp}(y)$, the second indicator is 1, and
$\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y}\right]$
is the probability that $\widehat{\mbox{ppp}}(\tilde{y})\leq\mbox{ppp}(y)$
from a draw of $\widetilde{K}$. When $\mbox{ppp}(\tilde{y})>\mbox{ppp}(y)$,
the second indicator is 0, and
$\mathbb{E}_{\widetilde{K}|\widetilde{Y}}\left[\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y}\right]$
is one minus the probability that
$\widehat{\mbox{ppp}}(\tilde{y})>\mbox{ppp}(y)$. As $m\rightarrow\infty$,
$\widehat{\mbox{ppp}}(\tilde{y})\rightarrow\mbox{ppp}(\tilde{y})$, and the
bias goes to zero.

The variance of $\widehat{\mbox{cppp}}(y)$ can be obtained using the law of
total variance

$\displaystyle\mathbb{V}[\widehat{\mbox{cppp}}(y)]=$
$\displaystyle\frac{1}{r}\left\\{\mathbb{E}_{\widetilde{Y}}\left[\mathbb{V}_{\widetilde{K}|\widetilde{Y}}(\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y})\right]\right.$
$\displaystyle+\left.\mathbb{V}_{\widetilde{Y}}\left[\mathbb{E}_{\widetilde{K}|\widetilde{Y}}(\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y})\right]\right\\}.$
(7)

The first term is the average variance, due to Monte Carlo sampling of
$\widetilde{K}$, of determining whether
$\mbox{ppp}(\tilde{y})\leq\mbox{ppp}(y)$. Hence this term will become small as
$\tilde{m}$ gets large. The second term is the variance across calibration
samples $\tilde{y}$ of the probability that
$\widehat{\mbox{ppp}}(\tilde{y})\leq\mbox{ppp}(y)$. As $\tilde{m}$ gets large,
this is the variance of a Bernoulli trial with success probability
$\mbox{cppp}(y)$, i.e. $\mbox{cppp}(y)(1-\mbox{cppp}(y))$.

##  4 Computational costs and accuracy

In this section, we study the costs and benefits of different computational
allocations for the purpose of estimating $\widehat{\mbox{cppp}}(y)$. To do
so, we skip the entire modeling exercise and assume that the null distribution
of the ppp statistic is $\mbox{ppp}(\widetilde{Y})\sim\mbox{Beta}(a,b)$, a
Beta distribution with parameters $a>0$, $b>0$. Although there are no actual
$y$, $y^{*(i)}$ or $\tilde{y}^{(j)}$ values, we retain them in the notation to
indicate the roles of different variables. Under this scenario,
$\mbox{cppp}(y)$ is the tail area of the Beta distribution determined by the
observed $\mbox{ppp}(y)$: $\mbox{cppp}(y)=F_{\mbox{Beta}}(\mbox{ppp}(y),a,b)$,
where $F_{\mbox{Beta}}(\cdot;a,b)$ indicates the cumulative density function
of a $\mbox{Beta}(a,b)$ (see Figure 3 for some examples). We also assume that
independent draws can be made from the posterior and hence the posterior
predictive distribution.

Figure 3: Examples of different scenarios for the null distribution of
$\mbox{ppp}(\widetilde{Y})$, for different choices of $\mbox{cppp}(y)$. Each
color shows a different null distribution. For a given null distribution, the
vertical lines show the ppp value (on the x-axis) corresponding to the given
cppp value for each sub-figure. For example, if (uncalibrated) ppp values
really follow a Beta(2,2) distribution, and the cppp is really 0.1, then the
ppp value will be about 0.2. The cppp procedure in effect uses Monte Carlo
draws (via posterior simulation and MCMC runs) from an overdispersed ppp
distribution (Beta-Binomial in these cases) to estimate the area under the
curve to the left of 0.2.

We are interested in how to allocate computational effort between calibration
replicates ($r$) and posterior sample size for each calibration replicate
($\tilde{m}$). We assume that posterior sampling is much more costly than any
other step so that the total computational cost is $c\approx r\tilde{m}$. For
a given computational cost, we consider different allocations to $r$ and
$\tilde{m}$. Specifically, scenarios are created as follows:

  1. 1.

Choose $a$ and $b$ to set the null distribution of
$\mbox{ppp}(\widetilde{Y})\sim\mbox{Beta}(a,b)$.

  2. 2.

Choose a value of $\mbox{cppp}(y)$ to be estimated by (5).

  3. 3.

Set the corresponding observed $\mbox{ppp}(y)$ as
$\mbox{ppp}(y)=F^{-1}_{\mbox{Beta}}(\mbox{cppp}(y);a,b)$. We assume this is
known via (4) with large number of MCMC samples for the original chain $m$.

  4. 4.

Choose $c$, $r$ and $\tilde{m}$ such that $c=r\tilde{m}$.

  5. 5.

For $j=1\ldots r$:

    1. (a)

Draw $\mbox{ppp}(\tilde{y}^{(j)})\sim\mbox{Beta(a, b)}$. This is the unknown
ppp for the $j$th calibration replicate.

    2. (b)

Draw
$\tilde{k}^{(j)}\sim\mbox{Binomial}(\tilde{m},\mbox{ppp}(\tilde{y}^{(j)}))$.
This is the result of $\tilde{m}$ independent posterior samples (with data
$\tilde{y}^{(j)}$) to estimate $\mbox{ppp}(\tilde{y}^{(j)})$ by
$\widehat{\mbox{ppp}}(\tilde{y}^{(j)})=\tilde{k}^{(j)}/\tilde{m}$.

  6. 6.

Use the $r$ draws $\tilde{k}^{(j)}$ to calculate $\widehat{\mbox{cppp}}(y)$
using (5).

This procedure aims at reconstructing the key quantities used in the Monte
Carlo estimation of $\mbox{cppp}(y)$ in a scenario where the null distribution
of the ppp statistics is known, to compare estimation results at different
computational costs. Since the marginal distribution of $\tilde{k}^{(j)}$ is
Beta-Binomial with parameters $(\tilde{m},a,b)$, we do not need to sample
values, but we can compute the bias and variance of the
$\widehat{\mbox{cppp}}(y)$ analytically. The expected value of (5) corresponds
to $F_{\mbox{Beta-Binomial}}(\tilde{m}\mbox{ppp}(y);\tilde{m},a,b)$, so that
the bias with respect the true $\mbox{cppp}(Y)$ is

$F_{\mbox{Beta-
Binomial}}(\tilde{m}\mbox{ppp}(y);\tilde{m},a,b)-F_{\mbox{Beta}}(\mbox{ppp}(y);a,b),$

while the variance corresponds to

$\mathbb{V}[\widehat{\mbox{cppp}}(y)]=\frac{1}{r}F_{\mbox{Beta-
Binomial}}(\tilde{m}\mbox{ppp}(y);\tilde{m},a,b)\left[1-F_{\mbox{Beta-
Binomial}}(\tilde{m}\mbox{ppp}(y);a,b)\right].$

We calculate the bias and variance in (6)-(3.1) for different combinations of
$(r,\tilde{m})$ and choices of $\\{a,b,\mbox{cppp}(y)\\}$ controlling the null
ppp distribution. Some representative results are shown in Figure 4. We
consider different values for the computational cost $c\in\\{5,000;20,000\\}$.
Values for the number of posterior samples $\tilde{m}$ of each calibration
replicate are fixed to $\\{10,20,50,100,200,500,1000\\}$, with values for the
number of calibration replicates $r$ defined correspondingly as $c/\tilde{m}$.
Plots in Figure 4 show absolute bias, standard error, and root-mean-squared
error (RMSE = $\sqrt{\mbox{bias}^{2}+\mbox{variance}}$) of
$\widehat{\mbox{cppp}}(y)$ under two different scenarios for the null ppp
distribution and four different values of $\mbox{cppp}(y)$.

Figure 4: Absolute bias (top row), standard error (mid row) and RMSE (bottom
row) for $\widehat{\mbox{cppp}}(y)$ considering different values of
$\mbox{cppp}(y)\in\\{0.01,0.05,0.20,0.50\\}$ and two different scenarios for
the null distribution of $\mbox{ppp}(\widetilde{Y})$. In the first column
$\mbox{ppp}(\widetilde{Y})\sim\mbox{Beta}(4,2)$ while in the second
$\mbox{ppp}(\widetilde{Y})\sim\mbox{Beta}(2,2)$.

Results from these experiments suggest that values of $\tilde{m}$ in the range
$50$–$200$ can be a good rule of thumb for choosing the number of posterior
samples for each calibration replicate ($\tilde{m}$) in order to minimize the
RMSE of $\widehat{\mbox{cppp}}(y)$. Notice that the bias does not depend on
the number of calibration replicates $r$, but it is instead controlled by
$\tilde{m}$. Very low values of $\tilde{m}$ can result in large bias, but even
small values such as $50$–$200$ can make the bias negligible (see top row of
Figure 4). Then the RMSE is dominated by the variance, and it is better to
allocate further computational effort to calibration replicates rather than
posterior samples $\tilde{m}$. In the case of MCMC samples rather than
independent samples, the values of $\tilde{m}$ in these results are to be
interpreted roughly as the average ESS value of the chain
$\\{\mathbb{I}^{(i)}\\},i=1,\ldots,m$ across calibration replicates. Hence,
depending on the mixing, one may need a larger number of MCMC samples to
achieve a corresponding ESS of $50$–$200$.

Even moderately large RMSE can be acceptable to determine whether a model will
not be rejected if $\mbox{cppp}(y)$ is not small, e.g., $0.20$ or $0.50$. If
$\tilde{m}$ is between $50$–$200$, then moderately large RMSE corresponds to a
relatively small number of calibration replicates $r$, for example, $r=100$.
Compare this to a naive implementation of $\mbox{cppp}(y)$. That might use
MCMC sample size of $m={10}^{3}$–${10}^{5}$ for the real data and the same
MCMC sample size for each calibration replicate ($\tilde{m}=m$). Most of the
computational effort beyond $\tilde{m}$ = $50$–$200$ samples contributes
little to the precision of $\widehat{\mbox{cppp}}(y)$, and the total
computational cost might be $1$–$3$ orders of magnitude lower by avoiding this
naive approach. However, for very small values of $\mbox{cppp}(y)$, for
example, close to common threshold values for statistical tests ($0.01,0.05$),
one may want to be conservative and choose higher values for $\tilde{m}$.

##  5 Estimating the standard error of $\widehat{\mbox{cppp}}(y)$

In practice, one needs to estimate whether $\widehat{\mbox{cppp}}(y)$ has been
approximated with sufficient precision to reach a clear interpretation. If
$\tilde{m}$ is sufficiently large to make the bias in (6) negligible, the
standard error of $\widehat{\mbox{cppp}}(y)$ reduces to
$\sqrt{\mathbb{V}[\widehat{\mbox{cppp}}(y)]}$. In this section, we focus on
how to estimate $\mathbb{V}[\widehat{\mbox{cppp}}(y)]$ in (3.1) in practice,
illustrating two possible estimators: plug-in and bootstrap. In either case,
the key challenge is to estimate
$\mathbb{V}[\widetilde{K}|\widetilde{Y}]=\mathbb{V}_{\widetilde{K}|\widetilde{Y}}(\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y})$
in the first term of (3.1) because it involves effective sample size with a
potentially a small MCMC sample.

###  5.1 Plug-in estimator

We obtain a plug-in estimator for the cppp variance by finding plug-in
estimators for the two terms in (3.1). The first term is the average variance,
due to Monte Carlo sampling of $\widetilde{K}$, of determining whether
$\mbox{ppp}(\tilde{y})\leq\mbox{ppp}(y)$. Notice that, given
$\widetilde{Y}=\tilde{y}$,
$\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}$ is a
Bernoulli random variable with probability of success
$\text{Pr}(\widetilde{K}\leq\tilde{m}\mbox{ppp}(y))=F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))$
where $F_{\widetilde{K}}(\cdot)$ is the cumulative density function of
$\widetilde{K}$. Then the variance of this indicator variable can be expressed
as:

$\mathbb{V}_{\widetilde{K}|\tilde{y}}\left(\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\widetilde{Y}\right)=F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))\left[1-F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))\right].$

As noted in Section 3, $\widetilde{K}$ is a random variable with mean
$\tilde{m}\mbox{ppp}(\tilde{y})$ and variance
$\tilde{\tau}\tilde{m}\mbox{ppp}(\tilde{y})(1-\mbox{ppp}(\tilde{y}))$, where
$\tilde{\tau}$ is the integrated autocorrelation time for an MCMC chain using
data $\tilde{y}$. Hence $\widetilde{K}$ is like a binomial random variable
with variance inflated by $\tilde{\tau}$. We use a normal approximation with a
continuity correction for the distribution of $\widetilde{K}$, so that

$F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))\approx\hat{F}_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))=F_{\mbox{N}}\left(\tilde{m}\mbox{ppp}(y)+\frac{1}{2};\tilde{m}\mbox{ppp}(\tilde{y}),\tilde{\tau}\tilde{m}\mbox{ppp}(\tilde{y})\big{(}1-\mbox{ppp}(\tilde{y})\big{)}\right)$
(8)

where $F_{\mbox{N}}(\cdot;\mu,\sigma^{2})$ denotes the cumulative distribution
function of a Normal distribution with mean $\mu$ and variance $\sigma^{2}$.
In practice, we can use $\widehat{\mbox{ppp}}(\tilde{y}^{(j)})$ as a plug-in
estimate for $\mbox{ppp}(\tilde{y})$ for each calibration replicate
$j=1,\ldots,r$ so that the first term of $\eqref{eq:cppp_variance}$ is
approximated via Monte Carlo as

$\mathbb{E}_{\widetilde{Y}}\left[\mathbb{V}[\widetilde{K}|\widetilde{Y}]\right]\approx\frac{1}{r}\sum_{j=1}^{r}\hat{F}_{\widetilde{K}}^{(j)}\left[1-\hat{F}_{\widetilde{K}}^{(j)}\right]$

where
$\hat{F}_{\widetilde{K}}^{(j)}=F_{\mbox{N}}\left(\tilde{m}\mbox{ppp}(y)+\frac{1}{2};\tilde{m}\widehat{\mbox{ppp}}(\tilde{y}^{(j)}),\hat{\tilde{\tau}}^{(j)}\tilde{m}\widehat{\mbox{ppp}}(\tilde{y}^{(j)})(1-\widehat{\mbox{ppp}}(\tilde{y}^{(j)}))\right)$.
The challenge here is to estimate the autocorrelation time
$\tilde{\tau}^{(j)}$ for each calibration replicate, as we would like to use a
small number of MCMC samples; we propose an estimator in 5.3.

For the second term in $\eqref{eq:cppp_variance}$, we have that the inner
expectation, given $\widetilde{Y}=\tilde{y}$, is

$\mathbb{E}_{\widetilde{K}|\tilde{y}}\left(\mathbb{I}\left\\{\widetilde{K}\leq\tilde{m}\mbox{ppp}(y)\right\\}|\tilde{y}\right)=\text{Pr}(\widehat{\mbox{ppp}}(\tilde{y})\leq\mbox{ppp}(y))=F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y)).$

Note that for large $\tilde{m}$ and/or $\mbox{ppp}(\tilde{y})$ far from
$\mbox{ppp}(y)$,
$F_{\widetilde{K}}(\tilde{m}\mbox{ppp}(y))\approx\mathbb{I}\left\\{\mbox{ppp}(\tilde{y})\leq\mbox{ppp}(y)\right\\}$.

Using the normal approximation for $F_{\widetilde{K}}(\cdot)$ and Monte Carlo
approximation of the outer variance, and combining the two terms of (3.1)
gives the final plug-in estimator as

$\displaystyle\widehat{\mathbb{V}}[\widehat{\mbox{cppp}}(y)]$
$\displaystyle=\frac{1}{r}\sum_{j=1}^{r}\hat{F}_{\widetilde{K}}^{(j)}\left[1-\hat{F}_{\widetilde{K}}^{(j)}\right]+\frac{1}{r}\sum_{j=1}^{r}\left[\left(\hat{F}_{\widetilde{K}}^{(j)}-\overline{F_{\widetilde{K}}}\right)^{2}\right]$
$\displaystyle=\overline{F_{\widetilde{K}}}(1-\overline{F_{\widetilde{K}}})$
(9)

where
$\overline{F_{\widetilde{K}}}=\sum_{j=1}^{r}\hat{F}_{\widetilde{K}}^{(j)}/r$.
Notice that, as $\tilde{m}\rightarrow\infty$:
$\hat{F}^{(j)}_{\widetilde{K}}\rightarrow\mbox{ppp}(\tilde{y})$ and
$\overline{F_{\widetilde{K}}}\rightarrow\widehat{\mbox{cppp}}(y)$.

###  5.2 Bootstrap estimators

A bootstrap procedure can be used to estimate the variance of
$\widehat{\mbox{cppp}}(y)$. The main idea is to obtain $b$ bootstrap estimates
$\\{\widehat{\mbox{cppp}}(y)^{(l)}\\}_{l=1,\ldots,b}$, and evaluate the sample
variance of these estimates. To do so we resample the calibration replicates
and MCMC samples used for the original $\mbox{cppp}(y)$ estimate, taking into
account the structure of the Monte Carlo estimation procedure in Figures 1-2.

The idea is the same as bootstrapping hierarchical data, where observations
are associated with groups, for example, students belonging to schools. A
bootstrap of this kind of data is typically performed in two stages to retain
the hierarchical structure of the data: first groups are resampled, then
observations are resampled within the groups. In our case, we first resample
the calibration replicates, and then the MCMC samples within each replicate.
In practice, instead of resampling all the MCMC outputs, we resample only the
derived output
$\left\\{\Delta(y^{*(i)},\theta^{(i)}|\tilde{y})\right\\},i=1,\ldots,\tilde{m}$,
for a calibration replicate $\tilde{y}$. Only these derived quantities are
needed to obtain the count $\tilde{k}$. In this step we need to retain serial
dependence between the MCMC samples, so we can either use moving block
bootstrap (MBB) (Mignani and Rosa, 1995) or exploit the normal approximation
for $\widetilde{K}$.

In the following, we summarize the steps to obtain these two bootstrap
estimates. For a fixed number of bootstrap samples $b$, each
$\widehat{\mbox{cppp}}(y)^{(l)}$, for $l=1,\ldots,b$, is obtained as follows:

  1. 1.

Sample with replacement $r$ calibration replicates from $S_{3}(y,r)$. We
denote with $S_{3}^{(l)}(y,r)$ the $l$-th bootstrap sample.

  2. 2.

For each calibration replicate $j$ in $S_{3}^{(l)}(y,r)$ perform one of the
following:

     * a)

Bootstrap-MBB: Calculate $\\{\tilde{k}^{(j)}\\}$ from a moving block bootstrap
(re)sample of the chain
$\left\\{\Delta(y^{*(i)},\theta^{(i)}|\tilde{y}^{(j)})\right\\},i=1,\ldots,m$.
Then use counts $\\{\tilde{k}^{(j)}\\}^{(l)}$ to obtain
$\widehat{\mbox{cppp}}(y)^{(l)}$ as in (5).

     * b)

Bootstrap-normal: draw $\\{\tilde{k}^{(j)}\\}$ from the normal approximation
with mean and variance as stated in (8).

  3. 3.

Use counts $\\{\tilde{k}^{(j)}\\}^{(l)}$, $j=1,\ldots,r$, to obtain
$\widehat{\mbox{cppp}}(y)^{(l)}$ as in (5).

Although the moving block bootstrap accounts for the original correlation
structure between MCMC samples, the number and sizes of blocks are additional
parameters to choose in comparison with the normal approximation.
Nevertheless, in both cases, one has to choose the number of bootstrap
replications which adds to the $\mbox{cppp}(y)$ computational cost.

###  5.3 A transfer estimator for the integrated autocorrelation time

For the normal approximation used in the first term of the plug-in estimate
and the boostrap-normal approach, we need an estimate of the integrated
correlation time $\tilde{\tau}$, or the corresponding
$\widetilde{\mbox{ESS}}=\tilde{m}/\tilde{\tau}$ for the Markov chain
$\left\\{\tilde{\mathbb{I}}^{(i)}\right\\}$, $i=1\ldots\tilde{m}$. The index
$(j)$ for the calibration replicate will be omitted in this section.

Each $\tilde{\mathbb{I}}^{(i)}$ is an indicator value from the posterior
predictive sample for calibration replicate $\tilde{y}$. That is,
$\tilde{\mathbb{I}}^{(i)}=\mathbb{I}\left\\{\Delta(y^{*(i)},\theta^{(i)}|\tilde{y})\geq
0\right\\}$, where
$S_{2}(\tilde{y},\tilde{m})=\left\\{(y^{*(i)},\theta^{(i)})\right\\}i=1,\ldots,\tilde{m}$
is the MCMC output from sampling $(y^{*(i)},\theta^{(i)})\sim
p(y,\theta|\tilde{y})$. However, methods for estimating $\tilde{\tau}$ or
$\widetilde{\mbox{ESS}}$ typically need long chains, whereas here we seek to
reduce computation by using relatively short chains for calibration
replicates.

On heuristic grounds, we propose to use information from the single long MCMC
chain for the real data to provide reasonable estimates of $\tilde{\tau}$ for
the calibration replicates. Specifically, we assume that the autocorrelation
structure of the chain
$\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)\right\\},i=1,\ldots,m$ will be
similar for different values of $y$. This is the chain of differences between
posterior predictive discrepancies and data discrepancies. It is important to
use this chain rather than the chain of indicators,
$\left\\{\tilde{\mathbb{I}}^{(i)}\right\\}$, because the correlation structure
of the latter will clearly depend on $\tilde{y}$. It is a chain of 0s and 1s
with probability $\mbox{ppp}(\tilde{y})$ of being 1, so for some values of
$\tilde{y}$ the $\tilde{\mathbb{I}}^{(i)}$ may be mostly 0s, for others mostly
1s, and so on. To reflect this dependence, we write $\tilde{\tau}$ as
$\tilde{\tau}(\mbox{ppp}(\tilde{y}))$.

The strategy will be as follows. For an estimated ppp for a calibration
replicate, $q=\widehat{\mbox{ppp}}(\tilde{y})$, we will use the shifted chain
$\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)-\Delta_{q}\right\\}$, where
$\Delta_{q}$ is the $q$ quantile of
$\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)\right\\}$. This satisfies
$\frac{1}{m}\sum\mathbb{I}\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)-\Delta_{q}\leq
0\right\\}=q$ and thus corresponds to a chain that yields an estimated ppp
equal to that of the calibration replicate.

Then, if the above assumption is reasonable, the indicators of the shifted
chain,

$\mathbb{I}\left\\{\Delta(y^{*(i)},\theta^{(i)}|y)\leq\Delta_{q}\right\\},$
(10)

should have integrated autocorrelation time close to that of
$\tilde{\tau}(q)$, i.e. of chain of indicators for the calibration replicate.
The ESS (and hence $\tilde{\tau}$) for the chain in (10) can be estimated by
any suitable method for a discrete Markov chain such as the batch mean
estimator of Flegal and Jones (2010), implemented in the mcmcse package
(Flegal et al., 2021, see function mcse).

##  6 Examples

In this section, we consider two examples that have been used in the
literature (Hjort et al., 2006). The first example uses a normal model for
data with outliers, while the second comprises two versions of a classic
example of a capture-recapture model. To estimate the models used in each
example we make use of the NIMBLE software (de Valpine et al., 2022) a
flexible R-based system for hierarchical models. The procedure for cppp is
implemented using NIMBLE’s algorithm programming system and can be used with
other models and discrepancies provided by the user.

These examples are difficult because the cppp for all cases are near or below
the traditional rejection threshold of 0.05. Since the most dramatic
computational gains of our proposed methods will be relevant when the cppp is
well above such a threshold (when a model is acceptable) we also include a
version of the capture-recapture example with simulated data.

We use these examples to i) assess how results from the simulated experiments
in Section 4 hold in real-data applications and ii) compare the different
estimators for the cppp standard error illustrated in Section 5. To do so, we
look at $\widehat{\mbox{cppp}}(y)$ and its standard error under different
scenarios, considering the computational cost, its allocation and mixing
properties of the original MCMC run. We show averages of
$\widehat{\mbox{cppp}}(y)$ and each kind of standard error estimate from $500$
runs of the whole procedure (Figures 5-7). In the Supplementary Materials, we
also report coverage for approximate confidence intervals at 95% (Table S1)
and estimates of the ppp distribution for each example (Figure S1).

Since we are working with real data examples, we do not have a ground truth to
compare with. Hence, we use brute force computing to obtain good estimates of
the truth for comparison. We obtain brute force estimates of cppp by choosing
large $r$ and $\tilde{m}$. We obtain brute force Monte Carlo estimates of the
standard error of cppp as the standard deviation of $\widehat{\mbox{cppp}}(y)$
from $500$ repeats of the entire cppp computation (illustrated in Figure 2),
for each scenario of $r$ and $\tilde{m}$ considered.

###  6.1 Normal model example using newcomb data

A standard example to illustrate posterior predictive p-values in the
literature uses the speed of light data (see for example Gelman et al., 2013a,
Section 6.3), that comprises $66$ measurements of the speed of light made by
Simon Newcomb in $1882$. This is a classic example of a dataset with outliers.
A standard model for the data would use a normal distribution with unknown
mean $\mu$ and variance $\sigma^{2}$, however, there are two extreme low
measurements not compatible with the model. As in Gelman et al. (2013a), we
consider as discrepancy measure $D(y,\mu)=|y_{(61)}-\mu|-|y_{(6)}-\mu|$, where
$y_{(i)}$ denotes the $i$th-ordered data point. Assuming a uniform prior for
$(\mu,\log\sigma)$, Gelman et al. (2013a) report
$\widehat{\mbox{ppp}}(y)=0.26$ using $200$ simulations from the posterior
predictive. We replicated the analysis considering the same model, using
$m=4,000$ samples after a burn-in of $1,000$ samples, obtaining
$\widehat{\mbox{ppp}}(y)=0.205$. This result is in line with Hjort et al.
(2006), who reported a value of $0.208$ using 1 million samples. Although the
authors did not compute a value for the $\widehat{\mbox{cppp}}(y)$, they
speculate that it falls on “the statistical borderline of surprise”; we in
fact estimate $\widehat{\mbox{cppp}}(y)=0.055$ using $r=\tilde{m}=1,000$.

Similarly to the experiments in Section 4, we consider $4$ values for the
total computational cost $c\in\\{5,000;10,000;20,000;50,000\\}$, fix values of
$\tilde{m}\in\\{50,100,200,500\\}$, and define $r=c/\tilde{m}$. We also
consider two different conditions for the mixing of the original MCMC. In the
first case, we have fairly good mixing for the model parameters, where
$\mbox{ESS}/m$ for $(\sigma,\mu)$ is $(0.2,1)$; in the second one, we used
random walk Metropolis-Hastings samplers for posterior sampling and purposely
tune them to induce bad mixing, so that the resulting ESS is 4 to 6 times
lower than in the first case, i.e. $\mbox{ESS}/m$ for $(\sigma,\mu)$ is (0.05,
0.15).

Figure 5: Comparison of estimates of $\widehat{\mbox{cppp}}(y)$ and its
standard error for different combinations of calibration replicates $r$ and
MCMC samples $\tilde{m}$ for fixed computational cost $c$, under good mixing
conditions (top row) and bad ones (bottom row). Error bars correspond to one
standard deviation, estimated using different methods. The “Monte Carlo”
(black) case shows brute force estimate of average $\widehat{\mbox{cppp}}(y)$
and its standard error. The dashed gray line shows brute force estimate of the
correct $\mbox{cppp}(y)$ while the shaded gray area is its Monte Carlo
standard deviation.

Figure 5 shows averages (from $500$ runs) of $\widehat{\mbox{cppp}}(y)$ and of
each standard error estimator, for different mixing conditions (subfigure
rows), fixed computational costs (subfigures), and combinations of calibration
replicates and MCMC samples ($r$ and $\tilde{m}$ values). Error bars represent
one standard deviation. The colored error bars show the different estimators
described in Section 5, while the black bar denotes the brute force Monte
Carlo estimate. Results for the bootstrap estimators are computed using $100$
replications. We also report, for each value of $\tilde{m}$, the average ESS
for the chain $\\{\mathbb{I}^{(i)}\\},i=1,\ldots,m$ across calibration
replicates, denoted as $\widehat{\widetilde{\mbox{ESS}}}$. We use the transfer
method described in Section 5.3 to estimate the ESS for each calibration
replicate and report the mean value averaged across calibration replicates and
the $500$ runs.

Under the good mixing conditions, we obtain comparable results to the brute
force values for the $\mbox{cppp}(y)$ estimate; in the bad mixing case, we
observe small but clear bias when either the number of MCMC samples or
calibration replicates is small (e.g., $\tilde{m}=50$ corresponding to
$\mbox{ESS}=20$ or $r\leq 50$). The different variance estimators seem to
perform similarly under both mixing conditions. However, there are relevant
differences between the estimators in terms of coverage (see Table S1).
Intervals based on the Plug-in or Bootstrap-MBB estimators are close to the
nominal 95% level when there are good compromises between ESS and the number
of calibration replicates ($\mbox{ESS}\geq 50$ and $r\geq 100$). Instead, the
bootstrap normal estimator leads more often to under-coverage, which becomes
more evident under bad mixing conditions. In settings with low ESS values and
a high number of calibration replicates, the coverage can be 0 due to biased
cppp estimates with low variance, which results in confidence intervals that
do not cover the true cppp.

###  6.2 Capture-recapture models using dipper data

Another example from the literature (Hjort et al., 2006; Nott et al., 2018),
considers capture-recapture models using data on the European Dipper (Cinclus
cinclus). Data consists of $n=294$ sighting histories collected over $k=7$
annual sighting occasions from $1980$ to $1987$. A simple capture-recapture
model can be described as a Hidden Markov Model, where each sighting is a
realization of a Bernoulli random variable
$Y_{i,t}\in\\{1=\mbox{seen},0=\mbox{not seen}\\}$, conditional to a binary
latent variable $X_{i,t}\in\\{1=\mbox{alive},0=\mbox{dead}\\}$. The model is
parametrized by the annual probability of survival
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{k-1})$ and the probability of being
seen $\boldsymbol{p}=(p_{2},\ldots,p_{k})$. The hierarchical model corresponds
to:

$\displaystyle X_{i,t}|X_{i,t-1}$
$\displaystyle\sim\mbox{Bernoulli}(\phi_{t-1}x_{i,t-1}),\quad t=2,\ldots,k,$
$\displaystyle Y_{i,t}|X_{i,t}$
$\displaystyle\sim\mbox{Bernoulli}(p_{t}x_{i,t})\quad t=2,\ldots,k$
$\displaystyle\phi_{t-1}$ $\displaystyle\sim\mbox{Uniform}(0,1)\quad
p_{t}\sim\mbox{Uniform}(0,1),\quad t=2,\ldots,k,$

where $X_{i,1}=Y_{i,1}=1$, for $i=1,\ldots,n$. A simpler version of the model
above assumes constant survival and capture probabilities over time, so that
$p_{2}=\ldots=p_{k}=p$ and $\phi_{1}=\ldots=\phi_{k-1}=\phi$. Under both
models, uniform and independent priors are chosen for the survival and capture
probabilities. Hjort et al. (2006) considers these two versions of the model,
following Brooks et al. (2000), referring to those as the _large model_ (T/T
model) and the _small model_ (C/C model) respectively.

To evaluate whether the models characterize the data well, Brooks et al.
(2000) use the Freeman-Tukey statistic (Freemen and Tukey, 1954) as a
discrepancy measure

$D(y,\theta)=\sum_{s=1}^{k-1}\sum_{t=2}^{k}(\sqrt{z_{st}}-\sqrt{e_{st}})^{2},$
(11)

where $z_{st}$ is the observed number of animals released at time $s$ and
captured at time $t$, while $e_{st}$ is the corresponding expected number.

For the original MCMC run we follow Brooks et al. (2000) and consider $10,000$
samples after $1,000$ burn-in. We obtain values for the observed ppp of
$0.064$ for the C/C model, and $0.083$ for the T/T, slightly differently from
Brooks et al. (2000) who reported $0.069$ and $0.086$ respectively. However,
we can attribute this small difference to Monte Carlo variation. Hjort et al.
(2006) reports different values, $0.060$ for the C/C model, and $0.075$ for
the T/T model; we obtain the same values when considering only $1,000$ samples
after burn-in. To compute the cppp, Hjort et al. (2006) use $500$ calibration
replicates without specifying the number of MCMC samples, obtaining cppp
values of $0.022$ and $0.002$ under the two models respectively. We obtain
slightly different values. Considering a naive estimate using $r=1,000$ and
$m=10,000$, the cppp is equal to $0.044$ for the C/C model and $0.010$ for the
$T/T$ model, similar to the Monte Carlo baseline. However, we reach the same
conclusion as in Hjort et al. (2006), “that the given data are much less
likely under the large model than under the small model, although the nominal
ppp values led Brooks et al. (2000) to the opposite conclusion.”

As in the previous example, we compare the performance of the standard
deviation estimators for the two models (Figure 6). Since for both models we
obtain significant cppps, i.e. values under the commonly used threshold of
$0.05$, we also simulated a case where the cppp is not significant, generating
data under the T/T model (Figure 7). This is important because cases where the
cppp turns out to be acceptable are the ones where the most computation can be
saved since rough estimates of cppp may be satisfactory.

Figure 6: Comparison of estimates of $\widehat{\mbox{cppp}}(y)$ and its
standard error for different combinations of calibration replicates $r$ and
MCMC samples $\tilde{m}$ for fixed computational cost $c$, under the C/C model
(top row) and T/T model (bottom row). Error bars correspond to one standard
deviation, estimated using different methods. The “Monte Carlo” (black) case
shows brute force estimate of average $\widehat{\mbox{cppp}}(y)$ and its
standard error. The dashed gray line shows brute force estimate of the correct
$\mbox{cppp}(y)$ while the shaded gray area is its Monte Carlo standard
deviation.

For both models, our approximate results are comparable to the brute force
values under some combinations of $m$ and $r$. However, we obtain biased point
estimates of the $\mbox{cppp}(y)$ when either the number of calibration
replicates or independent MCMC samples is low. In particular, for $r<50$ we
have biased estimates under the C/C model, while for the T/T model we obtain a
good point estimate when the ESS is at least $100$ ($\tilde{m}=500$) and
$r>50$.

This behavior can be related to the different shapes of the null
$\mbox{ppp}(Y)$ distributions under the two models (Figure S1). Under the C/C
model, the ppp statistic has a symmetric and almost flat distribution.
However, for the T/T model we observe a left-skewed distribution, similarly to
the two examples in Section 4. For similar values of the ESS under the two
models, we need a different number of calibration replicates to accurately
estimate the cppp. This is because when the null distribution of
$\mbox{ppp}(Y)$ is skewed and the $\mbox{cppp}(y)$ is very small, it is hard
to sample values of ppp smaller than the observed $\mbox{ppp}(y)$, (i.e. in
the tail of the distribution); hence a higher number of calibration replicates
is needed to achieve the same accuracy.

Considering the cppp standard error, we observe good performances of all
estimators at high computational cost ($c=50,0000$). The full nonparametric
bootstrap (bootstrap-MBB) is more robust to low ESS values than the bootstrap-
normal, at the cost of adding more steps in the computation. The plug-in
estimator performs quite well overall, making it the preferable option. The
same considerations apply to the performance of the intervals in terms of
coverage (Table S1).

Finally, in Figure 7 we report results for data simulated from the T/T model,
for which the brute force cppp is $0.23$. In this case, even results from the
lowest computational cost ($c=5,000$) are accurate enough to conclude that the
model will not be rejected. One hopes that models are often well chosen and
hence this substantial reduction in computational cost to conclude a model is
acceptable can be realized in many applications.

Figure 7: Comparison of estimates of $\widehat{\mbox{cppp}}(y)$ and its
standard error for different combinations of calibration replicates $r$ and
MCMC samples $\tilde{m}$ for fixed computational cost $c$ for a simulated
example using the T/T model. Error bars correspond to one standard deviation,
estimated using different methods. The “Monte Carlo” (black) case shows brute
force estimate of average $\widehat{\mbox{cppp}}(y)$ and its standard error.
The dashed gray line shows brute force estimate of the correct
$\mbox{cppp}(y)$ while the shaded gray area is its Monte Carlo standard
deviation.

##  7 Discussion

This paper gives computational methods to efficiently assess the goodness of
fit of Bayesian models via calibrated posterior predictive p-values. Our
proposal allows obtaining the cppp much faster than via naive implementation,
i.e., when the calibration procedure uses the same number of iterations as in
the original MCMC. The main takeaway is that good choices for the number of
calibration replicates vs. MCMC samples per replicate can reduce the
computational burden by orders of magnitude in model assessment without
significant loss in accuracy. These two quantities control the bias and
variance trade-off of the cppp estimate. A good rule of thumb for the number
of MCMC samples is to consider a number of iterations such that the ESS is
within $50$–$200$, making the bias negligible, and allocate further
computational effort to calibration replicates—which can be readily
parallelized—to reduce the variance. In practice, the greatest computational
savings will be obtained when one can determine that at model is acceptable
based on rough cppp with accurately quantified uncertainty, whereas the most
computation will be demanded if the cppp is near a decision threshold of
interest.

We also describe different methods to quantify the uncertainty associated with
the cppp estimate that rely on plug-in estimators or bootstrap procedures. We
find that the plug-in estimate for the cppp standard error performs generally
well in terms of accuracy, and is preferred because it does not add
computational effort as do the bootstrap procedures. This is important in that
accurate and fast uncertainty quantification for the cppp estimate allows for
immediate improvements in the procedure. When the cppp is not small, even a
rough estimate of the cppp obtained at low computational cost, is enough to
conclude that the model will not be rejected. As a consequence, we can make
the procedure adaptive, including stopping rules accounting for this
situation.

The proposed procedure to approximate the cppp is implemented using the NIMBLE
software, which is a flexible platform that can handle a wide range of models
and discrepancies. This software provides users with the tools to efficiently
implement our proposed procedure, allowing the use of the cppp methodology
beyond toy examples.

SUPPLEMENTARY MATERIAL

Supplementary figures and tables:

A .pdf file containing supplementary figures and tables.

## Supplementary figures and tables

computational cost | 5,000 | 10,000 | 20,000 | 50,000
---|---|---|---|---
Calibration replicates $r$ | 100 | 50 | 25 | 10 | 200 | 100 | 50 | 20 | 400 | 200 | 100 | 40 | 1000 | 500 | 250 | 100
MCMC samples $\tilde{m}$ | 50 | 100 | 200 | 500 | 50 | 100 | 200 | 500 | 50 | 100 | 200 | 500 | 50 | 100 | 200 | 500
Newcomb example - good mixing
Average ESS | 49 | 99 | 198 | 496 | 49 | 99 | 198 | 496 | 49 | 99 | 198 | 496 | 49 | 99 | 198 | 496
Plug-in | 0.988 | 0.970 | 0.774 | 0.568 | 0.998 | 0.958 | 0.952 | 0.708 | 1.000 | 0.982 | 0.930 | 0.886 | 1.000 | 0.984 | 0.982 | 0.938
Bootstrap - MBB | 0.988 | 0.982 | 0.890 | 0.556 | 0.970 | 0.982 | 0.976 | 0.784 | 0.940 | 0.978 | 0.970 | 0.924 | 0.780 | 0.926 | 0.978 | 0.966
Bootstrap - Normal | 0.942 | 0.960 | 0.844 | 0.594 | 0.874 | 0.974 | 0.960 | 0.740 | 0.628 | 0.984 | 0.966 | 0.892 | 0.000 | 0.962 | 0.968 | 0.946
Newcomb example - bad mixing
Average ESS | 20 | 41 | 81 | 204 | 20 | 41 | 82 | 203 | 20 | 41 | 82 | 203 | 20 | 41 | 82 | 203
Plug-in | 0.984 | 0.966 | 0.850 | 0.636 | 0.964 | 0.984 | 0.930 | 0.752 | 0.912 | 0.992 | 0.956 | 0.888 | 0.000 | 0.996 | 0.986 | 0.954
Bootstrap - MBB | 0.918 | 0.980 | 0.858 | 0.566 | 0.794 | 0.962 | 0.944 | 0.788 | 0.506 | 0.950 | 0.978 | 0.926 | 0.000 | 0.868 | 0.986 | 0.960
Bootstrap - Normal | 0.546 | 0.950 | 0.886 | 0.670 | 0.174 | 0.878 | 0.960 | 0.858 | 0.002 | 0.686 | 0.974 | 0.918 | 0.000 | 0.162 | 0.950 | 0.958
Dipper example - C/C model
Average ESS | 19 | 39 | 77 | 197 | 19 | 39 | 77 | 196 | 19 | 39 | 78 | 196 | 19 | 39 | 77 | 196
Plug-in | 0.996 | 0.970 | 0.884 | 0.610 | 0.998 | 0.996 | 0.928 | 0.724 | 1.000 | 1.000 | 0.964 | 0.836 | 1.000 | 1.000 | 0.984 | 0.940
Bootstrap - MBB | 0.962 | 0.976 | 0.842 | 0.462 | 0.924 | 0.958 | 0.970 | 0.700 | 0.802 | 0.968 | 0.986 | 0.904 | 0.248 | 0.974 | 0.984 | 0.950
Bootstrap - Normal | 0.720 | 0.980 | 0.908 | 0.636 | 0.256 | 0.898 | 0.962 | 0.866 | 0.014 | 0.724 | 0.976 | 0.910 | 0.000 | 0.262 | 0.950 | 0.968
Dipper example - T/T model
Average ESS | 13 | 25 | 49 | 120 | 13 | 25 | 49 | 119 | 13 | 25 | 49 | 120 | 13 | 25 | 49 | 120
Plug-in | 0.994 | 1.000 | 0.958 | 0.352 | 0.970 | 1.000 | 0.970 | 0.552 | 0.864 | 0.990 | 0.982 | 0.636 | 0.000 | 0.962 | 0.972 | 0.776
Bootstrap - MBB | 0.882 | 0.938 | 0.672 | 0.164 | 0.534 | 0.970 | 0.880 | 0.350 | 0.054 | 0.868 | 0.984 | 0.622 | 0.000 | 0.350 | 0.846 | 0.860
Bootstrap - Normal | 0.256 | 0.926 | 0.902 | 0.214 | 0.000 | 0.842 | 0.972 | 0.472 | 0.000 | 0.362 | 0.986 | 0.762 | 0.000 | 0.000 | 0.722 | 0.834
Simulated example - T/T model
Average ESS | 44 | 88 | 175 | 441 | 44 | 88 | 176 | 444 | 44 | 88 | 175 | 439 | 44 | 88 | 175 | 440
Plug-in | 0.964 | 0.962 | 0.918 | 0.912 | 0.986 | 0.968 | 0.932 | 0.932 | 0.998 | 0.984 | 0.958 | 0.946 | 1.000 | 0.998 | 0.980 | 0.962
Bootstrap - MBB | 0.968 | 0.970 | 0.936 | 0.900 | 0.974 | 0.964 | 0.948 | 0.928 | 0.986 | 0.968 | 0.962 | 0.934 | 1.000 | 0.992 | 0.972 | 0.950
Bootstrap - Normal | 0.952 | 0.958 | 0.916 | 0.874 | 0.964 | 0.964 | 0.944 | 0.916 | 0.968 | 0.954 | 0.960 | 0.922 | 1.000 | 0.982 | 0.962 | 0.934

Table S1:  Confidence interval coverage for the examples in Section 6 for each
scenario of computational cost, number of calibration replicates and MCMC
samples. We use normality approximation to calculate the confidence intervals
at 95% as $\widehat{\mbox{cppp}}(y)\pm
1.96\sqrt{\widehat{\mathbb{V}}[\widehat{\mbox{cppp}}(y)]}$

.

Figure S1: Histogram and density plot of estimated ppps obtained using
$r=m=1,000$ for all the examples in Section 6. The red line indicates the
observed value $\mbox{ppp}(y)$.

## References

  * Bayarri and Berger (2000) Bayarri, M., and Berger, J. O. (2000), “P values for composite null models (with discussion),” Journal of the American Statistical Association, 95, 1127–1142.
  * Bayarri and Castellanos (2007) Bayarri, M., and Castellanos, M. (2007), “Bayesian checking of the second levels of hierarchical models,” Statistical Science, 22, 322–343.
  * Box (1980) Box, G. E. (1980), “Sampling and Bayes’ inference in scientific modelling and robustness,” Journal of the Royal Statistical Society: Series A, 143, 383–404.
  * Brooks et al. (2000) Brooks, S. P., Catchpole, E. A., Morgan, B. J. et al. (2000), “Bayesian animal survival estimation,” Statistical Science, 15, 357–376.
  * de Valpine et al. (2022) de Valpine, P., Paciorek, C., Turek, D., Michaud, N., Anderson-Bergman, C., Obermeyer, F., Wehrhahn Cortes, C., Rodríguez, A., Temple Lang, D., Paganin, S., and Hug, J. (2022), “NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling,” .
  * de Valpine et al. (2017) de Valpine, P., Turek, D., Paciorek, C. J., Anderson-Bergman, C., Lang, D. T., and Bodik, R. (2017), “Programming With Models: Writing Statistical Algorithms for General Model Structures With NIMBLE,” Journal of Computational and Graphical Statistics, 26, 403–413.
  * Flegal et al. (2021) Flegal, J. M., Hughes, J., Vats, D., Dai, N., Gupta, K., and Maji, U. (2021), mcmcse: Monte Carlo Standard Errors for MCMC, Riverside, CA, and Kanpur, India, r package version 1.5-0.
  * Flegal and Jones (2010) Flegal, J. M., and Jones, G. L. (2010), “Batch means and spectral variance estimators in Markov chain Monte Carlo,” The Annals of Statistics, 38, 1034 – 1070.
  * Freemen and Tukey (1954) Freemen, M., and Tukey, J. W. (1954), “Transformations related to the angular and the square root,” The Annals of Mathematical Statistics, 21, 607–611.
  * Gelman et al. (2013a) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013a), Bayesian Data Analysis, CRC press.
  * Gelman et al. (1996) Gelman, A., Meng, X.-L., and Stern, H. (1996), “Posterior predictive assessment of model fitness via realized discrepancies,” Statistica Sinica, 733–760.
  * Gelman et al. (2013b) Gelman, A. et al. (2013b), “Two simple examples for understanding posterior p-values whose distributions are far from uniform,” Electronic Journal of Statistics, 7, 2595–2602.
  * Hjort et al. (2006) Hjort, N. L., Dahl, F. A., and Steinbakk, G. H. (2006), “Post-processing posterior predictive p values,” Journal of the American Statistical Association, 101, 1157–1174.
  * Johnson (2007) Johnson, V. E. (2007), “Bayesian model assessment using pivotal quantities,” Bayesian Analysis, 2, 719–733.
  * Meng (1994) Meng, X.-L. (1994), “Posterior predictive $p$-values,” The Annals of Statistics, 22, 1142–1160.
  * Mignani and Rosa (1995) Mignani, S., and Rosa, R. (1995), “The moving block bootstrap to assess the accuracy of statistical estimates in Ising model simulations,” Computer Physics Communications, 92, 203–213.
  * Nott et al. (2018) Nott, D. J., Drovandi, C. C., Mengersen, K., Evans, M. et al. (2018), “Approximation of Bayesian Predictive $p$-Values with Regression ABC,” Bayesian Analysis, 13, 59–83.
  * Robert and Casella (1999) Robert, C. P., and Casella, G. (1999), Monte Carlo statistical methods, vol. 2, Springer.
  * Robins et al. (2000) Robins, J. M., van der Vaart, A., and Ventura, V. (2000), “Asymptotic distribution of p values in composite null models,” Journal of the American Statistical Association, 95, 1143–1156.
  * Rubin (1984) Rubin, D. B. (1984), “Bayesianly Justifiable and Relevant Frequency Calculations for the Applied Statistician,” The Annals of Statistics, 12, 1151–1172.
  * Steinbakk and Storvik (2009) Steinbakk, G. H., and Storvik, G. O. (2009), “Posterior Predictive p-values in Bayesian Hierarchical Models,” Scandinavian Journal of Statistics, 36, 320–336.
  * Yuan and Johnson (2012) Yuan, Y., and Johnson, V. E. (2012), “Goodness-of-fit diagnostics for Bayesian hierarchical models,” Biometrics, 68, 156–164.