From 061702775a218739ec3a4cc5f5347a44844a2802 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 25 Mar 2025 11:47:53 -0500 Subject: [PATCH 1/6] add documentation page for each model --- docs/make.jl | 1 + docs/src/models.md | 115 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 docs/src/models.md diff --git a/docs/make.jl b/docs/make.jl index 6afdaa8b..7bcce664 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,6 +6,7 @@ makedocs( modules = [GLM], pages = [ "Home" => "index.md", + "Supported Models" => "models.md", "examples.md", "api.md", ], diff --git a/docs/src/models.md b/docs/src/models.md new file mode 100644 index 00000000..0ae81847 --- /dev/null +++ b/docs/src/models.md @@ -0,0 +1,115 @@ +# Supported Models + +Broadly, there are two types of models provided by this package: linear models and generalized linear models (GLM). +In fact, (regular) linear models are just a special case of GLM. + +## Linear Model + +Suppose we have a response variable $y_i$ and explanatory variables $x_{1i}, x_{2i}, ...$ +In its simplest form, the linear model uses the following: +```math +(Y_i\ \vert \ \mathbf{X}_i = \mathbf{x}_i)\ \sim\ \mathcal{N}(\mathbf{x}_i^\top \mathbb{\beta}, \sigma^2) +``` + + +### Weighted Linear Model + +This package also supports weighted .... + +## Generalized Linear Model (GLM) + +For GLMs, we have the same setup as before however the response $y_i$ is (usually) constrained in some way. +For instance, if we have count data $y_i$'s are non-negative integers. +GLMs model the parameters through a transformation of the mean that _links_ it to a linear predictor. +The function that operationalizes this is _link_ function (usually denoted by $g$). +```math +g(\mathbb{E}(Y_i | \mathbf{x}_i)) = \mathbf{x}_i^\top \mathbb{\beta} = \eta_i . +``` +The link function must be invertible, and the inverse link function is also sometimes referred to as the mean function as it translates the linear predictor into the mean, that is: +```math + \mathbb{E}(Y_i | \mathbf{x}_i) = g^{-1}(\eta_i) . +``` + +Thus, the link function allows us to go from the mean to the associated linear predictor and the mean function allows us to go from the linear predictor to the mean. + + +## Binary Response + +A binary response $y_i$ can be modeled as: +```math +y_i \sim \operatorname{Bernoulli}(p_i) , +``` +where the $p_i$ is the mean, i.e. $\mathbb{E}(Y_i | \mathbf{x}_i) = \mathbb{P}(Y_i = 1 | \mathbf{x}_i) = p_i$ is modeled with a linear predictor using a link function as: +```math +g(p_i) = \mathbf{x}_i^\top \mathbb{\beta}. +``` + +To fit a GLM for binary responses, use: +```julia +mod = glm(@formula(...), data, Bernoulli()) +``` + +Note that the above function call did not specify a link function. +That is because, the _canonical_ link function is used as the default. +Briefly, the canonical link function simplifies the model. + +### Logistic Regression: The canonical link function + +The _canonical_ link function for this model is the _logit_ function[^logit]: +$$ +g(p) = \log \left(\frac{p}{1-p}\right) . +$$ + +Using this, we get the relationship: +```math +\log \left(\frac{p_i}{1-p_i}\right) = \mathbf{x}_i^\top \mathbb{\beta}, +``` + +and the corresponding mean function is the _logistic_ function[^logistic] $\sigma(\eta) = (1 + e^{-\eta})^{-1}$, that is: +```math + p_i = \sigma(\mathbf{x}_i^\top \boldsymbol{\beta}) = \frac{1}{1 + e^{-\mathbf{x}_i^\top \boldsymbol{\beta}}} +``` + +The _logit_ link function is the default used for binary response is implemented as the `LogitLink` object, ans + + +### Probit Regression: Bernoulli GLM with the Probit link function + +Probit regression uses the cumulative distribution function of the standard normal distribution as the _mean_ function. +Let $\Phi(\eta)$ denote the standard normal CDF, then the link function is $\Phi^{-1}(\cdot)$, that is: +```math +\Phi^{-1}(p_i) = \mathbf{x}_i^\top \boldsymbol{\beta}, +``` +and the mean function is simply: +```math +p_i = \Phi(\mathbf{x}_i^\top \boldsymbol{\beta}) +``` + +The probit link is implemented as the `ProbitLink` type and can be performed using: +```julia +mod = glm(@formula(...), data, Bernoulli(), ProbitLink()) +``` + +### Poisson Regression for count data + +When we have count data, i.e. the observed $y_i$'s are non-negative integers we can use the following model: +$$ + ( Y_i\ \vert\ \mathbf{X}_i = mathbf{x}_i) \sim operatorname{Poisson}(\lambda_i). +$$ + +Since $\mathb{E}(Y_i\ \vert\ mathbf{x}_i) = \lambda_i$ this can be formulated as a GLM with _log_ as the link function: +```math +\log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} . +``` +The corresponding mean function is: +```math +\lambda_i = e^{\mathbf{x}_i^\top \boldsymbol{\beta} } . +``` +This can be implemented as: +```julia +mod = glm(@formula(...), data, Poisson(), LogLink()) +``` + + +[^logit]: https://en.wikipedia.org/wiki/Logit +[^logistic]: https://en.wikipedia.org/wiki/Logistic_function \ No newline at end of file From 5f6a229a798d33c2291c1c33a9a1b381ca4ecce4 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 25 Mar 2025 11:55:39 -0500 Subject: [PATCH 2/6] added three-components part --- docs/src/models.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docs/src/models.md b/docs/src/models.md index 0ae81847..01e088c4 100644 --- a/docs/src/models.md +++ b/docs/src/models.md @@ -32,6 +32,11 @@ The link function must be invertible, and the inverse link function is also some Thus, the link function allows us to go from the mean to the associated linear predictor and the mean function allows us to go from the linear predictor to the mean. +In short, a GLM consists of three key components[^GLMwiki]: + +1. The distribution +2. The linear predictor $(\eta = \mathbf{x}^\top \boldsymbol{\beta})$ +3. The link function $g(\mu) = \eta$ ## Binary Response @@ -112,4 +117,5 @@ mod = glm(@formula(...), data, Poisson(), LogLink()) [^logit]: https://en.wikipedia.org/wiki/Logit -[^logistic]: https://en.wikipedia.org/wiki/Logistic_function \ No newline at end of file +[^logistic]: https://en.wikipedia.org/wiki/Logistic_function +[^GLMwiki]: https://en.wikipedia.org/wiki/Generalized_linear_model#Model_components From 999c40921b2538a8ad093ac8104546048503b740 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 25 Mar 2025 12:09:26 -0500 Subject: [PATCH 3/6] fleshed out GLM intro --- docs/src/models.md | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/docs/src/models.md b/docs/src/models.md index 01e088c4..5c0b317e 100644 --- a/docs/src/models.md +++ b/docs/src/models.md @@ -22,6 +22,7 @@ For GLMs, we have the same setup as before however the response $y_i$ is (usuall For instance, if we have count data $y_i$'s are non-negative integers. GLMs model the parameters through a transformation of the mean that _links_ it to a linear predictor. The function that operationalizes this is _link_ function (usually denoted by $g$). +The link function maps the mean of the distribution to the linear predictor, that is: ```math g(\mathbb{E}(Y_i | \mathbf{x}_i)) = \mathbf{x}_i^\top \mathbb{\beta} = \eta_i . ``` @@ -29,8 +30,11 @@ The link function must be invertible, and the inverse link function is also some ```math \mathbb{E}(Y_i | \mathbf{x}_i) = g^{-1}(\eta_i) . ``` +Another (perhaps more intuitive) way to see this is that the _mean_ function maps the linear predictor to the space of possible mean values. +For instance, when we have binary response the mean must be in the interval $[0,1]$ however, the linear predictor can take any real value. +The mean function (say logistic function, more on that later) takes the linear predictor and squishes into the $(0, 1)$ interval. -Thus, the link function allows us to go from the mean to the associated linear predictor and the mean function allows us to go from the linear predictor to the mean. +In short, the link function allows us to go from the mean to the associated linear predictor and the mean function allows us to go from the linear predictor to the mean. In short, a GLM consists of three key components[^GLMwiki]: @@ -38,6 +42,39 @@ In short, a GLM consists of three key components[^GLMwiki]: 2. The linear predictor $(\eta = \mathbf{x}^\top \boldsymbol{\beta})$ 3. The link function $g(\mu) = \eta$ +### Supported Distribtions + +`GLM.jl` supports models with the following distributions: + +- Normal distribution +- Bernoulli +- Binomial +- Gamma +- Geometric +- Inverse Gaussian +- Negative Binomial +- Poisson + + +### Supported Link functions + +These can be combined with an appropriate link function from the following list: + +- Cauchit link +- complementary log log link +- Identity link +- Inverse Link (or reciprocal) +- Inverse square link +- Logit link +- Log link +- Negative Binomial link +- Power link +- Probit link +- Square root Link + +Note that not all combinations of distribution and link are appropriate. +For instance, + ## Binary Response A binary response $y_i$ can be modeled as: From dba3f7fa8ddbf0fc5190ae14c080193fd97ac773 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 25 Mar 2025 12:35:17 -0500 Subject: [PATCH 4/6] added formula table for links function and corresponding mean functions --- docs/src/models.md | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/docs/src/models.md b/docs/src/models.md index 5c0b317e..78803126 100644 --- a/docs/src/models.md +++ b/docs/src/models.md @@ -56,21 +56,26 @@ In short, a GLM consists of three key components[^GLMwiki]: - Poisson + + + ### Supported Link functions These can be combined with an appropriate link function from the following list: -- Cauchit link -- complementary log log link -- Identity link -- Inverse Link (or reciprocal) -- Inverse square link -- Logit link -- Log link -- Negative Binomial link -- Power link -- Probit link -- Square root Link +| Link | Link function | Mean function | +|:---------- |:---------- |:------------| +| Cauchit link | $\eta = \tan(\pi \times (\mu - \frac{1}{2}))$ | $\mu = \frac{1}{2}+ \frac{\tan^{-1}(\eta)}{\pi}$| +| complimentary log log link | $\eta = \log(-\log(1 - \mu))$ | $\mu = 1 - \exp(-\exp(\eta))$ | +| Identity link | $\eta = \mu$ | $\mu = \eta$ | +| Inverse Link (or reciprocal) | $\eta = \frac{1}{\mu}$ | $\mu = \frac{1}{\eta}$ | +| Inverse square link | $\eta = \frac{1}{\mu^2}$ | $\mu = \frac{1}{\eta^2}$ | +| Logit link | $\eta = \log(\frac{\mu}{1 - \mu})$ | $\mu = (1 + e^{-\eta})^{-1}$ | +| Log link | $\eta = \log(\mu)$ | $\mu = \exp(\eta)$ | +| Negative Binomial link (with parameter $\theta$) | $\eta = \log(\frac{\mu}{\mu + \theta})$ | $\mu = \frac{\theta \exp(\eta)}{1 - \exp(\eta)}$ | +| Power link (with parameter $k$) | $\eta = \mu^{k}$ | $\mu = \eta^{1/k}$ | +| Probit link | $\eta = \Phi^{-1}(\mu)$ | $\mu = \Phi(\eta)$ | +| Square root Link | $\eta = \sqrt{\mu}$ | $\mu = \sqrt{\eta}$ | Note that not all combinations of distribution and link are appropriate. For instance, From a4c46038cbf19d776d73e1e1ec0e4cfc07aa92d9 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 25 Mar 2025 15:21:10 -0500 Subject: [PATCH 5/6] added link to supported link types table to index.md --- docs/src/index.md | 3 ++- docs/src/models.md | 39 ++++++++++++++++++++++----------------- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 65fd104f..fa29d4f3 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -40,7 +40,7 @@ functions are Normal (IdentityLink) Poisson (LogLink) -Currently the available Link types are +Currently, the available Link types are: CauchitLink CloglogLink @@ -54,6 +54,7 @@ Currently the available Link types are ProbitLink SqrtLink +For more details, see [supported link functions](models.md#Supported-Link-functions). Note that the canonical link for negative binomial regression is `NegativeBinomialLink`, but in practice one typically uses `LogLink`. The `NegativeBinomial` distribution belongs to the exponential family only if θ (the shape diff --git a/docs/src/models.md b/docs/src/models.md index 78803126..1c07a070 100644 --- a/docs/src/models.md +++ b/docs/src/models.md @@ -38,16 +38,21 @@ In short, the link function allows us to go from the mean to the associated line In short, a GLM consists of three key components[^GLMwiki]: -1. The distribution -2. The linear predictor $(\eta = \mathbf{x}^\top \boldsymbol{\beta})$ +1. The linear predictor +2. The distribution family 3. The link function $g(\mu) = \eta$ +Different combinations of the family and link define the various GLMs: +```julia +glm(formula, family, link) +``` + ### Supported Distribtions -`GLM.jl` supports models with the following distributions: +`GLM.jl` supports models with following distribution families: -- Normal distribution -- Bernoulli +- Normal distribution (`family = Normal()`) +- Bernoulli (`family = Bernoulli()`) - Binomial - Gamma - Geometric @@ -63,19 +68,19 @@ In short, a GLM consists of three key components[^GLMwiki]: These can be combined with an appropriate link function from the following list: -| Link | Link function | Mean function | +| Link: Constructor | Link function | Mean function | |:---------- |:---------- |:------------| -| Cauchit link | $\eta = \tan(\pi \times (\mu - \frac{1}{2}))$ | $\mu = \frac{1}{2}+ \frac{\tan^{-1}(\eta)}{\pi}$| -| complimentary log log link | $\eta = \log(-\log(1 - \mu))$ | $\mu = 1 - \exp(-\exp(\eta))$ | -| Identity link | $\eta = \mu$ | $\mu = \eta$ | -| Inverse Link (or reciprocal) | $\eta = \frac{1}{\mu}$ | $\mu = \frac{1}{\eta}$ | -| Inverse square link | $\eta = \frac{1}{\mu^2}$ | $\mu = \frac{1}{\eta^2}$ | -| Logit link | $\eta = \log(\frac{\mu}{1 - \mu})$ | $\mu = (1 + e^{-\eta})^{-1}$ | -| Log link | $\eta = \log(\mu)$ | $\mu = \exp(\eta)$ | -| Negative Binomial link (with parameter $\theta$) | $\eta = \log(\frac{\mu}{\mu + \theta})$ | $\mu = \frac{\theta \exp(\eta)}{1 - \exp(\eta)}$ | -| Power link (with parameter $k$) | $\eta = \mu^{k}$ | $\mu = \eta^{1/k}$ | -| Probit link | $\eta = \Phi^{-1}(\mu)$ | $\mu = \Phi(\eta)$ | -| Square root Link | $\eta = \sqrt{\mu}$ | $\mu = \sqrt{\eta}$ | +| Cauchit link : `CauchitLink()` | $\eta = \tan(\pi \times (\mu - \frac{1}{2}))$ | $\mu = \frac{1}{2}+ \frac{\tan^{-1}(\eta)}{\pi}$| +| complimentary log log link : `CloglogLink()` | $\eta = \log(-\log(1 - \mu))$ | $\mu = 1 - \exp(-\exp(\eta))$ | +| Identity link : `IdentityLink()` | $\eta = \mu$ | $\mu = \eta$ | +| Inverse Link (or reciprocal) : `InverseLink()` | $\eta = \frac{1}{\mu}$ | $\mu = \frac{1}{\eta}$ | +| Inverse square link : `InverseSquareLink()` | $\eta = \frac{1}{\mu^2}$ | $\mu = \frac{1}{\eta^2}$ | +| Logit link : `LogitLink()` | $\eta = \log(\frac{\mu}{1 - \mu})$ | $\mu = (1 + e^{-\eta})^{-1}$ | +| Log link : `LogLink()` | $\eta = \log(\mu)$ | $\mu = \exp(\eta)$ | +| Negative Binomial link: `NegativeBinomialLink(θ)` | $\eta = \log(\frac{\mu}{\mu + \theta})$ | $\mu = \frac{\theta \exp(\eta)}{1 - \exp(\eta)}$ | +| Power link: `PowerLink(k)`| $\eta = \mu^{k}$ | $\mu = \eta^{1/k}$ | +| Probit link: `ProbitLink()` | $\eta = \Phi^{-1}(\mu)$ | $\mu = \Phi(\eta)$ | +| Square root Link: `SqrtLink()` |$\eta = \sqrt{\mu}$ | $\mu = \sqrt{\eta}$ | Note that not all combinations of distribution and link are appropriate. For instance, From 3ef731d59f9678c43a2b8fee82b00c15d9dce413 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Tue, 8 Apr 2025 09:52:14 -0500 Subject: [PATCH 6/6] added negative binomial to models.md and other stuff --- docs/src/models.md | 123 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 111 insertions(+), 12 deletions(-) diff --git a/docs/src/models.md b/docs/src/models.md index 1c07a070..25fc441d 100644 --- a/docs/src/models.md +++ b/docs/src/models.md @@ -70,20 +70,32 @@ These can be combined with an appropriate link function from the following list: | Link: Constructor | Link function | Mean function | |:---------- |:---------- |:------------| -| Cauchit link : `CauchitLink()` | $\eta = \tan(\pi \times (\mu - \frac{1}{2}))$ | $\mu = \frac{1}{2}+ \frac{\tan^{-1}(\eta)}{\pi}$| -| complimentary log log link : `CloglogLink()` | $\eta = \log(-\log(1 - \mu))$ | $\mu = 1 - \exp(-\exp(\eta))$ | -| Identity link : `IdentityLink()` | $\eta = \mu$ | $\mu = \eta$ | +| | $\eta = g(\mu)$ | $\mu = h(\eta) = g^{-1}(\eta)$ | +| Cauchit link: `CauchitLink()` | $\eta = \tan(\pi \times (\mu - \frac{1}{2}))$ | $\mu = \frac{1}{2}+ \frac{\tan^{-1}(\eta)}{\pi}$| +| complimentary log log link: `CloglogLink()` | $\eta = \log(-\log(1 - \mu))$ | $\mu = 1 - \exp(-\exp(\eta))$ | +| Identity link: `IdentityLink()` | $\eta = \mu$ | $\mu = \eta$ | | Inverse Link (or reciprocal) : `InverseLink()` | $\eta = \frac{1}{\mu}$ | $\mu = \frac{1}{\eta}$ | -| Inverse square link : `InverseSquareLink()` | $\eta = \frac{1}{\mu^2}$ | $\mu = \frac{1}{\eta^2}$ | -| Logit link : `LogitLink()` | $\eta = \log(\frac{\mu}{1 - \mu})$ | $\mu = (1 + e^{-\eta})^{-1}$ | -| Log link : `LogLink()` | $\eta = \log(\mu)$ | $\mu = \exp(\eta)$ | +| Inverse square link: `InverseSquareLink()` | $\eta = \frac{1}{\mu^2}$ | $\mu = \frac{1}{\eta^2}$ | +| Logit link: `LogitLink()` | $\eta = \log(\frac{\mu}{1 - \mu})$ | $\mu = (1 + e^{-\eta})^{-1}$ | +| Log link: `LogLink()` | $\eta = \log(\mu)$ | $\mu = \exp(\eta)$ | | Negative Binomial link: `NegativeBinomialLink(θ)` | $\eta = \log(\frac{\mu}{\mu + \theta})$ | $\mu = \frac{\theta \exp(\eta)}{1 - \exp(\eta)}$ | | Power link: `PowerLink(k)`| $\eta = \mu^{k}$ | $\mu = \eta^{1/k}$ | | Probit link: `ProbitLink()` | $\eta = \Phi^{-1}(\mu)$ | $\mu = \Phi(\eta)$ | | Square root Link: `SqrtLink()` |$\eta = \sqrt{\mu}$ | $\mu = \sqrt{\eta}$ | Note that not all combinations of distribution and link are appropriate. -For instance, +The following table summarizes the combinations + + +| Link\Dist | Normal | Bern. | Binom. | Gamma | Geom. | Inv-Gaussian | Neg-Binom | Pois. | +|:--|:------:|:----------:|:---------:|:------:|:----------:|:-----------------:|:--------:|:------:| +| Identity | `lm` | ❌ | ❌ | ? | ? | ? | ❌ | ❌ | +| Cauchit | ? | ✓ | ✓ | ? | ? | ? | ? | ? | +| Cloglog | ? | ✓ | ✓ | ? | ? | ? | ? | ? | +| Inverse | ? | ? | ? | ✅ | ? | ? | ? | ? | +| Logit | ? | ✅ | ✓ | ? | ? | ? | ? | ? | +| Log | ? | ? | ? | ? | ? | ? | ✓ | ✅ | +| NegBinomial | ? | ? | ? | ? | ? | ? | ✅ | ✓ | ## Binary Response @@ -144,12 +156,12 @@ mod = glm(@formula(...), data, Bernoulli(), ProbitLink()) ### Poisson Regression for count data -When we have count data, i.e. the observed $y_i$'s are non-negative integers we can use the following model: -$$ - ( Y_i\ \vert\ \mathbf{X}_i = mathbf{x}_i) \sim operatorname{Poisson}(\lambda_i). -$$ +When we have count data, i.e. the observed $y_i$'s are non-negative integers we can use the following model: +```math + ( Y_i\ \vert\ \mathbf{X}_i = \mathbf{x}_i) \sim \operatorname{Poisson}(\lambda_i). +``` -Since $\mathb{E}(Y_i\ \vert\ mathbf{x}_i) = \lambda_i$ this can be formulated as a GLM with _log_ as the link function: +Since $\mathbb{E}(Y_i\ \vert\ \mathbf{x}_i) = \lambda_i$ this can be formulated as a GLM with _log_ as the link function: ```math \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} . ``` @@ -163,6 +175,93 @@ mod = glm(@formula(...), data, Poisson(), LogLink()) ``` +### Binomial Regression + +**Note:** When performing binomial regression, the response variable must be a proportion, i.e. $0 \leq y_i \leq 1$. +If the data contains the number of successes and the number of trials, first generate the success proportion, and use that as the response variable. + +We want to model the number of successes from $n_i$ trials using a Binomial distribution as: +```math + Y_i \sim \operatorname{Binomial}(n_i, p_i) . +``` +However, we work with the proportion of successes and model the mean parameter $p_i$ which is: +```math +p_i = \mathbb{E}\left(\frac{Y_i}{n_i} \right) , +``` +using the link function: +```math + \log\left(\frac{p_i}{1 - p_i}\right) = \mathbf{x}_i^\top \boldsymbol{\beta} , +``` +and the corresponding mean function is the _logistic_ function $p_i = \sigma(\mathbf{x}_i^\top \boldsymbol{\beta})$. + +### Geometric GLM + +We wish to model count responses $Y_i \in \{0, 1, 2, ... \}$ using the Geometric distribution, that is: +```math + Y_i \sim \operatorname{Geometric}(p_i), +``` +which has the mean $\mu_i = \frac{1 - p_i}{p_i}$. +The canonical link function for Geometric GLM is the _log_ link: +```math + log(\mu_i) =\mathbf{x}_i^\top \boldsymbol{\beta}. +``` + + + +## Negative Binomial Models + +There are many parameterizations of the negative binomial distribution, but we use the following: +```math +p(Y = y) = \frac{\Gamma(y + \theta)}{ \Gamma(\theta) \times y!} \times \frac{\theta^\theta \mu^y}{(\mu + \theta)^{(y + \theta)}} , +``` +where $\mu$ is the mean of the distribution and $\theta$ is the number of successes when the experiment is stopped. +This parameterization can easily be derived by noting that $\mu = \frac{\theta(1 - p)}{p}$, and plugging in the corresponding value of $p$ in the more commonly used parameterization detailed on [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). +Also note that we use the $\Gamma$ function instead of factorial here, because we will allow the $\theta$ parameter to take noninteger positive values. + +We will model the mean as a function of the predictor variables. +The canonical link function for this distribution is the negative binomial link function that is: +```math + \log \left(\frac{\mu_i}{\mu_i + \theta}\right) = \mathbf{x}_i^{\top} \boldsymbol{\beta} +``` +but the log-link function is frequently used: +```math + \log(\mu_i) = \mathbf{x}_i^{\top} \boldsymbol{\beta} +``` + +Notice that the canonical link function takes $\theta$ as an argument. +This is because the negative binomial distribution belongs to the exponential family only if the number of successes parameter ($\theta$) is known. +As a result, a model that attempts to estimate both $\theta$ and $\boldsymbol{\beta}$ simultaneously is technically not a GLM. +However, if the $\theta$ parameter is specified beforehand the model is a GLM. + +This package supports both versions of the model, that is (1) the GLM version with the $\theta$ parameter is assumed to be known, and (2) the more general model where $\theta$ is estimated along with the $\beta$ parameter. + +## Negative Binomial Model with known $\theta$ + +This model can be fit as: +```julia +glm(frm, dt, NegativeBinomial(θ)) +``` +or by supplying the log-link: +```julia +glm(frm, dt, NegativeBinomial(θ), LogLink()) +``` + +## Negative binomial Model with Simultaneous Estimation of $\theta$ + +While this is technically not a GLM, the package does support this version using the `negbin` function which can be called as follows: +```julia +negbin(frm, dt) +``` +and optionally using a link function (usually the log-link) as: +```julia +negbin(frm, dt, link) +``` + +This function estimates both the coefficient vector $\boldsymbol{\beta}$ and the parameter $\theta$. +This is done in an iterative procedure in which the coefficient vector is estimated using a GLM assuming $\theta$ is fixed (and set to its previous estimate), and the $\theta$ is estimated by maximizing the log-likelihood assuming the coefficient vector is fixed (set to its previous estimate). +The first step in each iteration internally calls `glm`. + + [^logit]: https://en.wikipedia.org/wiki/Logit [^logistic]: https://en.wikipedia.org/wiki/Logistic_function [^GLMwiki]: https://en.wikipedia.org/wiki/Generalized_linear_model#Model_components