Robust Overdispersed Binomial Model Implementation

The point of departure for our estimation method is the fact that if the overdispersed binomial model of (1) and (2) is correctly specified then given a consistent estimate $\hat{\beta}$ for $\beta$ and hence $\hat{\pi}_i = [1 + \exp(-x'_i \hat{\beta})]^{-1}$, residuals of the form

$$r^*_i = \dfrac{y_i - m_i\hat{\pi}_i} {\sqrt{m_i\hat{\pi}_i (1- \hat{\pi}_i)}}$$ (7)

are approximately normal.²⁴ Indeed, a good moment estimator for $\sigma^2$ may be defined in terms of $r^*_i$ (McCullagh and Nelder, 1989, 126-127, eqn 4.23).

The LQD estimator focuses on the $\binom{h_{k}}{2}$ order statistic of the set $\{ \vert r^*_i - r^*_j \vert : i< j \}$ of absolute differences, where $\{ \vert r^*_i - r^*_j \vert : i< j \}$ has $\binom{n}{2}$ elements and $h_{k} = \lceil (n+k)/2 \rceil$. Following Croux et al. (1994) we use the notation

$$Q^*_{n} = \{ \vert r^*_i - r^*_j \vert : i< j \} _{\binom{h_{k}}{2}:\binom{n}{2}}$$ (8)

to denote that order statistic. To implement LQD we choose estimates $\hat{\beta}$ to minimize $Q^*_{n}$. Let $\hat{\beta}_{\text{LQD*}}$ designate the estimated coefficient vector and let $\hat{Q}^*_{n}$ designate the corresponding minimized value of $Q^*_{n}$. The LQD scale estimate is

$$\hat{\sigma} = \hat{Q}^*_{n} \dfrac{1}{\sqrt{2} \Phi^{-1}(5/8)},$$ (9)

where $\Phi^{-1}$ is the quantile function for the standard normal distribution (Rousseeuw and Croux, 1993, 1277).

The $\tanh$ estimator is based on the function

$$\psi(x) = \begin{cases}x , & \text{for } 0 \leq \vert x\vert \leq... ... \leq \vert x\vert \leq c 0 , & \text{for } c \leq \vert x\vert \end{cases}$$

where $p$, $c$, $d$, $A$ and $B$ are constants satisfying $0<p<c$ and other conditions.²⁵ Given a scale estimate $\hat{\sigma}$ and a vector of trial estimates $\hat{\beta}$, we compute the standardized residuals $r_i$ of (4) and then weights

$$w_i = \begin{cases}\psi(\vert r_i\vert)/\vert r_i\vert , & \text{for } r_i \neq 0 1 , & \text{for } r_i = 0 . \end{cases}$$

Observation $i$ is weighted by $w_i$ in what is otherwise the usual iteratively reweighted least squares algorithm to estimate $\beta$. The given scale value remains unchanged but the weights are updated to match the current coefficient estimates. Numerical convergence is required for both the $\hat{\beta}$ values and the weights. Because redescending $M$-estimators such as the $\tanh$ estimator have multiple solutions, starting values affect the results. We use $\hat{\beta}_{\text{LQD*}}$ to start the coefficients and use the LQD values $(r^*_i -\operatornamewithlimits{med}_i r^*_i)/\hat{\sigma}$ for an initial set of residuals ( $\operatornamewithlimits{med}_i r^*_i$ denotes the median of the $r^*_i$ values, $i=1,\dots,n$). To estimate the asymptotic variance of the coefficient estimates we use the sandwich estimator (White, 1994, 92) avar$(\hat{\beta}) = \hat{J}^{-1} \hat{I} \hat{J}^{-1}$ where $\hat{I}$ denotes the outer product of the score and $\hat{J}$ denotes the Hessian matrix.

The particular $\tanh$ estimator we use has $c=4.0$, $d=5.0$ and values for $p$, $A$ and $B$ as given in the indicated row of Table 2 in Hampel et al. (1981, 645).²⁶ The value of $c$ fixes the threshold for the magnitude of $r_i$ beyond which an observation is completely rejected by assigning it a weight $w_i=0$.