Unequal Cluster Sizes And Binary PQL

Why cluster size is not the whole story

In a binary generalized linear mixed model, clusters may contain different numbers of observations. It is tempting to assume that a cluster twice as large necessarily supplies twice as much useful information to a penalized quasi-likelihood (PQL) fit. That conclusion is generally too simple. For binary data, working weights depend on fitted probabilities, and a random intercept induces dependence among observations in a cluster. Consequently, both the outcome probabilities and the random-effects adjustment matter alongside the number of observations.

The issue is consequential because PQL estimates are constructed from a sequence of working linear mixed models. Unequal sizes can therefore alter the relative influence of clusters at each iteration. They may affect fixed-effect estimates, the amount of shrinkage applied to cluster effects, and the working evidence about random-effect variation. These effects do not have a universal direction: the result depends on covariate patterns, fitted probabilities, and the allocation of large and small clusters across groups.

Definition: Consider a binary random-intercept model with $Y_{ij}\mid b_j\sim\operatorname{Bernoulli}(\mu_{ij})$, link $g(\mu_{ij})=x_{ij}^{T}\beta+b_j$, and $b_j\sim N(0,\tau^2)$, independently across clusters. Here $j$ indexes clusters, $i=1,\ldots,n_j$ indexes observations within cluster $j$, and $n_j$ is its size. For binary logit PQL, the working weight at an iteration is $w_{ij}=\mu_{ij}(1-\mu_{ij})$. The sum $S_j=\sum_i w_{ij}$ is the cluster's working information about a common random-intercept displacement before accounting for its Gaussian penalty.

This definition assumes the stated conditional Bernoulli model, the specified link, normally distributed independent random intercepts, and a current set of fitted means used to form the PQL working quantities. It does not say that PQL is exact likelihood inference. As emphasized in the earlier material on binary PQL bias, the working response and weights approximate the nonlinear binary model; their values and the resulting estimates can be biased, especially when the approximation is poor.

Information in a cluster after random-intercept adjustment

For a given iteration, the working-response model has conditional variance $1/w_{ij}$. In a cluster with a common random intercept, the working covariance matrix is $V_j=W_j^{-1}+\tau^2\mathbf{1}\mathbf{1}^{T}$, where $W_j$ is diagonal with entries $w_{ij}$. Thus a cluster's contribution to information about a fixed effect that is constant across its observations, such as a cluster-level predictor, is not simply $S_j$.

Result: In this working model, the information in cluster $j$ for a common cluster-level mean is $$ I_j=\mathbf{1}^{T}V_j^{-1}\mathbf{1} =\frac{S_j}{1+\tau^2 S_j}. $$ The corresponding working variance of the cluster's weighted mean is $1/S_j+\tau^2$.

The formula is specific to a common cluster-level mean under a random-intercept working covariance. Fixed effects that vary within clusters have different information contributions, and general designs combine contributions through the full design matrix. Even so, this result gives a useful explanation for why increasing $n_j$ does not always increase effective information proportionally.

Why the information formula saturates

Let $a_j=W_j\mathbf{1}$, so that $a_j$ contains the working weights and $S_j=\mathbf{1}^{T}a_j$. Direct multiplication gives $$ V_j a_j=(W_j^{-1}+\tau^2\mathbf{1}\mathbf{1}^{T})W_j\mathbf{1} =\mathbf{1}+\tau^2 S_j\mathbf{1} =(1+\tau^2S_j)\mathbf{1}. $$ Therefore $V_j^{-1}\mathbf{1}=a_j/(1+\tau^2S_j)$. Multiplying on the left by $\mathbf{1}^{T}$ proves $I_j=S_j/(1+\tau^2S_j)$. The diagonal working variances alone would give information $S_j$; the shared random intercept reduces it.

When $\tau^2>0$, $I_j$ increases with $S_j$ but approaches $1/\tau^2$ as $S_j$ grows. Additional observations still help, but their marginal contribution to information about a common cluster-level mean diminishes. Moreover, $S_j$ itself depends on the fitted probabilities: for the logit link, weights are largest at $\mu_{ij}=0.5$ and smaller near 0 or 1. A large cluster with probabilities near the extremes can therefore have less working information than its raw size suggests.

Worked Example: Comparing Working Information Across Clusters

Suppose a binary logit PQL iteration has $\tau^2=0.5$. For illustration, consider a cluster of four observations with fitted probabilities all equal to 0.5 and another of eight observations with fitted probabilities all equal to 0.1. These are hypothetical fitted values, not observations from a real study.

For the four-observation cluster, each weight is $0.5(1-0.5)=0.25$, so $S_1=4(0.25)=1.00$. Its adjusted information is $$ I_1=\frac{1.00}{1+0.5(1.00)}=\frac{1}{1.5}=0.667. $$

For the eight-observation cluster, each weight is $0.1(1-0.1)=0.09$, so $S_2=8(0.09)=0.72$. Its adjusted information is $$ I_2=\frac{0.72}{1+0.5(0.72)}=\frac{0.72}{1.36}=0.529. $$

Despite having twice as many observations, the second cluster has less working information about a common cluster-level mean in this iteration. Both its lower weights and the random-intercept adjustment matter. This comparison does not establish which cluster contributes more to every fixed effect: a within-cluster covariate, for example, has a different information pattern.

How imbalance can move a fixed-effect estimate

For a cluster-level predictor, each cluster's working mean contributes according to its adjusted precision $I_j$, not merely as one equally weighted cluster and not as $n_j$ independent observations. If clusters with different predictor values also have different fitted probabilities or sizes, this weighting can change the fixed-effect estimate. PQL adds another layer: fitted probabilities and working weights are recalculated as the estimates change.

Worked Example: A Cluster-Level Contrast in a Working Fit

Consider four hypothetical clusters in one PQL working iteration. A binary predictor $x_j$ is constant within each cluster. Suppose $\tau^2=0.5$, every observation has working weight 0.25, and the weighted working-response mean for each cluster is as shown. With these weights, $S_j=0.25n_j$. The values below are illustrative working quantities, not actual binary outcomes.

Cluster	$x_j$	$n_j$	Working-response mean
1	0	4	0.0
2	0	8	1.0
3	1	2	1.5
4	1	6	2.0

The adjusted precisions are $I_j=S_j/(1+0.5S_j)$. For the four clusters these are $0.667$, $1.000$, $0.400$, and $0.857$, respectively. In this working comparison, the fitted group means are the precision-weighted averages: $$ \hat{\gamma}_0=\frac{0.667(0.0)+1.000(1.0)}{0.667+1.000}=0.600, \qquad \hat{\gamma}_1=\frac{0.400(1.5)+0.857(2.0)}{0.400+0.857}=1.841. $$

The working estimate of the predictor contrast is $\hat{\beta}_x=\hat{\gamma}_1-\hat{\gamma}_0=1.841-0.600=1.241$ working-response units. If instead each cluster mean received equal weight, the contrast would be $(1.5+2.0)/2-(0.0+1.0)/2=1.25$. The numerical difference is small here, but the calculation shows how cluster imbalance and random-intercept precision enter the working estimate. It is not a final log-odds estimate or an effect on the probability scale; the PQL algorithm updates the working quantities and coefficients.

Variance components and unequal cluster sizes

Cluster-specific random intercepts are also estimated with different precision. In the working Gaussian representation, let $q_j=\sum_i w_{ij}(z_{ij}-x_{ij}^{T}\beta)$, where $z_{ij}$ is the PQL working response. Conditional on the current fixed effects and variance, the penalized working estimate of the random intercept is $$ \hat b_j=\frac{\tau^2 q_j}{1+\tau^2S_j}. $$ The denominator reflects both the data's working information and the Gaussian penalty. For a given $q_j$, a larger $S_j$ means less shrinkage toward zero. But $q_j$, as well as $S_j$, depends on the working responses, covariates, and fitted means.

Worked Example: Different Shrinkage With Different Cluster Information

Suppose $\tau^2=0.5$ and two hypothetical clusters have the same working residual sum $q_j=0.40$, but different total weights: $S_1=0.50$ and $S_2=2.00$. The first cluster's working random-intercept estimate is $$ \hat b_1=\frac{0.5(0.40)}{1+0.5(0.50)} =\frac{0.20}{1.25}=0.160. $$ The second is $$ \hat b_2=\frac{0.5(0.40)}{1+0.5(2.00)} =\frac{0.20}{2.00}=0.100. $$

The second cluster has greater working information $S_j$, so the same residual sum corresponds to a smaller estimated random intercept under this formula. This is not evidence that larger clusters always have smaller random effects: their $q_j$ values need not be equal. Nor does one cluster's shrinkage calculation determine the fitted variance component $\tau^2$, which is estimated using information across the data and remains subject to PQL approximation.

Diagnostics, alternatives, and limitations

When cluster sizes are highly unequal, inspect the distribution of $n_j$ alongside fitted probabilities and working weights. In particular, check whether large clusters are concentrated in a predictor group, whether many fitted means are near zero or one, and whether a small number of clusters dominate the weighted equations. Track whether estimates or working weights change substantially over iterations, and assess sensitivity to influential clusters. These checks flag instability or uneven influence; they do not correct PQL approximation bias.

A common mistake is to reason, “This cluster has twice as many observations, so it must count twice as much in the fixed-effect fit.” That ignores probability-dependent weights and the shared random intercept. The correction is to distinguish raw size $n_j$, working information $S_j$, and adjusted information such as $S_j/(1+\tau^2S_j)$ for a common cluster-level mean. For other covariate patterns, inspect the corresponding full working information matrix rather than applying this scalar formula indiscriminately.

Equal weighting of cluster summaries answers a different descriptive question and is not automatically a better PQL estimator. Conversely, treating every observation as independent discards the modeled within-cluster dependence. PQL's computational attraction is that it repeatedly fits a working mixed model, but its estimates inherit the approximation's limitations. Likelihood-based integration is a distinct alternative for binary mixed models; it does not make unequal cluster sizes irrelevant, and its details are beyond this tutorial's focus.

Unequal sizes also do not imply a universal upward or downward shift in either fixed effects or $\tau^2$. The direction depends on how size relates to predictors, outcomes, fitted probabilities, and random-effect patterns. The working information calculations clarify possible mechanisms, but they are not a general bias correction or a guarantee of valid finite-sample inference.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For a logit PQL fit with $\mu=0.5$, what is the working weight for one observation? What is the total $S_j$ for six such observations?
With $S_j=1$ and $\tau^2=0.25$, calculate the adjusted information for a common cluster-level mean.
Why can a larger cluster have less working information than a smaller one?
In the shrinkage formula, what happens to $\hat b_j$ as $S_j$ increases while $q_j$ and $\tau^2$ stay fixed and positive?
Does the formula $S_j/(1+\tau^2S_j)$ describe every fixed effect in a mixed model? Explain.

Solutions. (1) The weight is $0.5(1-0.5)=0.25$, so $S_j=6(0.25)=1.50$. (2) $I_j=1/(1+0.25)=0.800$. (3) Its fitted probabilities may be closer to zero or one, reducing the weights and therefore $S_j$; random-intercept dependence further reduces information about a shared cluster-level mean. (4) The denominator $1+\tau^2S_j$ increases, so the positive estimate $\tau^2q_j/(1+\tau^2S_j)$ decreases toward zero. (5) No. It applies to a common cluster-level mean under the stated random-intercept working covariance; within-cluster covariates and general fixed-effect designs require their corresponding matrix information.