diff --git a/_subfiles/count-regression/_exr-prac-glm-interp.qmd b/_subfiles/count-regression/_exr-prac-glm-interp.qmd index 02554c9d2..dea7898be 100644 --- a/_subfiles/count-regression/_exr-prac-glm-interp.qmd +++ b/_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. ::: diff --git a/_subfiles/count-regression/_exr-prac-glm-score.qmd b/_subfiles/count-regression/_exr-prac-glm-score.qmd index 6356090a5..da12ff1a6 100644 --- a/_subfiles/count-regression/_exr-prac-glm-score.qmd +++ b/_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. diff --git a/_subfiles/count-regression/_sec-overdispersion.qmd b/_subfiles/count-regression/_sec-overdispersion.qmd index b460f3286..289e1c7f6 100644 --- a/_subfiles/count-regression/_sec-overdispersion.qmd +++ b/_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 . -When we encounter overdispersion, -we can try to reduce the residual variance -by adding more covariates. - ::: diff --git a/_subfiles/count-regression/_sec-quasipoisson.qmd b/_subfiles/count-regression/_sec-quasipoisson.qmd new file mode 100644 index 000000000..0faf45edd --- /dev/null +++ b/_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. diff --git a/_subfiles/count-regression/_sec_pois-reg-preds.qmd b/_subfiles/count-regression/_sec_pois-reg-preds.qmd index 6fcd4f9f8..9904c90c6 100644 --- a/_subfiles/count-regression/_sec_pois-reg-preds.qmd +++ b/_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 $$ diff --git a/_subfiles/count-regression/_sec_pois-reg_intro.qmd b/_subfiles/count-regression/_sec_pois-reg_intro.qmd index af9303f79..f6c890a4b 100644 --- a/_subfiles/count-regression/_sec_pois-reg_intro.qmd +++ b/_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. ::: diff --git a/_subfiles/count-regression/_sec_poisson_RRs.qmd b/_subfiles/count-regression/_sec_poisson_RRs.qmd index 5b09d92bc..ed1c2df00 100644 --- a/_subfiles/count-regression/_sec_poisson_RRs.qmd +++ b/_subfiles/count-regression/_sec_poisson_RRs.qmd @@ -2,7 +2,7 @@ Applying the [general procedure for interpreting a regression coefficient](Linea ::: notes Differences on the log-rate scale become ratios on the rate scale, -because +because exponentiating a difference yields a quotient: ::: $$\exp{a-b} = \frac{\exp{a}}{\exp{b}}$$ @@ -10,31 +10,42 @@ $$\exp{a-b} = \frac{\exp{a}}{\exp{b}}$$ (recall from [Algebra 2](math-prereqs.qmd#cor-exp-sum)) Therefore, according to this model, -**differences of $\delta$ in covariate $x_j$ correspond to rate ratios of $\exp{\b_j \cdot \delta}$**. +**a difference of $\delta \eqdef a - b$ +between two values $a$ and $b$ of covariate $x_j$ +corresponds to a rate ratio of $\exp{\b_j \cdot \delta}$**. -Specifically, letting $\vX_{-j}$ denote vector $\vX$ with element $j$ removed: +Specifically, let $\vX_{-j}$ denote the vector of all covariates except $X_j$, +and abbreviate the two expected counts being compared as $$ \ba -&\phantom{={}} \logf{\E{Y |\red{X_j = a}, \vX_{-j}=\vx_{-j},T=t}} \\ -&\phantom{={}} -\logf{\E{Y |\red{X_j = b}, \vX_{-j}=\vx_{-j},T=t}} \\ -&= \logf{t} + \b_0 + \b_1 x_1 + \lds + \red{\b_j (a)} + \lds + \b_p x_p \\ -&\phantom{={}} -\logf{t} - \b_0 - \b_1 x_1 - \lds - \red{\b_j (b)} - \lds - \b_p x_p \\ -&= \red{\b_j(a-b)} +\mu_a &\eqdef \E{Y \mid \red{X_j = a}, \vX_{-j} = \vx_{-j}, T = t} \\ +\mu_b &\eqdef \E{Y \mid \red{X_j = b}, \vX_{-j} = \vx_{-j}, T = t} \ea $$ -And accordingly, +$$ +\ba +&\phantom{={}} \logf{\mu_a} - \logf{\mu_b} \\ +&= \paren{\logf{t} + \b_0 + \b_1 x_1 + \lds + \red{\b_j a} + \lds + \b_p x_p} \\ +&\phantom{={}} - \paren{\logf{t} + \b_0 + \b_1 x_1 + \lds + \red{\b_j b} + \lds + \b_p x_p} \\ +&\quad \text{(substituting the linear predictor)} \\ +&= \red{\b_j a} - \red{\b_j b} \\ +&\quad \text{(canceling terms shared by both patterns)} \\ +&= \red{\b_j(a-b)} \\ +&\quad \text{(factoring out coefficient } \b_j\text{)} +\ea +$$ + +The rate ratio between the two covariate patterns +is therefore the exponential of that difference in log expectations: $$ \ba -\frac -{\E{Y |\red{X_j = a}, \vX_{-j} = \vx_{-j}, T = t} -} -{ -\E{Y |\red{X_j = b}, \vX_{-j}=\vx_{-j},T=t} -} -= -\exp{\red{\b_j(a-b)}} +&\phantom{={}} \frac{\mu_a}{\mu_b} \\ +&= \exp{\logf{\mu_a} - \logf{\mu_b}} \\ +&\quad \text{(by identity } \frac{u}{v} = \exp{\logf{u} - \logf{v}}\text{)} \\ +&= \exp{\red{\b_j(a-b)}} \\ +&\quad \text{(substituting the difference in log expectations)} \ea $$ diff --git a/_subfiles/count-regression/_sec_poisson_dx.qmd b/_subfiles/count-regression/_sec_poisson_dx.qmd index 21ff25aac..f4377cfd8 100644 --- a/_subfiles/count-regression/_sec_poisson_dx.qmd +++ b/_subfiles/count-regression/_sec_poisson_dx.qmd @@ -2,27 +2,40 @@ #### Observation residuals -$$e \eqdef y - \hat y$$ +$$e_i \eqdef y_i - \ey_i$$ #### Pearson residuals -$$r \eqdef \frac{e}{\hse{e}} \approx \frac{e}{\sqrt{\hat y}}$$ +$$r_i \eqdef \frac{e_i}{\hse{e_i}} \approx \frac{e_i}{\sqrt{\ey_i}}$$ #### Standardized Pearson residuals -$$r_p \eqdef \frac{r}{\sqrt{1-h}}$$ -where $h$ is the "leverage" (which we will continue to leave undefined). +$$r_{p,i} \eqdef \frac{r_i}{\sqrt{1-h_i}}$$ +where $h_i$ is the +[leverage](logistic-regression.qmd#def-leverage-glm) of observation $i$: +the $i$-th diagonal element of the weighted hat matrix. +That definition is stated for a logistic model, +so it uses binomial weights; +for Poisson regression the corresponding weight is $\emu_i$. +That definition also indexes covariate patterns, +which here we take to be the individual observations. #### Deviance residuals $$ -d_k \eqdef \signt(y - \hat y) - \left\{\sqrt{2[\ell_{\text{full}}(y) - \ell(\hat\beta; y)]}\right\} +d_i \eqdef \signt(y_i - \ey_i) +\sqrt{2\sb{\llik_{\text{full}}(y_i) - \llik(\evb; y_i)}} $$ :::{.callout-note} -$$\signt(x) \eqdef \frac{x}{|x|}$$ +$$ +\signt(x) \eqdef +\begin{cases} +\frac{x}{|x|} & x \neq 0\\ +0 & x = 0 +\end{cases} +$$ In other words: * $\signt(x) = -1$ if $x < 0$ @@ -32,7 +45,7 @@ In other words: ::::{.content-hidden} ```{r} -plot(sign,xlim = c(-1,1), xlab = "x", ylab = "sign(x)") +plot(sign, xlim = c(-1, 1), xlab = "x", ylab = "sign(x)") ``` :::: diff --git a/_subfiles/count-regression/_sec_poisson_inference.qmd b/_subfiles/count-regression/_sec_poisson_inference.qmd index 5ddde4031..7381bbd86 100644 --- a/_subfiles/count-regression/_sec_poisson_inference.qmd +++ b/_subfiles/count-regression/_sec_poisson_inference.qmd @@ -3,30 +3,30 @@ A Wald 95% confidence interval for a single coefficient $\beta_j$ is: $$ -\beta_j \in \sb{\hat\beta_j \pm \ciradf{\hat\beta_j}} +\beta_j \in \sb{\eb_j \pm \ciradf{\eb_j}} $$ where $z_{1-\alpha/2} \approx 1.96$ for $\alpha = 0.05$. Because the log-rate scale is related to the rate scale by exponentiation, -we obtain a confidence interval for the rate ratio $e^{\beta_j}$ +we obtain a confidence interval for the rate ratio $\exp{\beta_j}$ by exponentiating both endpoints: $$ -e^{\beta_j} \in +\exp{\beta_j} \in \sb{ - \exp{\hat\beta_j - \ciradf{\hat\beta_j}},\; - \exp{\hat\beta_j + \ciradf{\hat\beta_j}} + \exp{\eb_j - \ciradf{\eb_j}},\; + \exp{\eb_j + \ciradf{\eb_j}} } $$ ### Hypothesis tests for regression coefficients -To test $H_0: \beta_j = \beta_0$ against a one- or two-sided alternative, +To test $H_0: \beta_j = \beta_{j,0}$ against a one- or two-sided alternative, compute the Wald $z$-statistic: $$ -z = \frac{\hat \beta_j - \beta_0}{\hse{\hat\beta_j}} +z = \frac{\eb_j - \beta_{j,0}}{\hse{\eb_j}} $$ and compare $z$ (one-sided) or $|z|$ (two-sided) to the tails of the @@ -40,13 +40,11 @@ To compare a smaller model $M_0$ (with $p_0$ parameters) to a larger model $M_1$ (with $p_1 > p_0$ parameters), use the likelihood ratio test statistic: $$ -G^2 = 2\bigl[\hat\ell_1 - \hat\ell_0\bigr] +G^2 = 2\sb{\est{\llik}_1 - \est{\llik}_0} $$ -where $\hat\ell_1$ and $\hat\ell_0$ are the maximized log-likelihoods +where $\est{\llik}_1$ and $\est{\llik}_0$ are the maximized log-likelihoods of $M_1$ and $M_0$ respectively. -(Here the subscripts index the two *models*; they are unrelated to the -scalar null value $\beta_0$ used in the Wald test.) Under $H_0$ that the additional $p_1 - p_0$ parameters are all zero, $G^2 \dsim \chi^2_{p_1 - p_0}$. diff --git a/_subfiles/count-regression/_sec_zero-inflation-moments.qmd b/_subfiles/count-regression/_sec_zero-inflation-moments.qmd new file mode 100644 index 000000000..fe232d354 --- /dev/null +++ b/_subfiles/count-regression/_sec_zero-inflation-moments.qmd @@ -0,0 +1,105 @@ +::: {#exr-zinf-moments} + +Derive the expected value and variance of $Y$, +conditional on $\vX=\vx$ and $T=t$, +as functions of $\pi$ and $\mu_0$. +::: + +::: {.solution} + +**Expected value.** +By the Law of Total Expectation +(conditioning on $Z$, within the subpopulation $\{\vX=\vx, T=t\}$): + +$$ +\ba +\Expp[Y \mid \vX=\vx, T=t] +&= \pi \, \Expp[Y \mid Z=1, \vX=\vx, T=t] \\ +&\phantom{={}} + (1-\pi) \, \Expp[Y \mid Z=0, \vX=\vx, T=t] \\ +&\quad \text{(by Law of Total Expectation; } Z \ind T \mid \vX\text{)} \\ +&= 0 \cdot \pi + \mu_0 (1-\pi) \\ +&\quad \text{(substituting the conditional means)} \\ +&= (1-\pi) \mu_0 \\ +&\quad \text{(simplifying arithmetic)} +\ea +$$ + +The substitution $\Expp[Y \mid Z=0, \vX=\vx, T=t] = \mu_0$ follows immediately +from the definition of $\mu_0$. + +**Variance.** +Within this derivation, +write $\Expp[\,\cdot \mid Z]$ and $\Var{\cdot \mid Z}$ +for the moments conditional on $Z$ **and** on $\vX=\vx, T=t$; +the outer operators keep their conditioning explicit. +By the Law of Total Variance: + +$$ +\ba +\Var{Y \mid \vX=\vx, T=t} +&= \Expp\sb{\Var{Y \mid Z} \mid \vX=\vx, T=t} \\ +&\phantom{={}} + \Var{\Expp[Y \mid Z] \mid \vX=\vx, T=t} \\ +&\quad \text{(by Law of Total Variance)} +\ea +$$ + +For the expected conditional variance term, note that + +$$\Var{Y \mid Z=1} = 0 \quad\text{and}\quad \Var{Y \mid Z=0} = \mu_0$$ + +since the $Z=0$ arm is Poisson, so: + +$$ +\ba +\Expp\sb{\Var{Y \mid Z} \mid \vX=\vx, T=t} +&= \pi \, \Var{Y \mid Z=1} + (1-\pi) \, \Var{Y \mid Z=0} \\ +&\quad \text{(expectation over } Z\text{; } Z \ind T \mid \vX\text{)} \\ +&= 0 \cdot \pi + \mu_0 (1-\pi) \\ +&\quad \text{(substituting the conditional variances)} \\ +&= (1-\pi)\mu_0 \\ +&\quad \text{(simplifying arithmetic)} +\ea +$$ + +For the variance of conditional expectation term, +$\Expp[Y \mid Z]$ takes value 0 (with probability $\pi$) +or $\mu_0$ (with probability $1-\pi$), so: + +$$ +\ba +\Var{\Expp[Y \mid Z] \mid \vX=\vx, T=t} +&= \pi \paren{0 - (1-\pi)\mu_0}^2 + (1-\pi) \paren{\mu_0 - (1-\pi)\mu_0}^2 \\ +&\quad \text{(variance over } Z\text{; } Z \ind T \mid \vX\text{)} \\ +&= \pi(1-\pi)^2 \mu_0^2 + (1-\pi)\pi^2 \mu_0^2 \\ +&\quad \text{(expanding squared terms)} \\ +&= \pi(1-\pi)\mu_0^2 \sb{(1-\pi) + \pi} \\ +&\quad \text{(factoring)} \\ +&= \pi(1-\pi)\mu_0^2 \\ +&\quad \text{(since } (1-\pi) + \pi = 1\text{)} +\ea +$$ + +Combining both terms gives: + +$$ +\ba +\Var{Y \mid \vX=\vx, T=t} +&= (1-\pi)\mu_0 + \pi(1-\pi)\mu_0^2 \\ +&\quad \text{(summing expected variance and variance of expectation)} \\ +&= (1-\pi)\mu_0 \paren{1 + \pi\mu_0} \\ +&\quad \text{(factoring out } (1-\pi)\mu_0\text{)} +\ea +$$ + +The logistic model puts $\pi$ strictly between 0 and 1, +since $\expit$ never attains its limits, +and $\mu_0 = t \exp{\eta(\vx)} > 0$ whenever $t > 0$. +Then $1 + \pi\mu_0 > 1$, so + +$$(1-\pi)\mu_0 (1+\pi\mu_0) > (1-\pi)\mu_0 = \Expp[Y \mid \vX=\vx, T=t]$$ + +and a zero-inflated count model is overdispersed +relative to a Poisson model with the same mean, +at every covariate pattern with positive exposure. + +::: diff --git a/_subfiles/count-regression/_sec_zero-inflation.qmd b/_subfiles/count-regression/_sec_zero-inflation.qmd index b8347fbbd..aae1bc5c4 100644 --- a/_subfiles/count-regression/_sec_zero-inflation.qmd +++ b/_subfiles/count-regression/_sec_zero-inflation.qmd @@ -5,157 +5,110 @@ which we model using logistic regression: $$ \ba -\P(Z=1 \mid \vX=\vx) -&\eqdef \pi(\vx)\\ -&= \expit(\gamma_0 + \gamma_1 x_1 + \dots) +\pi(\vx) +&\eqdef \P(Z=1 \mid \vX=\vx)\\ +&= \expit(\gamma_0 + \gamma_1 x_1 + \dots + \gamma_p x_p) \ea $$ +The model makes $Z$ depend on the covariates alone, +not on the exposure magnitude $T$, +so $\P(Z=1 \mid \vX=\vx, T=t) = \pi(\vx)$ for every $t$ +--- that is, $Z \ind T \mid \vX$. + According to this model, if $Z=1$, then $Y$ will always be zero, regardless of $\vX$ and $T$: $$\P(Y=0 \mid Z=1, \vX=\vx, T=t) = 1$$ -Otherwise (if $Z=0$), $Y$ will have a Poisson distribution, +Otherwise (if $Z=0$), $Y$ follows a Poisson distribution, conditional on $\vX$ and $T$, as in a standard Poisson regression model. +Throughout this section, abbreviate +$\pi \eqdef \pi(\vx)$ +and $\mu_0 \eqdef \Expp[Y \mid Z=0, \vX=\vx, T=t]$. +Since the $Z=0$ arm is an ordinary Poisson regression model, +$\mu_0 = t \exp{\eta(\vx)}$ by @eq-mean-poisson, +which is strictly positive whenever $t > 0$. + Even though we never observe $Z$, we can estimate the parameters $\gamma_0, \dots, \gamma_p$ via maximum likelihood: $$ -\begin{aligned} -\P(Y=y \mid \vX=\vx, T=t) &= \P(Y=y, Z=1 \mid \vX=\vx, T=t) \\ -&\quad + \P(Y=y, Z=0 \mid \vX=\vx, T=t) -\end{aligned} +\ba +\P(Y=y \mid \vX=\vx, T=t) +&= \P(Y=y, Z=1 \mid \vX=\vx, T=t) \\ +&\phantom{={}} + \P(Y=y, Z=0 \mid \vX=\vx, T=t) \\ +&\quad \text{(by Law of Total Probability)} +\ea $$ -(by the Law of Total Probability) where $$ -\begin{aligned} +\ba \P(Y=y, Z=z \mid \vX=\vx, T=t) -&= \P(Y=y \mid Z=z, \vX=\vx, T=t)\,\P(Z=z \mid \vX=\vx) -\end{aligned} +&= \P(Y=y \mid Z=z, \vX=\vx, T=t) \P(Z=z \mid \vX=\vx) \\ +&\quad \text{(with } Z \ind T \mid \vX\text{)} +\ea $$ ---- +{{< slidebreak >}} ::: {#exr-zinf-pmf} Expand $\P(Y=0 \mid \vX=\vx, T=t)$, $\P(Y=1 \mid \vX=\vx, T=t)$, and $\P(Y=y \mid \vX=\vx, T=t)$ -into expressions involving $\P(Z=1 \mid \vX=\vx)$ -and $\P(Y=y \mid Z=0, \vX=\vx, T=t)$. +into expressions involving $\pi$ and $\mu_0$. ::: ::: {.solution} -Let $\pi \eqdef \P(Z=1 \mid \vX=\vx)$ -and $\mu_0 \eqdef \Expp[Y \mid Z=0, \vX=\vx, T=t]$. - **$\P(Y=0)$:** $Y=0$ occurs either because $Z=1$ (always zero) or because $Z=0$ and the Poisson draw equals 0: $$ \ba \P(Y=0 \mid \vX=\vx, T=t) -&= \P(Z=1 \mid \vX=\vx) + \P(Z=0 \mid \vX=\vx)\,\P(Y=0 \mid Z=0, \vX=\vx, T=t)\\ -&= \pi + (1-\pi)\,e^{-\mu_0} +&= \pi \, \P(Y=0 \mid Z=1, \vX=\vx, T=t) \\ +&\phantom{={}} + (1-\pi) \, \P(Y=0 \mid Z=0, \vX=\vx, T=t) \\ +&\quad \text{(by Law of Total Probability; } Z \ind T \mid \vX\text{)} \\ +&= \pi \cdot 1 + (1-\pi) \exp{-\mu_0} \\ +&\quad \text{(substituting the two conditional PMFs)} \\ +&= \pi + (1-\pi) \exp{-\mu_0} \\ +&\quad \text{(simplifying arithmetic)} \ea $$ **$\P(Y=1)$:** $Z=1$ can never produce $Y=1$, so: -$$ -\P(Y=1 \mid \vX=\vx, T=t) -= (1-\pi)\,\P(Y=1 \mid Z=0, \vX=\vx, T=t) -= (1-\pi)\,\mu_0 e^{-\mu_0} -$$ - -**$\P(Y=y)$ for $y \geq 1$:** Identical reasoning gives - -$$ -\P(Y=y \mid \vX=\vx, T=t) -= (1-\pi)\,\frac{\mu_0^y e^{-\mu_0}}{y!} -$$ - -::: - -{{< slidebreak >}} - -::: {#exr-zinf-moments} - -Derive the expected value and variance of $Y$, conditional on $\vX$ and $T$, -as functions of $\pi \eqdef \P(Z=1 \mid \vX=\vx)$ -and $\mu_0 \eqdef \Expp[Y \mid Z=0, \vX=\vx, T=t]$. -::: - -::: {.solution} - -Let $\pi \eqdef \P(Z=1 \mid \vX=\vx)$ -and $\mu_0 \eqdef \Expp[Y \mid Z=0, \vX=\vx, T=t]$. - -**Expected value.** By the Law of Total Expectation -(conditioning on $Z$, within the subpopulation $\{\vX=\vx, T=t\}$): - $$ \ba -\Expp[Y \mid \vX=\vx, T=t] -&= \Expp[Y \mid Z=1, \vX=\vx, T=t]\,\pi \\ -&\quad + \Expp[Y \mid Z=0, \vX=\vx, T=t]\,(1-\pi)\\ -&= 0 \cdot \pi + \mu_0(1-\pi)\\ -&= (1-\pi)\,\mu_0 +\P(Y=1 \mid \vX=\vx, T=t) +&= \pi \, \P(Y=1 \mid Z=1, \vX=\vx, T=t) \\ +&\phantom{={}} + (1-\pi) \, \P(Y=1 \mid Z=0, \vX=\vx, T=t) \\ +&\quad \text{(by Law of Total Probability; } Z \ind T \mid \vX\text{)} \\ +&= \pi \cdot 0 + (1-\pi) \mu_0 \exp{-\mu_0} \\ +&\quad \text{(since } \P(Y=1 \mid Z=1, \vX=\vx, T=t)=0\text{)} \\ +&= (1-\pi) \mu_0 \exp{-\mu_0} \\ +&\quad \text{(simplifying arithmetic)} \ea $$ -The substitution $\Expp[Y \mid Z=0, \vX=\vx, T=t] = \mu_0$ follows immediately -from the definition of $\mu_0$. - -**Variance.** By the Law of Total Variance. -To reduce clutter, we suppress the $(\vX=\vx, T=t)$ conditioning in the -intermediate steps: every expectation and variance is taken within -the subpopulation $\{\vX=\vx, T=t\}$, and we restore the explicit conditioning -in the final line. - -$$ -\Var{Y} = \Expp[\Var{Y \mid Z}] + \Var{\Expp[Y \mid Z]} -$$ - -For the first term, -since $\Var{Y \mid Z=1}=0$ and $\Var{Y \mid Z=0}=\mu_0$ (Poisson): - -$$ -\Expp[\Var{Y \mid Z}] = 0 \cdot \pi + \mu_0(1-\pi) = (1-\pi)\mu_0 -$$ - -For the second term, $\Expp[Y \mid Z]$ takes the value 0 (with prob $\pi$) -or $\mu_0$ (with prob $1-\pi$), so: +**$\P(Y=y)$ for $y \geq 1$:** Identical reasoning gives: $$ \ba -\Var{\Expp[Y \mid Z]} -&= \pi(0 - (1-\pi)\mu_0)^2 + (1-\pi)(\mu_0 - (1-\pi)\mu_0)^2\\ -&= \pi(1-\pi)^2\mu_0^2 + (1-\pi)\pi^2\mu_0^2\\ -&= \pi(1-\pi)\mu_0^2[\,(1-\pi)+\pi\,]\\ -&= \pi(1-\pi)\mu_0^2 +\P(Y=y \mid \vX=\vx, T=t) +&= \pi \, \P(Y=y \mid Z=1, \vX=\vx, T=t) \\ +&\phantom{={}} + (1-\pi) \, \P(Y=y \mid Z=0, \vX=\vx, T=t) \\ +&\quad \text{(by Law of Total Probability; } Z \ind T \mid \vX\text{)} \\ +&= \pi \cdot 0 + (1-\pi) \frac{\mu_0^y \exp{-\mu_0}}{y!} \\ +&\quad \text{(since } \P(Y=y \mid Z=1, \vX=\vx, T=t)=0 \text{ for } y \geq 1\text{)} \\ +&= (1-\pi) \frac{\mu_0^y \exp{-\mu_0}}{y!} \\ +&\quad \text{(simplifying arithmetic)} \ea $$ -Combining: - -$$ -\Var{Y \mid \vX=\vx, T=t} = (1-\pi)\mu_0 + \pi(1-\pi)\mu_0^2 -= (1-\pi)\mu_0\bigl(1 + \pi\mu_0\bigr) -$$ - -Since -$$ -(1-\pi)\mu_0(1+\pi\mu_0) \geq (1-\pi)\mu_0 = \Expp[Y \mid \vX=\vx, T=t] -$$ -for any $\pi > 0$, -zero-inflated count models always exhibit overdispersion -relative to a Poisson model with the same mean. - ::: diff --git a/chapters/count-regression.qmd b/chapters/count-regression.qmd index d2a8fff0a..86ca82dee 100644 --- a/chapters/count-regression.qmd +++ b/chapters/count-regression.qmd @@ -62,6 +62,10 @@ This content is adapted from: ### Models for zero-inflated counts {{< include _subfiles/count-regression/_sec_zero-inflation.qmd >}} +{{< slidebreak >}} + +{{< include _subfiles/count-regression/_sec_zero-inflation-moments.qmd >}} + # Over-dispersion {{< include _subfiles/count-regression/_sec-overdispersion.qmd >}} @@ -79,22 +83,25 @@ Most notably, the [negative binomial model](probability.qmd#sec-nb-dist). ::: -When diagnostics reveal overdispersion, -the standard Poisson assumption $\Var{Y \mid \vX} = \Expp[Y \mid \vX]$ -is violated, -leading to artificially narrow standard errors +Diagnostics that reveal overdispersion are telling us +that the data violate the standard Poisson assumption + +$$\Var{Y \mid \vX=\vx, T=t} = \Expp[Y \mid \vX=\vx, T=t]$$ + +Ignoring that violation leads to artificially narrow standard errors and inflated type I error rates. The [negative binomial distribution](probability.qmd#sec-nb-dist) serves as a natural generalization of the Poisson distribution for count outcomes. It introduces an overdispersion parameter $\rho$ that allows the conditional variance to exceed the mean -($\Var{Y \mid \vX} = \mu + \mu^2 / \rho$). +($\Var{Y \mid \vX=\vx, T=t} = \mu(\vx, t) + \mu(\vx, t)^2 / \rho$). We still model $\mu(\vx, t) = t \exp{\eta(\vx)}$ as before, preserving the rate-ratio interpretation for regression coefficients. Furthermore, negative binomial models can be combined with zero-inflation -to account for both structural zeros -and variance expansion in count data. +to account for both structural zeros and variance expansion in count data, +with the negative binomial serving as the conditional distribution +for the count component. --- @@ -106,24 +113,7 @@ and variance expansion in count data. ## Quasipoisson regression -Another flexible approach to address overdispersion -— without changing the underlying Poisson point estimates — -is quasi-likelihood estimation via the quasipoisson model family. -Rather than specifying a complete likelihood function -and fitting by maximum likelihood, -quasipoisson models specify only the mean-variance relationship -$\Var{Y \mid \vX} = \theta \mu(\vx, t)$, -where $\theta$ is a dispersion parameter estimated from 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 method-of-moments approach provides robust standard error estimates -and valid $p$-values when overdispersion is multiplicative, -though it does not specify a full parametric distribution -for prediction intervals or model likelihood comparisons. - -See `?quasipoisson` in R for implementation details. +{{< include _subfiles/count-regression/_sec-quasipoisson.qmd >}} # More on count regression