Simoiu 2017 - Threshold Test Generative Model

Summary

Each stop of a race-$r$ driver by department $d$ draws a signal $p\sim \text{Beta}({\varphi }_{rd},{\lambda }_{rd})$, a search occurs iff $p\ge t_{rd}$, and a hit occurs with probability $p$. Search and hit counts are binomial with probabilities $p_{rd}=1-I_{t_{rd}}({\varphi }_{rd},{\lambda }_{rd})$ and $q_{rd}={\varphi }_{rd}(1-I_{t_{rd}}({\mu }_{rd},{\lambda }_{rd}+1))/(1-I_{t_{rd}}({\varphi }_{rd},{\lambda }_{rd}))$. Identification comes from additive race and department structure on $\varphi ,\lambda$ plus hierarchical priors on thresholds.

Overview

Observed per stop: race $r_i$, department $d_i$, search $S_i\in \{0,1\}$, hit $H_i\le S_i$. Goal: race-specific thresholds per department. Bayesian benefits: uncertainty without bootstrapping, and pooling across departments for rarely stopped groups.

Main Content

Generative process (Sec. 2.3)

Signal: $p_i\sim \text{Beta}({\varphi }_{r_i{d}_i},{\lambda }_{r_i{d}_i})$ with mean $\varphi \in (0,1)$ and total count $\lambda >0$ (i.e. $\alpha =\varphi \lambda$, $\beta =(1-\varphi )\lambda$). Search: $S_i=1$ iff $p_i\ge t_{r_i{d}_i}$. Hit: $H_i\sim \text{Bernoulli}(p_i)$ if $S_i=1$, else $H_i=0$. Structure (Eqs. 1-2): ${\varphi }_{rd}={\text{logit}}^{-1}({\varphi }_r+{\varphi }_d)$, ${\lambda }_{rd}=\exp ({\lambda }_r+{\lambda }_d)$, with ${\varphi }_d={\lambda }_d=0$ for the largest department. This gives $2(D+R-1)$ parameters for $DR$ signal distributions.

Priors

${\varphi }_r,{\lambda }_r\sim N(0,2)$; ${\varphi }_d\sim N({\mu }_{\varphi },{\sigma }_{\varphi })$, ${\mu }_{\varphi }\sim N(0,2)$, ${\sigma }_{\varphi }\sim N^{+}(0,2)$; ${\lambda }_d$ likewise; thresholds $t_{rd}\sim {\text{logit}}^{-1}(N({\mu }_{t_r},{\sigma }_{t_r}))$ with ${\mu }_{t_r}\sim N(0,2)$, ${\sigma }_{t_r}\sim N^{+}(0,2)$ per race. $N^{+}$ is half-normal. Two SDs of $N(0,2)$ on ${\varphi }_r$ span roughly a 2%-98% chance of contraband; the same range holds for thresholds; $\lambda$ can differ across departments by a factor of about 50. Results hold with $N(0,4)$ priors (slower).

Sufficient statistics and exact likelihood (Eq. 3)

Aggregate $S_{rd}={\sum }_{i\in T_{rd}}{S}_i$, $H_{rd}={\sum }_{i\in T_{rd}}{H}_i$, with $n_{rd}$ stops fixed. Then

$$S_{rd}\sim \text{Bin}(n_{rd},p_{rd}),p_{rd}=1-I_{t_{rd}}({\varphi }_{rd},{\lambda }_{rd});H_{rd}\sim \text{Bin}(S_{rd},q_{rd}),q_{rd}=\mathbb{E}[p\mid p\ge t_{rd}]={\varphi }_{rd}\frac{1-I_{t_{rd}}({\mu }_{rd},{\lambda }_{rd}+1)}{1-I_{t_{rd}}({\varphi }_{rd},{\lambda }_{rd})},$$

with $I_x(\varphi ,\lambda )$ the beta CDF in mean/count parametrisation and ${\mu }_{rd}=({\varphi }_{rd}{\lambda }_{rd}+1)/({\lambda }_{rd}+1)$. The likelihood thus has $2DR$ binomial terms, independent of the number of stops.

Proof. $\mathbb{E}[p\mid p\ge t]=\frac{1}{1-F(t)}{\int }_t^{1}pf(p)dp$, and $p\text{Beta}(p;\alpha ,\beta )=\frac{\alpha }{\alpha +\beta }\text{Beta}(p;\alpha +1,\beta )$, giving $\varphi (1-F(t;\alpha +1,\beta ))/(1-F(t;\alpha ,\beta ))$, whose first argument pair in mean/count form is $(\mu ,\lambda +1)$. $\square$

Parameter counting (why the model is identified)

Per race-department pair there are three unknowns $(t,\varphi ,\lambda )$ but only two observed rates (search, hit), so $3DR$ unknowns face $2DR$ data. The additive structure reduces the total to $2(D+R-1)+DR$, which is $<2DR$ when $R\ge 3$ and $D\ge 5$. Hierarchical priors add regularisation. (Identification is by structure plus priors, not by data alone; see Simoiu 2017 - Limitations Extensions and Identification Ridge.)

Examples

Checking Eq. 3 numerically

With $\varphi =0.3,\lambda =25,t=0.2$ the formula gives $q=0.320578$; direct conditional expectation of a $\text{Beta}(7.5,17.5)$ above 0.2 gives $0.320578$ as well. Interpretation: the closed form is exact, so no simulation of individual stops is needed.

from scipy import stats
phi, lam, t = 0.3, 25, 0.2
a, b = phi*lam, (1-phi)*lam
q = phi * (1 - stats.beta(a+1, b).cdf(t)) / (1 - stats.beta(a, b).cdf(t))

PyMC-style sketch (untested).

with pm.Model():
    phi_r = pm.Normal("phi_r", 0, 2, shape=R); lam_r = pm.Normal("lam_r", 0, 2, shape=R)
    mu_t = pm.Normal("mu_t", 0, 2, shape=R); sd_t = pm.HalfNormal("sd_t", 2, shape=R)
    t = pm.math.invlogit(pm.Normal("t_raw", mu_t[:,None], sd_t[:,None], shape=(R, D)))
    phi = pm.math.invlogit(phi_r[:,None] + phi_d[None,:]); lam = pm.math.exp(lam_r[:,None] + lam_d[None,:])
    # beta CDF via pt.betainc(a, b, x): p = 1 - betainc(phi*lam, (1-phi)*lam, t)
    pm.Binomial("S", n, p, observed=S); pm.Binomial("H", S, q, observed=H)

Connections

Hierarchical Models (shrinkage of $t_{rd}$ toward race means), Efficient MCMC (the paper samples with Stan NUTS, 5 chains of 2,500 warm-up + 2,500 draws, $\hat{R}<1.05$), a related fast variant by Pierson et al. (2017) is cited.

See Also

Simoiu 2017 - Infra-Marginality and the Threshold Test - Overview