Generalized linear models and Fisher scoring

Fix a finite integer $n\ge 1$ and a dispersion value $\varphi$. For each $1\le i\le n$, let ${\mathcal{Y}}_i\subseteq \mathbb{R}$ carry its Borel sigma algebra and let ${\nu }_i$ be a sigma-finite positive measure on it. Let $J\subseteq \mathbb{R}$ be an open interval, let $b:J\to \mathbb{R}$ be twice continuously differentiable, let $a_i(\varphi )>0$, and let $c_i(\cdot ,\varphi )$ be measurable. An explicit one-parameter exponential-family response contract supplies, for each ${\theta }_i\in J$, a density or mass function relative to ${\nu }_i$,

$f_i(y_i;{\theta }_i,\varphi )=\exp \text{\hspace{0.17em}⁣}\{\frac{y_i{\theta }_i-b({\theta }_i)}{a_i(\varphi )}+c_i(y_i,\varphi )\}\mathrm{.}$

The contract requires $f_i\ge 0$, $\int f_i(\cdot ;{\theta }_i,\varphi )d{\nu }_i=1$ for every ${\theta }_i\in J$, and a fixed support $S_i\in {\mathcal{A}}_i$ such that $f_i>0$ on $S_i$ and $f_i=0$ off $S_i$, up to a ${\nu }_i$-null set independent of ${\theta }_i$. When the mean and variance identities are used, it also requires a finite second moment and locally uniform integrable envelopes for the first and second ${\theta }_i$ derivatives of $f_i$ on $S_i$.

For a vector $\theta =({\theta }_1,\dots ,{\theta }_n)\in J^n$, the product sampling law is the finite product ${⨂}_{i=1}^n{P}_{i,{\theta }_i,\varphi }$, where $P_{i,{\theta }_i,\varphi }$ has density $f_i$ relative to ${\nu }_i$. The coordinate responses are independent in the sense of Independent random elements. Their product base measure is $\nu ={⨂}_{i=1}^n{\nu }_i$ on ${⨂}_{i=1}^n{\mathcal{A}}_i$. The next lemma verifies that the joint density is ${\prod }_i{f}_i$.

The design or conditioning convention, observation units, and response/index roles are part of the surrounding scientific model contract. A sampling-law specification is not evidence that the scientific model is empirically adequate.

Every scientific application also records a Scientific statistical model contracts: it names the observation space and record, scientific and analysis units, index meanings, response and predictor maps, target, observation mechanism, support, law, conditioning, and adequacy assumptions. The contract does not itself imply independence or common support; those are explicit hypotheses of this response-family definition. Here $i$ is only a component index until the application assigns it a scientific unit or time meaning.

In the setup of Explicit exponential-family response contracts, let $P_{\theta }={⨂}_{i=1}^n{P}_{i,{\theta }_i,\varphi }$ and $\nu ={⨂}_{i=1}^n{\nu }_i$. Then the joint sampling law has density

$q_{\theta }(y_1,\dots ,y_n)={\prod }_{i=1}^n{f}_i(y_i;{\theta }_i,\varphi )$

relative to $\nu$. At any $y$ for which every factor $f_i(y_i;{\theta }_i,\varphi )$ is positive, the log likelihood is $\log q_{\theta }(y)={\sum }_i\log f_i(y_i;{\theta }_i,\varphi )$.

Fix $i$ and $\varphi$ in Explicit exponential-family response contracts, and write $a=a_i(\varphi )$. Under its fixed-support, twice differentiable density, local integrable derivative-envelope, and finite-second-moment hypotheses,

${\mathbb{E}}_{{\theta }_i,\varphi }[Y_i]=b^{\mathrm{\prime }}({\theta }_i),{\mathrm{Var}}_{{\theta }_i,\varphi }(Y_i)=a{b}^{\mathrm{\prime }\mathrm{\prime }}({\theta }_i)\mathrm{.}$

Fix $n\ge 1$, $p\ge 1$, a finite design matrix $X\in {\mathbb{R}}^{n\times p}$ with rows $x_i^{\mathsf{T}}$, a dispersion value $\varphi$, and an open parameter domain $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$. A generalized linear model (GLM) on this design consists of the response family in Explicit exponential-family response contracts, a differentiable one-to-one link $g:\mathcal{M}\to g(\mathcal{M})$ on a specified mean domain $\mathcal{M}$, and a specified map $\tau :\mathcal{M}\to J$ satisfying $b^{\mathrm{\prime }}(\tau (\mu ))=\mu$. For $\beta \in \mathrm{\Theta }$, put

${\eta }_i=x_i^{\mathsf{T}}\beta ,{\mu }_i=g^{-1}({\eta }_i),{\theta }_i=\tau ({\mu }_i),$

and require ${\eta }_i\in g(\mathcal{M})$ and ${\theta }_i\in J$ for every $i$. The canonical link is the choice $g=\tau$, equivalently ${\theta }_i={\eta }_i$. When ${\theta }_i={\eta }_i$ and the component scales are common, write $a_i(\varphi )=a(\varphi )>0$; this common-scale restriction is part of every formula that uses $a(\varphi )$.

The joint sampling density is the product density from The product response law has the product density. For an observed $y=(y_1,\dots ,y_n)$ with $f_i(y_i;{\theta }_i,\varphi )>0$ for every $i$, define its log likelihood by

$\mathrm{\ell }(\beta ;y,\varphi )={\sum }_{i=1}^n\{\frac{y_i{\theta }_i-b({\theta }_i)}{a_i(\varphi )}+c_i(y_i,\varphi )\}\mathrm{.}$

For a differentiable log likelihood, the score is $U(\beta ;y,\varphi ):={\mathrm{\nabla }}_{\beta }\mathrm{\ell }(\beta ;y,\varphi )$. When its entries have finite second moments, the expected score outer-product information is $I(\beta ,\varphi ):={\mathbb{E}}_{\beta ,\varphi }[U{U}^{\mathsf{T}}]$. This page treats $X$ as fixed or conditions on its realized value; it supplies no random-design asymptotic conclusion.

Let $(Y_i{)}_{i=1}^n$ follow Generalized linear model with explicit conditioning with a common scale $a(\varphi )>0$, canonical link ${\theta }_i={\eta }_i=x_i^{\mathsf{T}}\beta$, and twice differentiable $b$. For an observed $y$ with $f_i(y_i;{\eta }_i,\varphi )>0$ for every $i$,

${\mathrm{\nabla }}_{\beta }\mathrm{\ell }(\beta ;y,\varphi )=\frac{X^{\mathsf{T}}(y-\mu )}{a(\varphi )},{\mathrm{\nabla }}_{\beta }^2\mathrm{\ell }(\beta ;y,\varphi )=-\frac{X^{\mathsf{T}}\mathrm{diag}(b^{\mathrm{\prime }\mathrm{\prime }}({\eta }_i))X}{a(\varphi )}\mathrm{.}$

Consider Generalized linear model with explicit conditioning at an admissible finite $\beta$ with common scale $a(\varphi )>0$. Assume each ${\eta }_i\mapsto {\theta }_i$ is differentiable, each ${\eta }_i\mapsto {\mu }_i$ is differentiable with ${\dot{\mu }}_i:=d{\mu }_i\mathrm{/}d{\eta }_i\mathrm{\ne }0$, $V_i:=b^{\mathrm{\prime }\mathrm{\prime }}({\theta }_i)>0$, the component scores have finite second moments, and the common-support, differentiability, and local integrable derivative-envelope conditions from Explicit exponential-family response contracts hold. Put

$w_i=\frac{{\dot{\mu }}_i^2}{a(\varphi )V_i},z_i={\eta }_i+\frac{y_i-{\mu }_i}{{\dot{\mu }}_i}\mathrm{.}$

For every observed $y$ with $f_i(y_i;{\theta }_i,\varphi )>0$ for each $i$, the score and expected score outer-product information are

$U(\beta ;y,\varphi )=X^{\mathsf{T}}{\left(\frac{{\dot{\mu }}_i(y_i-{\mu }_i)}{a(\varphi )V_i}\right)}_{i=1}^n,I(\beta ,\varphi )=X^{\mathsf{T}}\mathrm{diag}(w_i)X\mathrm{.}$

In particular, if $X$ has full column rank, then $I$ is positive definite.

Fix an observed response vector $y$ for which the component densities are positive throughout both models' parameter sets, and two stated GLMs with the same response family, base support, and log-likelihood function $\mathrm{\ell }(\beta ;y)$. They are nested here only when a stated embedding identifies the fitted parameter set with a subset ${\mathrm{\Theta }}_{\mathrm{f}\mathrm{i}\mathrm{t}}\subseteq {\mathrm{\Theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ of the saturated reference parameter set. Require both sets to be nonempty and both suprema to be finite real numbers. The likelihood deviance is

$D(y)=2\{{\sup }_{\beta \in {\mathrm{\Theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}}}\mathrm{\ell }(\beta ;y)-{\sup }_{\beta \in {\mathrm{\Theta }}_{\mathrm{f}\mathrm{i}\mathrm{t}}}\mathrm{\ell }(\beta ;y)\}\mathrm{.}$

The definition makes no claim about deviance residuals, goodness-of-fit diagnostics, or a reference distribution. The two full-text GLM notes cited above support likelihood setup but do not supply a deviance-diagnostics treatment; the nonnegativity result below is proved only from the stated set inclusion.

Suppose the canonical-link GLM has an open convex parameter domain $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$, full-column-rank design $X$, common scale $a(\varphi )>0$, and $b^{\mathrm{\prime }\mathrm{\prime }}({\eta }_i)>0$ at every predictor vector in $X\mathrm{\Theta }$. For a fixed response vector $y$ whose component densities are positive throughout $\mathrm{\Theta }$, $\mathrm{\ell }(\cdot ;y,\varphi )$ is strictly concave on $\mathrm{\Theta }$. It has at most one stationary point; if a stationary point exists, it is the unique global maximizer. This conclusion does not assert that a finite maximizer exists.

Inputs are an observed response vector $y$ with $f_i(y_i;{\theta }_i,\varphi )>0$ for every $i$, a finite design matrix $X$, an admissible finite current parameter $\beta$, and a GLM with common scale $a(\varphi )>0$. Require each ${\dot{\mu }}_i=d{\mu }_i\mathrm{/}d{\eta }_i$ to be finite and nonzero and each $V_i=b^{\mathrm{\prime }\mathrm{\prime }}({\theta }_i)$ to be finite and positive.

1. Compute ${\eta }_i=x_i^{\mathsf{T}}\beta$ and ${\mu }_i=g^{-1}({\eta }_i)$.
2. Compute the working response and positive weight $z_i={\eta }_i+\frac{y_i-{\mu }_i}{{\dot{\mu }}_i},w_i=\frac{{\dot{\mu }}_i^2}{a(\varphi )V_i}\mathrm{.}$
3. Form $W^{1\mathrm{/}2}=\mathrm{diag}(\sqrt{w_i})$ using Square roots exist: a unique $\sqrt{a}\ge 0$ with $(\sqrt{a}{)}^2=a$; the positives are $\{x^2:x\mathrm{\ne }0\}$, and put $A=W^{1\mathrm{/}2}X$, $r=W^{1\mathrm{/}2}z$. By The Moore--Penrose pseudoinverse $A^{+}$ as the solution of the four Penrose equations and Every finite real or complex matrix has a unique Moore--Penrose pseudoinverse, form $A^{+}$. Return ${\beta }_{\mathrm{n}\mathrm{e}\mathrm{w}}=A^{+}r$, the unique minimum-Euclidean-norm member of the nonempty weighted least-squares minimizer set by For every right-hand side $b$, $A^{+}b$ is the unique least-squares solution of minimum Euclidean norm.

The Moore--Penrose convention makes the procedure single-valued even when $X$ is rank deficient; the underlying least-squares minimizer need not then be unique. Existence and the minimum-norm property follow from Every finite real or complex matrix has a unique Moore--Penrose pseudoinverse and For every right-hand side $b$, $A^{+}b$ is the unique least-squares solution of minimum Euclidean norm. The output need not lie in $\mathrm{\Theta }$, and this definition makes no claim that repeated updates terminate or converge to a maximizer.

For every input satisfying the Definition of One minimum-norm IRLS step, every intermediate quantity in its three instructions is uniquely determined and finite. Its output ${\beta }_{\mathrm{n}\mathrm{e}\mathrm{w}}\in {\mathbb{R}}^p$ is the unique minimum-Euclidean-norm member of the nonempty set ${\mathrm{argmin}}_{\gamma \in {\mathbb{R}}^p}{\sum }_{i=1}^n{w}_i(z_i-x_i^{\mathsf{T}}\gamma {)}^2\mathrm{.}$ This holds for every rank of $X$, including $X=0$. The weighted normal matrix $H=X^{\mathsf{T}}WX$ is positive semidefinite and has kernel $\ker H=\ker X$; its singularity does not obstruct this output. No assertion that the output lies in the model parameter domain, or that repeated updates converge, is made.

This is the explicit well-definedness discharge for the local algorithm. The proof uses the complete locally proved square-root, pseudoinverse-existence, and minimum-norm least-squares suppliers named above. The weighted objective is checked directly, and no Fisher-scoring equivalence or convergence theorem is assumed to prove existence of the procedure's output.

For the deviance in Likelihood deviance for nested models on common support, $D(y)\ge 0$.

Under the admissibility hypotheses of One minimum-norm IRLS step, assume also that each map ${\eta }_i\mapsto {\theta }_i$ is differentiable, each map ${\eta }_i\mapsto {\mu }_i$ is differentiable with finite nonzero derivative ${\dot{\mu }}_i$, and $V_i:=b^{\mathrm{\prime }\mathrm{\prime }}({\theta }_i)>0$. Each response has finite second moment, its density has the fixed common support and locally uniform integrable envelopes for its first and second ${\theta }_i$ derivatives required by Explicit exponential-family response contracts, and each component score has finite second moment. If $I(\beta )=X^{\mathsf{T}}WX$ is invertible, let $U(\beta )$ and $I(\beta )$ be the score and expected score outer-product information from General-link score and Fisher information. Then the weighted least-squares output satisfies

${\beta }_{\mathrm{n}\mathrm{e}\mathrm{w}}=\beta +I(\beta {)}^{-1}U(\beta )\mathrm{.}$