Robit Regression

Definition

Robit regression (Liu 2006) replaces the normal cumulative distribution function (CDF) link in probit with a Student-t CDF: $\Pr(y_i=1\mid x_i,\beta)=F_\nu(x_i'\beta)$ , where $F_\nu$ is the t-CDF with $\nu$ degrees of freedom, location zero, and scale one. As $\nu\to\infty$ , robit recovers probit exactly; robit(7) closely approximates the logistic link (max quantile-quantile (Q-Q) deviation 0.0006 with scale $\sigma=1.5484$ ). The key advantage over probit and logistic is bounded influence: the maximum likelihood estimator (MLE) influence function is bounded by $\nu$ in absolute value, making the MLE resistant to outlying observations.

Key Ideas

Robustness of MLE: The score factor $f_\nu(x'\hat\beta)/F_\nu(x'\hat\beta)\cdot x'\hat\beta$ is bounded (approaches $-\nu$ as $x'\hat\beta\to-\infty$ , 0 as $x'\hat\beta\to+\infty$ ). For logistic, influence is unbounded; for probit, it diverges as $(x'\hat\beta)^2$ . Bounded influence = no single observation can dominate the MLE.
Scale mixture structure: The complete-data model $\tau_i|\theta\sim\text{Gamma}(\nu/2,\nu/2)$ , $z_i|(\tau_i,\theta)\sim N(x_i'\beta,1/\tau_i)$ , $y_i=\mathbf{1}[z_i>0]$ marginalizes to the t-link. The Gamma weights $\tau_i$ downweight observations whose latent $z_i$ is far from the boundary, enabling outlier identification: small $\hat\tau_i$ flags observation $i$ as influential.
Robit(7) $\approx$ logistic: Scale $\sigma=1.5484$ is chosen to minimize the maximum Q-Q distance between $F_7(x/\sigma)$ and $F_\text{logistic}(x)$ over the 0.001–0.999 range; the max distance is 0.0006. Robit with 5 or 6 df balances robustness and logistic approximation for practitioners who want both.
Expectation-Maximization (EM)-type MLE algorithms:
- EM (known $\nu$ ): E-step computes $\hat\tau_i$ and $\hat z_i$ analytically as ratios of t-CDF tail probabilities; M-step is weighted least squares (WLS) $\hat\beta=(\sum\hat\tau_i x_ix_i')^{-1}(\sum\hat\tau_i x_i\hat z_i)$ .
- ECME (Expectation/Conditional Maximization Either; unknown $\nu$ ): CM-step 1 = WLS for $\beta$ ; CM-step 2 = profile likelihood search over $\nu$ by half-interval method on eq. (11).
- PX-EM (Parameter-Expanded EM; Liu-Rubin-Wu 1998): Expand model with working parameters $(\alpha,\sigma)$ ; compute $\hat\alpha=\bar\tau$ and $\hat\sigma^2=n^{-1}(S_{\tau zz}-S_{\tau xz}'S_{\tau xx}^{-1}S_{\tau xz})$ ; update $\beta^{(t+1)}=(\hat\alpha/\hat\sigma)\hat\beta^*$ . Converges dramatically faster than EM.
Bayesian data augmentation (DA) algorithms:
- Basic DA: Prior $\beta\sim t_p(0,S_0^{-1},\nu_0)$ hierarchically represented via $\tau_0\sim\text{Gamma}(\nu_0/2,\nu_0/2)$ . I-step: draw $z_i$ from truncated $t(\mu_i,1,\nu)$ (by inversion), $\tau_i$ from $\text{Gamma}((\nu+1)/2,\;(\nu+(z_i-\mu_i)^2)/2)$ , $\tau_0$ from its prior. P-step: draw $\beta\sim N(\hat\beta,(\tau_0S_0+S_{\tau xx})^{-1})$ .
- E-DA 1: After P-step, rescale $\hat\beta$ by $\chi_n/s^{1/2}$ where $s^2=\sum\tau_i(z_i-x_i'\hat\beta)^2+\hat\beta'\tau_0S_0\hat\beta$ . Integrates out the common scale of $z_i$ ; converges faster than DA.
- E-DA 2: Further rescale by $(w/\chi^2_{\nu_0+n\nu})^{1/2}$ where $w=\sum_{i=0}^n\nu_i\tau_i$ . Integrates out both the $z$ -scale and the $\tau$ -scale. Fastest convergence.
Outlier detection: Fitted weights $\hat\tau_i\ll1$ identify influential observations. In Finney (1947) vaso-constriction data ( $n=39$ ), $\hat\nu\approx0.11$ (near-Cauchy, quasi-separation), and observations 4, 18, 24 are flagged — consistent with Pregibon (1982) resistant fits.

How It Works

EM E-Step (Known $\nu$ )

$\hat\tau_i = \frac{\nu+1}{\nu}\cdot\frac{\Pr(t_{\nu+2}<-\sqrt{1+2/\nu}\,\mu_i)}{\Pr(t_\nu<-\mu_i)} \quad (y_i=0)$ with sign flip for $y_i=1$ . The pseudo-response is: $\hat z_i = \mu_i + (2y_i-1)\frac{f_\nu(\mu_i)}{y_i-(2y_i-1)\Pr(t_{\nu+2}<-\sqrt{1+2/\nu}\,\mu_i)}$ Then M-step: $\hat\beta=\bigl(\sum\hat\tau_i x_ix_i'\bigr)^{-1}\bigl(\sum\hat\tau_i x_i\hat z_i\bigr)$ .

E-DA 2 Cycle (One Bayesian Iteration)

Draw $z_i\mid\beta,y_i$ from truncated $t(\mu_i,1,\nu)$ by inversion of $F_\nu$ .
Draw $\tau_i\sim\text{Gamma}\bigl((\nu+1)/2,\;(\nu+(z_i-\mu_i)^2)/2\bigr)$ .
Draw $\tau_0\sim\text{Gamma}(\nu_0/2,\nu_0/2)$ .
Compute $\hat\beta=(\tau_0S_0+S_{\tau xx})^{-1}S_{\tau xz}$ ; draw $w\sim\text{Gamma}((\nu_0+n\nu)/2,1/2)$ , $s^2\sim\text{Gamma}(n/2,w/2)$ .
Set $\beta^{\text{new}}=sw^{-1/2}\hat\beta$ (E-DA 2 rescaling).

Why It Matters

Provides a tractable robust alternative to logistic and probit: t-link with small $\nu$ downweights outliers automatically via WLS weights, with no user-set tuning.
PX-EM achieves much faster MLE convergence than basic EM, extending the Liu-Rubin-Wu (1998) framework to binary regression.
Bayesian DA is simpler than Holmes-Held (2006) exact logistic sampler (standard Gamma and truncated-t draws only, no Generalized Inverse Gaussian (GIG) rejection), making it easier to extend to multivariate binary responses.
Fitted weights $\hat\tau_i$ are a natural outlier-detection diagnostic, produced as a free by-product of the EM E-step.

Open Questions

Robit vs. exact logistic: Robit(7) approximates the logistic to within 0.0006 but is not exact. For strict Bayesian logistic inference, Holmes-Held (2006) Kolmogorov-Smirnov (KS) scale mixture is preferable.
Degrees of freedom estimation: ECME profile likelihood search for $\nu$ can be noisy in small samples; placing a prior on $\nu$ and marginalizing is preferable in Bayesian settings.
Polychotomous extension: Liu (2006) mentions the multivariate binary extension briefly; a complete algorithm comparable to Holmes-Held's polychotomous Gibbs is not provided.