Bias Reduction for Grouped Binomial Responses

Grouped Binomial Responses in a Logistic Model

A grouped binomial response records a number of successes out of a known number of trials at each covariate pattern. For example, a row might report 2 infections among 40 patients with the same exposure status. Treating that row as one binary response would discard the trial count; treating it as 40 separate records is possible, but often unnecessary. The grouped representation retains the information needed for logistic regression and Firth’s adjustment, provided the trial count is included correctly.

Let $Y_i$ be the number of successes among $m_i$ trials in row $i$, and let $x_i$ be its $p$-dimensional covariate vector, including an intercept if required. Assume the rows are independent, $m_i$ is fixed and positive, and the trials within row $i$ are independent Bernoulli trials with a common success probability $p_i$. In the logistic model, $p_i=\{1+\exp(-x_i^\mathsf{T}\beta)\}^{-1}$. Assume the design matrix has full column rank so the Fisher information is nonsingular at finite $\beta$.

Definition: The grouped-binomial logistic model has $Y_i\sim\operatorname{Binomial}(m_i,p_i)$, independently across rows, with $\operatorname{logit}(p_i)=x_i^\mathsf{T}\beta$. Write $q_i=1-p_i$. Its grouped-row Fisher weight is $w_i=m_i p_iq_i$, and its Fisher information is $I(\beta)=\sum_i w_i x_ix_i^\mathsf{T}$.

The ordinary grouped log-likelihood, apart from its displayed binomial coefficient, is

\ell(\beta)=\sum_i\left\{\log {m_i\choose Y_i}+Y_i\log p_i+(m_i-Y_i)\log q_i\right\}.

The binomial coefficient does not depend on $\beta$. The ordinary score is therefore $\sum_i x_i(Y_i-m_ip_i)$. Firth’s method adds one-half the log determinant of the Fisher information to the log-likelihood. The adjustment is not made by adding an artificial success to every row: its size depends on the model’s leverage.

Deriving the Grouped Adjusted Score

Theorem: For a full-rank grouped-binomial logistic model with finite $\beta$, the score of the Jeffreys-penalized log-likelihood $\ell^*(\beta)=\ell(\beta)+\tfrac12\log|I(\beta)|$ has the form $$ U^*(\beta)=\sum_i x_i\left\{Y_i-m_ip_i+h_i\left(\tfrac12-p_i\right)\right\}, \qquad h_i=w_i x_i^\mathsf{T}I(\beta)^{-1}x_i. $$ Here $h_i$ is the leverage of grouped row $i$, using its full trial-count weight $w_i=m_ip_iq_i$.

To derive the result, differentiate the ordinary likelihood. Since $\partial p_i/\partial\beta=p_iq_ix_i$, its score is $\sum_i x_i(Y_i-m_ip_i)$, and the negative expected derivative of that score is $I(\beta)=\sum_i m_ip_iq_ix_ix_i^\mathsf{T}$. For the penalty, use the matrix identity $\partial\log|I|/\partial\beta_j=\operatorname{tr}\{I^{-1}(\partial I/\partial\beta_j)\}$, valid when $I$ is differentiable and nonsingular. Because $\partial(p_iq_i)/\partial\beta_j=p_iq_i(1-2p_i)x_{ij}$, we obtain

\frac12\frac{\partial\log|I|}{\partial\beta_j} =\frac12\sum_i m_ip_iq_i(1-2p_i)x_{ij}\,x_i^\mathsf{T}I^{-1}x_i =\sum_i h_i\left(\tfrac12-p_i\right)x_{ij}.

Adding these penalty derivatives to the ordinary score gives the stated adjusted score. Thus the adjustment uses both the row’s leverage and how far its fitted probability is from one-half. It vanishes for a row at $p_i=1/2$, but the full adjusted score need not vanish there because other rows may contribute.

Grouped Rows and Individual Trials

A grouped row represents $m_i$ Bernoulli trials sharing the same covariates. If the data are expanded into those trials, each trial in row $i$ has weight $p_iq_i$, while the information matrix remains the same: summing the $m_i$ individual contributions gives $m_ip_iq_ix_ix_i^\mathsf{T}$. This gives a precise relationship between the two leverage scales.

Result: In the expanded trial-level representation, each trial in grouped row $i$ has leverage $h_i^{\mathrm{trial}}=p_iq_i x_i^\mathsf{T}I^{-1}x_i$. Consequently, $h_i=m_i h_i^{\mathrm{trial}}$. The grouped-row leverage is the sum of the leverages of the $m_i$ trials it represents; it is not the leverage of one Bernoulli observation.

To verify this, the weighted design matrix for the expanded data has one row $x_i^\mathsf{T}$ for each trial in group $i$, each with weight $p_iq_i$. The diagonal entry of the weighted hat matrix for any one such trial is its weight times $x_i^\mathsf{T}I^{-1}x_i$. Summing those identical entries over $m_i$ trials gives $m_ip_iq_i x_i^\mathsf{T}I^{-1}x_i=h_i$. This also explains why grouped-row leverage can exceed one without contradiction: it aggregates multiple trial-level leverages. The total leverage across all expanded trials is $\operatorname{tr}(I^{-1}I)=p$, the number of model coefficients.

Worked Example: One Group With No Successes

Suppose 10 trials share only an intercept, and none succeeds: $Y=0$, $m=10$. Write $p=\operatorname{logit}^{-1}(\beta)$. The ordinary score is $Y-mp=-10p$. With an intercept-only model, $I=mpq$, so the sole row’s leverage is $h=mpq(1/I)=1$. The adjusted score equation is

-10p+\left(\tfrac12-p\right)=0 \quad\Longrightarrow\quad 11p=\tfrac12 \quad\Longrightarrow\quad \hat p_F=\frac1{22}.

Thus $\hat\beta_F=\log\{(1/22)/(21/22)\}=-\log 21$, a finite estimate. More generally, the intercept-only adjusted equation is $Y-mp+\tfrac12-p=0$, giving $\hat p_F=(Y+\tfrac12)/(m+1)$. Here the grouped row has leverage 1, while each of its 10 trials has leverage $1/10$. This example shows why the trial count belongs in the information even though there is only one grouped row.

A Sparse Saturated Example

In a saturated two-group logistic model, each group has its own fitted probability. Each grouped row has leverage 1: with two rows and two group-specific coefficients, the model fits one parameter per row. The adjusted score for each group therefore reduces to $Y_i-m_ip_i+\tfrac12-p_i=0$.

Worked Example: Sparse Treatment and Control Counts

Suppose the control group has 0 successes in 4 trials and the treatment group has 1 success in 3 trials. Let the model be $\operatorname{logit}(p_i)=\beta_0+\beta_1 z_i$, where $z_i=0$ for control and $z_i=1$ for treatment. This two-parameter model is saturated for the two observed groups.

Control probability. The control adjusted score equation is $0-4p_C+\tfrac12-p_C=0$. Hence $p_C=\tfrac1{10}$, with control odds $p_C/(1-p_C)=1/9$.

Treatment probability. The treatment equation is $1-3p_T+\tfrac12-p_T=0$. Hence $p_T=\tfrac38$, with treatment odds $3/5$.

Coefficients. The intercept is the control log-odds, $\hat\beta_{0,F}=\log(1/9)=-\log 9$. The treatment coefficient is the difference in log-odds: $\hat\beta_{1,F}=\log(3/5)-\log(1/9)=\log(27/5)$. The fitted probabilities are strictly between zero and one despite a zero-success group. Each grouped row has leverage 1, but the individual trial leverages are $1/4$ for each control trial and $1/3$ for each treatment trial. Summing within each group recovers row leverage 1.

Checking an Adjusted Score With a Numerical Step

For a model that is not saturated, the row leverage depends on the entire design and fitted probabilities. A single evaluation of the adjusted score is still useful: it checks that the trial counts enter both the information and the leverage correctly. The following calculation is one modified IRLS-style scoring step from a specified starting value, not a claim that the step is already the final estimate.

Worked Example: Three Grouped Rows at the Starting Value

Take three rows with $z_i=-1,0,1$, covariates $x_i=(1,z_i)^\mathsf{T}$, trial counts $m_i=(2,4,2)$, and successes $Y_i=(0,2,2)$. Start at $\beta=(0,0)^\mathsf{T}$, so every fitted probability is $1/2$. The grouped weights are $w_i=m_i/4=(1/2,1,1/2)$, and

I=X^\mathsf{T}WX= \begin{pmatrix}2&0\\0&1\end{pmatrix}, \qquad I^{-1}= \begin{pmatrix}1/2&0\\0&1\end{pmatrix}.

For $z=-1$ and $z=1$, $x_i^\mathsf{T}I^{-1}x_i=3/2$, so each row has leverage $h_i=(1/2)(3/2)=3/4$. For $z=0$, the quadratic form is $1/2$ and the leverage is $1(1/2)=1/2$. These sum to 2, the number of coefficients. The corresponding trial-level leverages are $3/8,1/8,3/8$, respectively, after dividing each row’s leverage by its trial count.

At this starting value, every Firth adjustment $h_i(1/2-p_i)$ is zero. The ordinary residuals $Y_i-m_ip_i$ are $(-1,0,1)$, so

U^*= \sum_i x_i(Y_i-m_ip_i) = \begin{pmatrix}0\\2\end{pmatrix}, \qquad I^{-1}U^*= \begin{pmatrix}0\\2\end{pmatrix}.

The scoring step therefore leaves the intercept at 0 and increases the slope by 2. The check is internally consistent: the intercept score sums the residuals to zero, while the slope score weights them by $z_i$, giving 2. Away from this starting value, the leverage adjustment generally contributes to both score components.

Why Grouping Matters and a Common Pitfall

Grouping is convenient when many trials share the same covariates, but the grouping must preserve the model information. Under the stated independent-trial model, replacing $m_i$ identical-covariate Bernoulli observations by their total $Y_i$ gives the grouped binomial likelihood. The expanded likelihood and grouped likelihood can differ by factors that do not depend on $\beta$; the score and Fisher information agree. Accordingly, the Firth adjustment agrees too when the grouped weight $m_ip_iq_i$ is used.

A common mistake is to calculate $w_i=p_iq_i$ for a grouped row as if the entire row were one Bernoulli observation. That omits $m_i$, changes the information matrix, and produces an incorrect leverage and adjustment. For example, in the 10-trial intercept-only example, this mistake would still yield leverage 1 for its only row, but it would incorrectly use information $pq$ instead of $10pq$ in the model’s information and scoring calculations. The apparent agreement in row leverage is not evidence that the trial count can be ignored: the score’s ordinary term, information, uncertainty, and trial-level leverage all depend on it.

The other frequent error is to interpret $h_i$ as the leverage of one trial. It is the aggregate for the entire grouped row. For row-level diagnostics that compare groups, $h_i$ is the appropriate quantity in the grouped adjusted score. For a comparison with an expanded Bernoulli data set, use $h_i/m_i$ for each trial in that row.

Key takeaway: In grouped-binomial Firth logistic regression, use $m_ip_iq_i$ in the Fisher information and $h_i=m_ip_iq_i\,x_i^\mathsf{T}I^{-1}x_i$ in the adjusted score. This grouped-row leverage is the sum of the $m_i$ equal trial-level leverages, not the leverage of one trial.

Check Your Understanding

Use the grouped model and leverage definitions to answer these questions.

For a grouped row with $m_i=8$, what is the relationship between its row leverage and the leverage of each represented trial?
In an intercept-only model with 6 trials and 2 successes, calculate the Firth-adjusted fitted probability.
Why does the grouped information contribution contain $m_i$, even though the row contributes only one observed count?
In the three-row example, why do the grouped-row leverages sum to 2?
What error results from treating a grouped count as one Bernoulli response when calculating the information matrix?