Improve narrative flow and math derivations in count regression chapter

_subfiles/count-regression/_exr-prac-glm-interp.qmd

@@ -14,11 +14,11 @@ where $x_i$ is a binary indicator ($x_i = 0$ or $x_i = 1$).
 
 **(a)** Express $\mu_i$ as a function of $x_i$.
 
-**(b)** Interpret $e^{\beta_0}$.
+**(b)** Interpret $\exp{\beta_0}$.
 
-**(c)** Interpret $e^{\beta_1}$.
+**(c)** Interpret $\exp{\beta_1}$.
 
-**(d)** If $\hat\beta_0 = 1.2$ and $\hat\beta_1 = 0.5$,
+**(d)** If $\eb_0 = 1.2$ and $\eb_1 = 0.5$,
 compute the estimated mean event count for $x_i = 0$ and $x_i = 1$.
 :::
 
@@ -27,25 +27,25 @@ compute the estimated mean event count for $x_i = 0$ and $x_i = 1$.
 **(a)**
 
 $$
-\mu_i = e^{\beta_0 + \beta_1 x_i} = e^{\beta_0} \cdot (e^{\beta_1})^{x_i}
+\mu_i = \exp{\beta_0 + \beta_1 x_i} = \exp{\beta_0} \cdot (\exp{\beta_1})^{x_i}
 $$
 
-For $x_i = 0$: $\mu_0 = e^{\beta_0}$.
-For $x_i = 1$: $\mu_1 = e^{\beta_0 + \beta_1}$.
+For $x_i = 0$: $\mu_0 = \exp{\beta_0}$.
+For $x_i = 1$: $\mu_1 = \exp{\beta_0 + \beta_1}$.
 
 **(b)**
 
-$e^{\beta_0}$ is the expected mean count when $x_i = 0$ (the reference group).
+$\exp{\beta_0}$ is the expected mean count when $x_i = 0$ (the reference group).
 
 **(c)**
 
 $$
-e^{\beta_1}
+\exp{\beta_1}
 = \frac{\mu_1}{\mu_0}
-= \frac{e^{\beta_0+\beta_1}}{e^{\beta_0}}
+= \frac{\exp{\beta_0+\beta_1}}{\exp{\beta_0}}
 $$
 
-$e^{\beta_1}$ is the **rate ratio** (or count ratio):
+$\exp{\beta_1}$ is the **rate ratio** (or count ratio):
 the multiplicative factor by which the expected count changes
 when $x_i$ increases from 0 to 1.
 
@@ -55,15 +55,15 @@ If $\beta_1 > 0$, the group with $x_i = 1$ has a higher expected count.
 
 For $x_i = 0$:
 $$
-\hat\mu_0 = e^{1.2} \approx 3.32
+\emu_0 = \exp{1.2} \approx 3.32
 $$
 
 For $x_i = 1$:
 $$
-\hat\mu_1 = e^{1.2 + 0.5} = e^{1.7} \approx 5.47
+\emu_1 = \exp{1.2 + 0.5} = \exp{1.7} \approx 5.47
 $$
 
-The estimated rate ratio is $e^{0.5} \approx 1.65$,
+The estimated rate ratio is $\exp{0.5} \approx 1.65$,
 meaning the group with $x_i = 1$ has about 65% more events on average.
 
 :::

_subfiles/count-regression/_exr-prac-glm-score.qmd

@@ -16,7 +16,7 @@ and
 $\deriv{\beta_1}\ell = 0$.
 
 **(c)** Interpret the score equations:
-what condition on the fitted values $\hat\mu_i$ do they imply?
+what condition on the fitted values $\emu_i$ do they imply?
 :::
 
 ::: {.solution}
@@ -33,7 +33,7 @@ $$
 $$
 
 where $\logf{\mu_i} = \beta_0 + \beta_1 x_i$,
-so $\mu_i = e^{\beta_0 + \beta_1 x_i}$.
+so $\mu_i = \exp{\beta_0 + \beta_1 x_i}$.
 
 **(b)**
 
@@ -60,19 +60,23 @@ $$
 
 **(c)**
 
-The first equation says $\sum_i y_i = \sum_i \hat\mu_i$:
+These equations hold at the maximum likelihood estimate $\evb$,
+where $\mu_i$ takes its fitted value
+$\emu_i \eqdef \exp{\eb_0 + \eb_1 x_i}$.
+
+The first equation says $\sum_i y_i = \sum_i \emu_i$:
 the total fitted count equals the total observed count.
 
-The second equation says $\sum_i x_i y_i = \sum_i x_i \hat\mu_i$:
+The second equation says $\sum_i x_i y_i = \sum_i x_i \emu_i$:
 the fitted counts are balanced against observed counts,
 weighted by $x_i$.
 
 More generally,
-these score equations say that the **residuals $(y_i - \hat\mu_i)$
+these score equations say that the **residuals $(y_i - \emu_i)$
 are [orthogonal](math-prereqs.qmd#def-orthogonal-vectors)
 to each predictor column**:
-for each predictor $j$, the residual vector $(\vy - \hat{\vec{\mu}})$ satisfies
-$\tp{\vx_{(j)}}(\vy - \hat{\vec{\mu}}) = 0$,
+for each predictor $j$, the residual vector $(\vy - \est{\vec{\mu}})$ satisfies
+$\tp{\vx_{(j)}}(\vy - \est{\vec{\mu}}) = 0$,
 where $\vx_{(j)} = (x_{1j}, \ldots, x_{nj})$ is the column of $j$-th
 predictor values across observations.
 This system of equations is the GLM analogue of the OLS normal equations.

_subfiles/count-regression/_sec-overdispersion.qmd

@@ -3,10 +3,11 @@
 ::: notes
 
 The Poisson distribution model **forces** the conditional variance
-to equal the conditional mean ($\Var{Y \mid \vX} = \Expp[Y \mid \vX]$).
+to equal the conditional mean
+($\Var{Y \mid \vX=\vx, T=t} = \Expp[Y \mid \vX=\vx, T=t]$).
 In practice, observational count data frequently exhibit variance
 substantially larger than the mean
-(or occasionally, smaller, termed underdispersion).
+(or occasionally smaller, termed underdispersion).
 
 :::
 
@@ -15,25 +16,47 @@ substantially larger than the mean
 :::: {#def-overdispersion}
 #### Overdispersion
 
-A random variable $Y$ is **overdispersed**
-relative to a parametric model $\P(Y=y \mid \vX)$ if
-its conditional empirical variance in a dataset exceeds
-the theoretical variance imposed by the fitted model $\est{\P}(Y=y \mid \vX)$.
+Write $m_{\P}(\vx, t)$ and $v_{\P}(\vx, t)$
+for the conditional mean and variance
+that a model $\P(Y=y \mid \vX=\vx, T=t)$ specifies;
+for a Poisson model, $v_{\P}(\vx, t) = m_{\P}(\vx, t) = \mu(\vx, t)$.
+Let $\P$ specify $Y$'s conditional mean correctly,
+so that $m_{\P}(\vx, t) = \Expp[Y \mid \vX=\vx, T=t]$
+for every $\vx$ and $t$.
+Then $Y$ is **overdispersed** relative to $\P$ if
+its conditional variance exceeds the one $\P$ specifies
+at some covariate pattern and exposure:
+
+$$\Var{Y \mid \vX=\vx, T=t} > v_{\P}(\vx, t)$$
+
+The same-mean requirement is what makes this a statement about dispersion:
+without it, any model that simply understates the mean
+would look overdispersed.
+In practice we detect overdispersion by comparing
+the conditional empirical variance in a dataset
+against $v_{\est{\P}}(\vx, t)$, the variance a *fitted* model predicts.
 
 ::::
 
 In Poisson regression, unmodeled heterogeneity, clustering,
-or omitted predictors cause overdispersion.
+and omitted predictors all inflate the conditional variance.
+Where the mean model remains correct,
+that inflation is overdispersion in the sense of @def-overdispersion.
+An omitted predictor may instead --- or additionally --- misspecify the mean,
+and to whatever extent it does,
+the model is failing the correct-mean requirement @def-overdispersion imposes,
+which is a different problem requiring a different remedy.
 When overdispersion is present but ignored,
-the point estimates $\evb$ remain unbiased,
+the point estimates $\evb$ remain consistent,
 but the standard errors produced by standard maximum likelihood estimation
 are severely underestimated.
 This underestimation leads to overly narrow confidence intervals
 and inflated false-positive (type I error) rates during hypothesis testing.
 
 When overdispersion is detected via residual diagnostics
-(such as a mean deviance or Pearson $\chi^2$ statistic
-significantly exceeding 1),
+(such as the deviance or Pearson $\chi^2$ statistic,
+divided by its residual degrees of freedom,
+substantially exceeding 1),
 practitioners can address it by incorporating missing predictors,
 using quasipoisson estimation, or fitting a negative binomial regression model.
 
@@ -44,8 +67,4 @@ c.f.
 @vittinghoff2e §8.1.5;
 and <https://en.wikipedia.org/wiki/Overdispersion>.
 
-When we encounter overdispersion,
-we can try to reduce the residual variance
-by adding more covariates.
-
 :::

_subfiles/count-regression/_sec-quasipoisson.qmd

@@ -0,0 +1,48 @@
+Another way to handle overdispersion ---
+rather than switching to the negative binomial distributional family ---
+is the quasipoisson approach.
+It assumes less *model structure* than the negative binomial does:
+rather than committing to a complete probability distribution
+and estimating by maximum likelihood,
+it assumes only the mean-variance relationship
+$\Var{Y \mid \vX=\vx, T=t} = \theta \mu(\vx, t)$,
+for a dispersion parameter $\theta$.
+It pairs that single assumption with a moment-based *inference method*,
+estimating $\theta$ from the Pearson residuals.
+
+While point estimates for regression coefficients $\evb$
+remain identical to standard Poisson regression,
+their estimated standard errors are scaled by $\sqrt{\eth}$.
+This approach provides valid standard errors and $p$-values
+when overdispersion is multiplicative.
+That validity comes from the assumed mean-variance relationship itself.
+A sandwich (robust) variance estimator also uses the residuals,
+but it does not require the model's variance function to be correct;
+it accumulates the squared residuals across observations
+instead of scaling the variance formula the model supplied.
+Its validity therefore does not depend on that variance function
+being the right one,
+only on the mean model being approximately correct
+[@vittinghoff2e, §4.7.3.6 for linear regression;
+@vittinghoff2e, §8.3.1 for the generalized linear model case].
+The quasipoisson scaling has no such guarantee:
+its single $\eth$ is estimated under the assumption of proportionality,
+so if the variance is not proportional to the mean,
+the scaled standard errors are simply wrong.
+
+Robust standard errors are not a free improvement, though.
+For linear regression, @vittinghoff2e [§4.7.3.6] reports simulations
+in which robust standard errors can be too small
+in samples as large as 250 observations,
+and recommends the more conservative HC3 estimator at those sizes.
+The specific remedy is a linear-model construction,
+but the caution about small samples is worth carrying over
+to the sample sizes many epidemiological studies actually have.
+
+The quasipoisson approach is simpler to implement
+than the negative binomial model,
+but provides less information than a full negative binomial likelihood:
+it does not specify a full parametric distribution
+for prediction intervals or model likelihood comparisons.
+
+See `?quasipoisson` in R for implementation details.

_subfiles/count-regression/_sec_pois-reg-preds.qmd

@@ -1,11 +1,15 @@
-
 $$
-\begin{aligned}
-\hat y
-&\eqdef \hat{\Expp}[Y|\vX= \vx,T=t]\\
-&=\hat\mu(\vx, t)\\
-&=\hat\lambda(\vx) \cdot t\\
-&=\exp{\hat\eta(\vx)} \cdot t\\
-&=\exp{\vx'\hat{\vec{\beta}}} \cdot t
-\end{aligned}
+\ba
+\ey
+&\eqdef \eExp{Y \mid \vX = \vx, T = t} \\
+&\quad \text{(definition of estimated conditional expectation)} \\
+&= \emu(\vx, t) \\
+&\quad \text{(estimated mean count function)} \\
+&= \el(\vx) \cdot t \\
+&\quad \text{(substituting the estimated rate relation)} \\
+&= \exp{\est{\eta}(\vx)} \cdot t \\
+&\quad \text{(substituting the inverse link function)} \\
+&= \exp{\eb_0 + \eb_1 x_1 + \dots + \eb_p x_p} \cdot t \\
+&\quad \text{(substituting estimated linear predictor } \est{\eta}(\vx)\text{)}
+\ea
 $$

_subfiles/count-regression/_sec_pois-reg_intro.qmd

@@ -5,26 +5,28 @@ This chapter presents models for
 [count data](data.qmd#sec-count-vars) outcomes.
 With covariates,
 the event rate $\lambda$
-becomes a function of the covariates
-$\vX = (X_1, \dots,X_n)$.
+becomes a function of the covariate vector
+$\vX = \tp{(X_1, \dots, X_p)} \in \reals^p$.
 Typically, count data models use
 a $\logf{}$ link function,
 and thus an $\exp{}$ inverse-link function.
-That is:
+Specifically, the model relates the expected outcome count
+to the event rate and linear predictor as:
 :::
 
 $$
-\begin{aligned}
-\Expp[Y | \vX = \vx, T = t] &= \mu(\vx,t)
-\\ \mu(\vx,t) &= \lambda(\vx)\cdot t
-\\ \lambda(\vx) &= \exp{\eta(\vx)}
-\\ \eta(\vx) &= \vx'\tilde \beta = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p
-\end{aligned}
+\ba
+\Expp[Y \mid \vX = \vx, T = t] &= \mu(\vx,t) \\
+\mu(\vx,t) &= \lambda(\vx)\cdot t \\
+\lambda(\vx) &= \exp{\eta(\vx)} \\
+\eta(\vx) &= \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p
+\ea
 $${#eq-mean-poisson}
 
 ::: notes
-$T=t$ is called the [exposure magnitude](probability.qmd#def-exposure)
-and has a special role in this model.
+The term $T=t$ represents the [exposure magnitude](probability.qmd#def-exposure)
+(such as person-years or observation time)
+and plays a structural role in scaling rates to expected counts.
 :::
 
 ---
@@ -49,26 +51,32 @@ in Binomial models.
 ---
 
 ::: notes
-We can also think of $t$ as a special part of the linear component:
+We can also express the exposure magnitude $t$
+directly as a component of the linear predictor:
 :::
 
 $$
-\begin{aligned}
-\logf{\Expp[Y | \vX = \vx,T=t]}
-&= \logf{\mu(\vx)}\\
-&=\logf{\lambda(\vx) \cdot t}\\
-&=\logf{\lambda(\vx)} + \log{t}\\
-&=\logf{\exp{\eta(\vx)}} + \log{t}\\
-&=\eta(\vx) + \log{t}\\
-&=\vx'\tilde\beta + \log{t}\\
-&=(\beta_0 +\beta_1 x_1+\dots + \beta_p x_p) + \log{t}\\
-\end{aligned}
+\ba
+\logf{\Expp[Y \mid \vX = \vx, T = t]}
+&= \logf{\mu(\vx, t)} \\
+&\quad \text{(definition of the conditional mean)} \\
+&= \logf{\lambda(\vx) \cdot t} \\
+&\quad \text{(substituting the rate relationship)} \\
+&= \logf{\lambda(\vx)} + \logf{t} \\
+&\quad \text{(logarithmic product rule)} \\
+&= \logf{\exp{\eta(\vx)}} + \logf{t} \\
+&\quad \text{(substituting the rate function)} \\
+&= \eta(\vx) + \logf{t} \\
+&\quad \text{(by inverse relationship of log and exp)} \\
+&= (\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p) + \logf{t} \\
+&\quad \text{(expanding the linear predictor } \eta(\vx)\text{)}
+\ea
 $$
 
 ::: notes
 In contrast with the other covariates (represented by $\vX$),
 $t$ enters this expression with a $\log{}$ transformation
-and without a corresponding $\beta$ coefficient;
+and without an estimated $\beta$ coefficient;
 in other words, $\logf{t}$ is an [offset term](probability.qmd#def-offset).
 :::
 
@@ -82,14 +90,14 @@ What are the units of $\mu$ in @eq-mean-poisson?
 
 :::{#sol-eq-mean-poisson}
 
-$\mu$ is the mean of $Y$,
-and $Y$ is a count,
-so $\mu$ is also a count;
+$\mu$ is the expected value of $Y$.
+Because $Y$ represents a count of events,
+$\mu$ is expressed in units of event counts;
 for example:
 
 - 3.1 cyclones,
-- 10.23 ER visits
-- 15.01 infections
+- 10.23 ER visits,
+- 15.01 infections.
 
 :::
 
@@ -103,11 +111,11 @@ What are the units of $\lambda$ in @eq-mean-poisson?
 
 :::{#sol-eq-rate-poisson}
 $\lambda = \mu/t$,
-so $\lambda$ is a rate of counts per unit of $t$.
+so $\lambda$ is an event rate per unit of exposure $t$.
 For example:
 
-- 3.1 cyclones *per year*
-- 2.023 ER visits per 10 person-years
-- 15.01 infections per 1000 person-years at risk
+- 3.1 cyclones *per year*,
+- 2.023 ER visits per 10 person-years,
+- 15.01 infections per 1000 person-years at risk.
 
 :::