Worked Quasi-Poisson Analysis With Exposure

Why Exposure Matters

Counts are only comparable as measures of event frequency when the amounts of time, people, or other opportunity for events are comparable. Ten events over six patient-years and seven events over three patient-years do not represent rates of $10$ and $7$ events per patient-year: each count must be related to its exposure. A quasi-Poisson model does this while allowing the variance to be a multiple of the mean.

We will compare event rates in a reference group and a programme group. The observations have known positive exposures $E_i$, counts $Y_i$, and group indicator $x_i$, where $x_i=0$ denotes the reference group and $x_i=1$ the programme group. Assume observations are independent, and specify the mean and variance by

E[Y_i]=\mu_i,\qquad \operatorname{Var}(Y_i)=\phi\mu_i,\qquad \log(\mu_i)=\log(E_i)+\beta_0+\beta_1x_i.

Here $\phi>0$ is the dispersion parameter. The term $\log(E_i)$ is an offset: its coefficient is fixed at one, so it adjusts the mean for exposure rather than introducing another coefficient to estimate. Equivalently, the model says that the expected count is exposure multiplied by a group-specific rate.

Model and Rate Interpretation

Definition: In this two-group log-link quasi-Poisson model, the fitted rate for group $x$ is $\lambda_x=\exp(\beta_0+\beta_1x)$, in events per unit of exposure. Thus $\lambda_0=\exp(\beta_0)$, $\lambda_1=\exp(\beta_0+\beta_1)$, and the programme-to-reference rate ratio is $\exp(\beta_1)$. The quasi-Poisson variance specification is $\operatorname{Var}(Y_i)=\phi\mu_i$; it specifies the first two moments, not a full probability distribution.

The coefficient $\beta_1$ is a log rate ratio, not a difference in event counts and not the programme rate itself. The fitted expected count for an observation with exposure $E_i$ in group $x_i$ is $E_i\hat\lambda_{x_i}$. In this tutorial, rate estimates come from the quasi-score equations, and uncertainty accounts for estimated dispersion.

Worked Example: Fitting Rates and a Rate Ratio

Suppose the observed counts and exposures are as follows. Exposure is measured in patient-years.

Group	Exposure $E_i$	Count $Y_i$
Reference	1	0
Reference	2	5
Reference	1	0
Reference	2	5
Programme	0.5	0
Programme	0.5	0
Programme	1	2
Programme	1	5

The reference group has total count $10$ and total exposure $6$, so its fitted rate is $10/6=5/3\approx1.667$ events per patient-year. The programme group has total count $7$ and total exposure $3$, so its fitted rate is $7/3\approx2.333$ events per patient-year. The estimated rate ratio is

\frac{\hat\lambda_1}{\hat\lambda_0} =\frac{7/3}{5/3}=\frac75=1.4.

Thus the programme group’s fitted event rate is $1.4$ times the reference rate, or $40\%$ higher. The fitted coefficients in reference-indicator coding are $\hat\beta_0=\log(5/3)\approx0.5108$ and $\hat\beta_1=\log(7/5)\approx0.3365$. The fitted programme rate is $\exp(\hat\beta_0+\hat\beta_1)=7/3$, not $\hat\beta_1$.

Why the Fitted Rates Are Count Totals Divided by Exposure Totals

For this model, the group-specific fitted rates have a simple form. The calculation follows directly from the quasi-score for the log link. For a coefficient vector $\beta$, the score contribution for observation $i$, apart from the common factor $1/\phi$, is $x_i(Y_i-\mu_i)$, where $x_i$ here denotes its design row. With separate group means, the score equation for each group is the sum of its residuals being zero.

Result: If each of two groups has positive total exposure and its mean is modelled by its own log rate, the quasi-score fit satisfies $$ \hat\lambda_g=\frac{\sum_{i:\,x_i=g}Y_i}{\sum_{i:\,x_i=g}E_i}, \qquad g\in\{0,1\}. $$ This statement concerns the fitted mean parameters; it does not require a Poisson distribution.

Derivation. Write $\alpha_g=\log(\lambda_g)$ for the log rate in group $g$. For an observation in that group, $\mu_i=E_i\exp(\alpha_g)$ and $\partial\mu_i/\partial\alpha_g=\mu_i$. The quasi-score for $\alpha_g$, apart from the nonzero common factor $1/\phi$, is therefore

\sum_{i:\,x_i=g}\frac{\partial\mu_i}{\partial\alpha_g} \frac{Y_i-\mu_i}{\mu_i} = \sum_{i:\,x_i=g}(Y_i-\mu_i).

Setting this score to zero gives $\sum_{i:x_i=g}Y_i=\sum_{i:x_i=g}E_i\exp(\alpha_g)$. Since $\exp(\alpha_g)$ is constant within the group and total exposure is positive, division by that exposure gives the stated rate estimate. In the example, the two group score equations give $10=6\hat\lambda_0$ and $7=3\hat\lambda_1$, yielding the rates already calculated. $\square$

Estimating Dispersion With Pearson Residuals

The rate fit does not by itself specify how variable the counts are around their fitted means. A common estimate of dispersion is the Pearson statistic divided by its residual degrees of freedom. For a fitted mean $\hat\mu_i$, the Pearson residual is $(Y_i-\hat\mu_i)/\sqrt{\hat\mu_i}$, and its squared value is the observation’s contribution to the statistic. This estimate requires positive fitted means, as here.

Definition: For $n$ independent observations and $p$ estimated regression coefficients, the Pearson dispersion estimate is $$ \hat\phi=\frac{X_P^2}{n-p}, \qquad X_P^2=\sum_{i=1}^n\frac{(Y_i-\hat\mu_i)^2}{\hat\mu_i}, $$ provided $n>p$. It is an estimate of the variance multiplier, not a rate or a count.

Worked Example: Calculating Pearson Dispersion

For the fitted rates above, the reference fitted means are $5/3,10/3,5/3,10/3$, because each is exposure times $5/3$. The reference-group Pearson contributions are

\frac{(0-5/3)^2}{5/3}+ \frac{(5-10/3)^2}{10/3}+ \frac{(0-5/3)^2}{5/3}+ \frac{(5-10/3)^2}{10/3} =\frac53+\frac56+\frac53+\frac56=5.

The programme fitted means are $7/6,7/6,7/3,7/3$. The first two observations have count zero, so together they contribute $7/6+7/6=7/3$. The other two contributions are

\frac{(2-7/3)^2}{7/3}=\frac1{21}, \qquad \frac{(5-7/3)^2}{7/3}=\frac{64}{21}.

The programme contribution is $7/3+1/21+64/21=114/21=38/7$. There are $n=8$ observations and $p=2$ estimated coefficients, so the residual degrees of freedom are $8-2=6$. Consequently,

X_P^2=5+\frac{38}{7}=\frac{73}{7}, \qquad \hat\phi=\frac{73/7}{6}=\frac{73}{42}\approx1.738.

The estimated dispersion is above one, indicating more variability than the Poisson variance would specify at these fitted means. It is still an estimate from a small data set, not proof that a particular alternative count distribution generated the observations.

Uncertainty for a Fitted Rate

For a log-link quasi-Poisson fit, the estimated covariance matrix of the regression coefficients is the usual inverse weighted information matrix multiplied by $\hat\phi$. In this example, it can be understood especially directly by first using the two group log rates $\alpha_0=\log\lambda_0$ and $\alpha_1=\log\lambda_1$. With the fitted means, the information matrix for these parameters, before multiplying the covariance by dispersion, is diagonal: its $g$th diagonal entry is the sum of fitted means in group $g$. The score equations make this sum equal to the group’s observed total count $Y_g$. Thus the estimated variance of the fitted log rate in group $g$ is $\hat\phi/Y_g$.

To express the same result in the regression coefficients, $\alpha_0=\beta_0$ and $\alpha_1=\beta_0+\beta_1$. Therefore the fitted programme log rate is $\beta_0+\beta_1$, and its estimated variance is $\hat\phi/Y_1$. The log rate ratio is $\beta_1=\alpha_1-\alpha_0$; since the two group log rates have zero estimated covariance in this group-specific parameterization, its estimated variance is $\hat\phi(1/Y_0+1/Y_1)$.

Result: In a two-group model with separate log rates, independent observations, and the log-link quasi-Poisson variance function, the model-based estimated variances are $$ \widehat{\operatorname{Var}}(\hat\alpha_g)=\frac{\hat\phi}{Y_g}, \qquad \widehat{\operatorname{Var}}(\hat\beta_1) =\hat\phi\left(\frac1{Y_0}+\frac1{Y_1}\right), $$ when both group totals $Y_g$ are positive. These are estimated covariance results, not claims of exact finite-sample sampling distributions.

For completeness, the information calculation uses $\partial\mu_i/\partial\alpha_g=\mu_i$ and variance function $V(\mu_i)=\mu_i$. The information contribution for $\alpha_g$ is $(\partial\mu_i/\partial\alpha_g)^2/V(\mu_i)=\mu_i$. Contributions from different groups do not overlap, so the information is diagonal. Summing within group gives $\sum_i\hat\mu_i=Y_g$, by the score equation derived above. Inverting this diagonal matrix and multiplying by $\hat\phi$ gives $\hat\phi/Y_g$. For the log rate ratio, the contrast is $\alpha_1-\alpha_0$; its variance is the sum of the two diagonal variances. $\square$

Worked Example: An Approximate Interval for the Programme Rate

The estimated programme log rate is $\log(7/3)$, and its estimated standard error is

\operatorname{SE}\{\log(\hat\lambda_1)\} =\sqrt{\frac{\hat\phi}{Y_1}} =\sqrt{\frac{73/42}{7}} =\sqrt{\frac{73}{294}} \approx0.498296.

Using the approximate $t$ critical value $t_{0.975,6}=2.447$, the margin on the log-rate scale is $2.447\sqrt{73/294}\approx1.21933$. Keeping precision through the calculation gives

\log(7/3)\pm1.21933 \approx0.84730\pm1.21933 =[-0.37203,\ 2.06663].

The exponential function is increasing, so exponentiating both endpoints gives the approximate 95% interval for the programme rate:

[\exp(-0.37203),\exp(2.06663)] \approx[0.689,\ 7.898] \quad\text{events per patient-year}.

This interval is wide: with only seven programme events and an estimated dispersion based on eight observations, substantial uncertainty remains. The $t$-based interval is an approximate quasi-GLM procedure, not an exact interval guaranteed by a fully specified count distribution.

Common Pitfall and Practical Interpretation

A frequent mistake is to compare the event totals and conclude that the programme group has a lower rate because it has fewer events ($7$ rather than $10$). That comparison ignores exposure: the programme group has half the total exposure. Dividing by exposure gives $7/3\approx2.333$ events per patient-year, compared with $10/6\approx1.667$ in the reference group, and hence an estimated rate ratio of $1.4$. Another error is to call $\hat\beta_1\approx0.3365$ the rate ratio. It is the log rate ratio; exponentiating it gives the rate ratio $1.4$.

The fitted rates and their ratio address the size and direction of the group difference. The dispersion estimate adjusts the estimated covariance for variability beyond the Poisson variance. In reporting results, state the exposure unit, the estimated rates, the rate ratio, and whether an uncertainty interval is approximate. A quasi-Poisson fit provides moment-based inference; it does not, on its own, specify a complete distribution for future counts.

Check Your Understanding

Use the model and calculations in this tutorial to check both the interpretation and the arithmetic.

If a group records $12$ events over $8$ patient-years, what is its fitted rate under a group-specific log-link model?
In the example, why is the programme-to-reference rate ratio $7/5$, rather than $7/10$?
What does $\hat\phi=73/42$ describe, and how was its denominator determined?
Why does the estimated variance of the programme log rate use the programme event total $Y_1=7$?
Why must the endpoints of a log-rate interval be exponentiated before they can be reported as rates?