Regression Analysis Note

Matrix Theory

Matrix Theory Note

Matrix Theory

Basic Theories

Matrix

Definition

$A_{m\times n}=\left[\begin{array}{ccc}{a}_{11}& \dots & a_{1n}\\ ⋮& \ddots & ⋮\\ a_{m1}& \dots & a_{mn}\end{array}\right]$

Operations

Addition, Subtraction

$A_{m\times n}\pm B_{m\times n}=\left[\begin{array}{ccc}{a}_{11}\pm b_{11}& \dots & a_{1n}\pm b_{1n}\\ ⋮& \ddots & ⋮\\ a_{m1}\pm b_{m1}& \dots & a_{mn}\pm b_{mn}\end{array}\right]$

Scalar Multiplication

$k{A}_{m\times n}=\left[\begin{array}{ccc}k{a}_{11}& \dots & k{a}_{1n}\\ ⋮& \ddots & ⋮\\ k{a}_{m1}& \dots & k{a}_{mn}\end{array}\right]$

Transpose

$A^{⊺}$ $a_{ij}\stackrel{⊺}{\to }{a}_{ji}$, $A_{m\times n}\stackrel{\top }{\to }{A}_{n\times m}$

Trace

Determinant

Matrix Multiplication

$A_{m\times p}\cdot B_{p\times n}=\left[\begin{array}{ccc}{a}_{11}& \dots & a_{1p}\\ ⋮& \ddots & ⋮\\ a_{m1}& \dots & a_{mp}\end{array}\right]\left[\begin{array}{ccc}{b}_{11}& \dots & b_{1n}\\ ⋮& \ddots & ⋮\\ b_{p1}& \dots & b_{pn}\end{array}\right]={\left[\begin{array}{ccc}{a}_{11}+b_{11}+\cdots +a_{1p}+b_{p1}& \dots & a_{11}+b_{1n}+\cdots +a_{1p}+b_{pn}\\ ⋮& \ddots & ⋮\\ a_{m1}+b_{11}+\cdots +a_{mp}+b_{p1}& \dots & a_{m1}+b_{1n}+\cdots +a_{mp}+b_{pn}\end{array}\right]}_{m\times n}$ $A_i^{⊺}\cdot B_j={\sum }_{k=1}^p{a}_{ik}\cdot b_{kj}=c_{ij}$ where $1\le i\le m$, $1\le j\le n$

Vector-Matrix Multiplication

$A\mathbf{x}=[A_1,A_2,\dots ,A_n][x_1,x_2,\dots ,x_n{]}^{⊺}=x_1{A}_1+x_2{A}_2+\cdots +x_n{A}_n={\sum }_{k=1}^n{x}_k{A}_k$

Link to original

Trace

Definition

$\mathrm{tr}(A)={\sum }_{i=1}^n{a}_{ii}=a_{11}+a_{22}+\cdots +a_{nn}$

The sum of elements on the main diagonal

Facts

$\mathrm{tr}(A+B)=\mathrm{tr}(A)+\mathrm{tr}(B)$

$\mathrm{tr}(cA)=c\mathrm{tr}(A)$

$\mathrm{tr}(A^{⊺})=\mathrm{tr}(A)$

$\mathrm{tr}(AB)=\mathrm{tr}(BA)$

$\mathrm{tr}(A{B}^{⊺})=\mathrm{tr}(B{A}^{⊺})$

$\mathrm{tr}(ABCD)=\mathrm{tr}(BCDA)=\mathrm{tr}(CDAB)=\mathrm{tr}(DABC)$^[Cyclic property]

$\mathrm{tr}(\mathbf{A})={\sum }_{i=1}^n{\lambda }_i$ The Trace of the matrix is equal to the sum of the eigenvalues of the matrix.

Proof Since the matrix $\mathbf{A}-\lambda \mathbf{I}$ only has $\lambda$ term on diagonal, and the calculation of cofactor deletes a row and a column, the coefficient of ${\lambda }^{n-1}$ of $\text{det}(\mathbf{A}-\lambda \mathbf{I})$ is always $-(a_{11}+a_{22}+\cdots +a_{nn})$. Also, since ${\lambda }_i$ are the solution of the Characteristic Polynomial, the expression is factorized as $\text{det}(\mathbf{A}-\lambda \mathbf{I})=(\lambda -{\lambda }_1)(\lambda -{\lambda }_2)\dots (\lambda -{\lambda }_n)$ and the coefficient of ${\lambda }^{n-1}$ become a $-({\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n)$ Therefore, $-(a_{11}+a_{22}+\cdots +a_{nn})=-({\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n)$ and $\text{Trace}(\mathbf{A})={\sum }_{i=1}^n{\lambda }_i$.

Link to original

Link to original

Determinant

Definition

$\text{det}(\mathbf{A})=|\mathbf{A}|$

Computation

$2\times 2$ matrix

$\left|\begin{array}{cc}a& b\\ c& d\end{array}\right|=ad-bc$

$n\times n$ matrix

Laplace Expansion

Definition

Let $C$ be a cofactor matrix of $\mathbf{A}$

Along the $j$th column gives

$\text{det}(\mathbf{A})={\sum }_{i=1}^n{a}_{ij}{C}_{ij}={\sum }_{i=1}^n{a}_{ij}(-1{)}^{i+j}{M}_{i,j}$

Along the $i$th row gives

$\text{det}(\mathbf{A})={\sum }_{j=1}^n{a}_{ij}{C}_{ij}={\sum }_{j=1}^n{a}_{ij}(-1{)}^{i+j}{M}_{i,j}$

Link to original

Facts

$\text{det}\mathbf{A}\ne 0\iff \exists {\mathbf{A}}^{-1}$ $\text{det}\mathbf{A}=0\iff \nexists {\mathbf{A}}^{-1}$

The volume of box^[parallelepiped] $\in {\mathbb{R}}^n$ is expressed by a determinant

$\text{det}(\mathbf{I})=1$

$\text{det}(\mathbf{A})$ changes sign, when two rows are exchanged (permutated)

$\text{det}(\mathbf{A})$ depends on linearly on the 1st row.

$\text{det}(t\mathbf{A})=t^n\text{det}(\mathbf{A})$

$\text{det}(\mathbf{AB})=\text{det}(\mathbf{A})\text{det}(\mathbf{B})$

If two row or column vectors in $\mathbf{A}$ are equal, then $\text{det}(\mathbf{A})=0$

Gauss Elimination does not change $\text{det}(\mathbf{A})$

If a matrix $\mathbf{A}$ has zero row or column vectors, then $\text{det}(\mathbf{A})=0$

If a matrix $\mathbf{A}$ is diagonal or triangular, then the determinant of $\mathbf{A}$ is the product of diagonal elements $\text{det}(\mathbf{A})={\prod }_{i=1}^n{a}_{ii}$ where $a_{ii}$ is the $i,i$ element of the matrix $\mathbf{A}$.

$\text{det}({\mathbf{A}}^{-1})=\text{det}(\mathbf{A}{)}^{-1}$

$\text{det}({\mathbf{A}}^{⊺})=\text{det}(\mathbf{A})$

$|\mathbf{A}|={\prod }_{i=1}^n{\lambda }_i$ The Determinant of the matrix is equal to the product of the eigenvalues of the matrix.

Link to original

$|I_m+AB|=|I_n+BA|$ where the dimensions of the matrices are $m\times n$ and $n\times m$ respectively.

$|A\pm \mathbf{u}{\mathbf{u}}^{⊺}|=|A|(1\pm {\mathbf{u}}^{⊺}{A}^{-1}\mathbf{u})$

Link to original

Idempotent Matrix

Definition

$A^2=A$ A matrix whose squared matrix is itself.

Facts

Eigenvalues of idempotent matrix is either $0$ or $1$.

Let $A$ be a symmetric matrix. Then, $A$ has $r$ eigenvalues equal to $1$ and the rest zero $\iff A^2=A,\mathrm{rank}(A)=r$.

Let $A$ be an idempotent matrix, then $\mathrm{rank}(A)=\mathrm{tr}(A)$

If $A$ is idempotent, then $I-A$ also idempotent.

Link to original

Inverse Matrix

Rank of Matrix

Definition

$\mathrm{rank}(\mathbf{A})$

The rank of matrix is the number of linearly independent column vectors, or the number of non-zero pivots in Gauss Elimination

Facts

$\mathrm{rank}(AB)\le \min (\mathrm{rank}(A),\mathrm{rank}(B))$

For the non-singular matrices $P$ and $Q$, $\mathrm{rank}(PAQ)=\mathrm{rank}(A)$

Rank–Nullity Theorem

$\mathrm{rank}(A)=\mathrm{rank}(A^{⊺})=\mathrm{rank}(A^{⊺}A)=\mathrm{rank}(A{A}^{⊺})$

If a matrix $A$ is Symmetric Matrix, then $\mathrm{rank}(A)$ is the number of non-zero eigenvalues.

Link to original

Inverse Matrix

Definition

${\mathbf{A}}^{-1}$ satisfies $\mathbf{A}{\mathbf{A}}^{-1}={\mathbf{A}}^{-1}\mathbf{A}=\mathbf{I}$ for given matrix $\mathbf{A}$

Computation

2 x 2 matrix

${\mathbf{A}}^{-1}=\frac{1}{|\mathbf{A}|}\left[\begin{array}{cc}d& -b\\ -c& a\end{array}\right]$ where the matrix $\mathbf{A}=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$, and $|\mathbf{A}|$ is a Determinant of the matrix

Using Cofactor

$\mathbf{A}\mathrm{adj}(\mathbf{A})={\mathbf{AC}}^{⊺}=\text{det}(\mathbf{A})\mathbf{I}$ where the matrix $\mathbf{C}=\left[\begin{array}{cccc}{C}_{11}& C_{12}& \cdots & C_{1n}\\ C_{21}& C_{22}& \cdots & C_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ C_{n1}& C_{n2}& \cdots & C_{nn}\end{array}\right]$ is the cofactor matrix, which is formed by all the cofactors of a given matrix $\mathbf{A}$. Then, by the definition of inverse matrix, ${\mathbf{A}}^{-1}=\frac{1}{\text{det}(\mathbf{A})}{\mathbf{C}}^{⊺}$

Link to original

Sherman–Morrison Formula

Definition

Suppose $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is an invertible matrix an d $\mathbf{v},\mathbf{u}\in {\mathbb{R}}^n$. Then, $\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺}$ is invertible $\iff 1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}\ne 0$. In this case, $(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺}{)}^{-1}={\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}}$

Proof

Multiplying $(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺})$ to the RHS gives $\left({\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}}\right)(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺})=\mathbf{I}+{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺})}{1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}}\text{(1)}$

Since ${\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}$ is a scalar, The numerator of the last term can be expressed as ${\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺})={\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}+({\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}){\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}=(1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}){\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}$

So, $\text{(1)}=\mathbf{I}+{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}-\frac{(1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}){\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}}{1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}}=\text{(1)}=\mathbf{I}+{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}-{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}=\mathbf{I}$ $\therefore (\mathbf{A}+\mathbf{u}{\mathbf{v}}^{⊺}{)}^{-1}={\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^{⊺}{\mathbf{A}}^{-1}\mathbf{u}}$

Examples

Updating Fitted Least Square Estimator

Sherman–Morrison formula can be used for updating fitted Least Square estimator

Let $\hat{\mathbf{x}}=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{b}$ be the least square estimator of $\mathbf{Ax}=\mathbf{b}$ and the matrices with a new data be ${\mathbf{A}}_a=\left[\begin{array}{c}\mathbf{A}\\ {\mathbf{a}}^{⊺}\end{array}\right]$ and ${\mathbf{b}}_a=\left[\begin{array}{c}\mathbf{b}\\ b\end{array}\right]$ Then, ${\hat{\mathbf{x}}}_{\text{new}}=({\mathbf{A}}_a^{⊺}{\mathbf{A}}_a{)}^{-1}{\mathbf{A}}_a^{⊺}{\mathbf{b}}_a={\left(\left[\begin{array}{cc}{\mathbf{A}}^{⊺}& \mathbf{a}\end{array}\right]\left[\begin{array}{c}\mathbf{A}\\ {\mathbf{a}}^{⊺}\end{array}\right]\right)}^{-1}\left[\begin{array}{cc}{\mathbf{A}}^{⊺}& \mathbf{a}\end{array}\right]\left[\begin{array}{c}\mathbf{b}\\ b\end{array}\right]=({\mathbf{A}}^{⊺}\mathbf{A}+\mathbf{a}{\mathbf{a}}^{⊺}{)}^{-1}({\mathbf{A}}^{⊺}\mathbf{b}+b\mathbf{a})$

Let $\mathbf{P}:=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}$ for convenience Then, ${\mathbf{P}}_a:=({\mathbf{A}}^{⊺}\mathbf{A}+\mathbf{a}{\mathbf{a}}^{⊺}{)}^{-1}=\mathbf{P}-\frac{\mathbf{Pa}{\mathbf{a}}^{⊺}\mathbf{P}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}-\frac{({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}\mathbf{a}{\mathbf{a}}^{⊺}({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}}{1+{\mathbf{a}}^{⊺}({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}\mathbf{a}}$ by Sherman-Morrison formula

So, ${\hat{\mathbf{x}}}_{\text{new}}=({\mathbf{A}}^{⊺}\mathbf{A}+\mathbf{a}{\mathbf{a}}^{⊺}{)}^{-1}({\mathbf{A}}^{⊺}\mathbf{b}+b\mathbf{a})=({\mathbf{A}}^{⊺}\mathbf{A}{)}^{-1}{\mathbf{A}}^{⊺}\mathbf{b}-\frac{\mathbf{Pa}{\mathbf{a}}^{⊺}\mathbf{P}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}{\mathbf{A}}^{⊺}\mathbf{b}+b{\mathbf{P}}_a\mathbf{a}=\hat{\mathbf{x}}-\frac{\mathbf{Pa}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}{\mathbf{a}}^{⊺}\hat{\mathbf{x}}+b{\mathbf{P}}_a\mathbf{a}$ where ${\mathbf{P}}_a\mathbf{a}=\mathbf{Pa}-\frac{\mathbf{Pa}{\mathbf{a}}^{⊺}\mathbf{Pa}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}=\frac{\mathbf{Pa}+\mathbf{Pa}{\mathbf{a}}^{⊺}\mathbf{Pa}-\mathbf{Pa}{\mathbf{a}}^{⊺}\mathbf{Pa}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}=\frac{\mathbf{Pa}}{1+{\mathbf{a}}^{⊺}\mathbf{Pa}}$ by expansion $\therefore {\hat{\mathbf{x}}}_{\text{new}}=\hat{\mathbf{x}}-{\mathbf{P}}_a\mathbf{a}{\mathbf{a}}^{⊺}\hat{\mathbf{x}}+b{\mathbf{P}}_a\mathbf{a}=\hat{\mathbf{x}}+{\mathbf{P}}_a\mathbf{a}(b-{\mathbf{a}}^{⊺}\hat{\mathbf{x}})$

So, ${\hat{\mathbf{x}}}_{\text{new}}$ can be obtained without additional inverse matrix calculation.

Facts

Sherman–Morrison formula is a special case of the Woodbury Formula

Link to original

Matrix Inversion Lemma

Definition

$(\mathbf{A}+\mathbf{UCV}{)}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{U}({\mathbf{C}}^{-1}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}{)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}$ where the size of matrices are $\mathbf{A}$ is $n\times n$, $\mathbf{C}$ is $k\times k$, and $\mathbf{V}$ is $k\times n$ and $\mathbf{A},\mathbf{C},({\mathbf{C}}^{-1}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U})$ should be invertible.

The inverse of a rank-k correction of some matrix can be computed by doing a rank-k correction to the inverse of the original matrix

Proof

Start by the matrix $\left[\begin{array}{cc}\mathbf{A}& \mathbf{U}\\ \mathbf{V}& \mathbf{C}\end{array}\right]$

The $LD\mathbf{U}$ decomposition of the original matrix become $\left[\begin{array}{cc}\mathbf{A}& \mathbf{U}\\ \mathbf{V}& \mathbf{C}\end{array}\right]=\left[\begin{array}{cc}I& 0\\ \mathbf{V}{\mathbf{A}}^{-1}& I\end{array}\right]\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& \mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\end{array}\right]\left[\begin{array}{cc}I& {\mathbf{A}}^{-1}\mathbf{U}\\ 0& I\end{array}\right]$ Inverting both sides gives^[Diagonalizable Matrix] $\begin{array}{rl}{\left[\begin{array}{cc}\mathbf{A}& \mathbf{U}\\ \mathbf{V}& \mathbf{C}\end{array}\right]}^{-1}& ={\left[\begin{array}{cc}I& {\mathbf{A}}^{-1}\mathbf{U}\\ 0& I\end{array}\right]}^{-1}{\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& \mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\end{array}\right]}^{-1}{\left[\begin{array}{cc}I& 0\\ \mathbf{V}{\mathbf{A}}^{-1}& I\end{array}\right]}^{-1}\\ & =\left[\begin{array}{cc}I& -{\mathbf{A}}^{-1}\mathbf{U}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{\mathbf{A}}^{-1}& 0\\ 0& {\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\end{array}\right]\left[\begin{array}{cc}I& 0\\ -\mathbf{V}{\mathbf{A}}^{-1}& I\end{array}\right]\\ & =\left[\begin{array}{cc}{\mathbf{A}}^{-1}+{\mathbf{A}}^{-1}\mathbf{U}{\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}& -{\mathbf{A}}^{-1}\mathbf{U}{\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\\ -{\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}& {\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\end{array}\right](1)\end{array}$

The $\mathbf{U}DL$ decomposition of the original matrix become $\left[\begin{array}{cc}\mathbf{A}& \mathbf{U}\\ \mathbf{V}& \mathbf{C}\end{array}\right]=\left[\begin{array}{cc}I& \mathbf{U}{\mathbf{C}}^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}& 0\\ 0& \mathbf{C}\end{array}\right]\left[\begin{array}{cc}I& 0\\ {\mathbf{C}}^{-1}\mathbf{V}& I\end{array}\right]$ Again inverting both sides, $\begin{array}{rl}{\left[\begin{array}{cc}\mathbf{A}& \mathbf{U}\\ \mathbf{V}& \mathbf{C}\end{array}\right]}^{-1}& ={\left[\begin{array}{cc}I& 0\\ {\mathbf{C}}^{-1}\mathbf{V}& I\end{array}\right]}^{-1}{\left[\begin{array}{cc}\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}& 0\\ 0& \mathbf{C}\end{array}\right]}^{-1}{\left[\begin{array}{cc}I& \mathbf{U}{\mathbf{C}}^{-1}\\ 0& I\end{array}\right]}^{-1}\\ & =\left[\begin{array}{cc}I& 0\\ -{\mathbf{C}}^{-1}\mathbf{V}& I\end{array}\right]\left[\begin{array}{cc}{\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}& 0\\ 0& {\mathbf{C}}^{-1}\end{array}\right]\left[\begin{array}{cc}I& -\mathbf{U}{\mathbf{C}}^{-1}\\ 0& I\end{array}\right]\\ & =\left[\begin{array}{cc}{\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}& -{\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}\mathbf{U}{\mathbf{C}}^{-1}\\ -{\mathbf{C}}^{-1}\mathbf{V}{\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}& {\mathbf{C}}^{-1}+{\mathbf{C}}^{-1}\mathbf{V}{\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}\mathbf{U}{\mathbf{C}}^{-1}\end{array}\right](2)\end{array}$

The first-row first-column entry of RHS of (1) and (2) above gives the Woodbury formula ${\left(\mathbf{A}-\mathbf{U}{\mathbf{C}}^{-1}\mathbf{V}\right)}^{-1}={\mathbf{A}}^{-1}+{\mathbf{A}}^{-1}\mathbf{U}{\left(\mathbf{C}-\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}\right)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}$

Link to original

Moore-Penrose Inverse

Definition

For a $p\times q$ matrix $A$, the Moore-Penrose inverse of $A$ is a $q\times p$ matrix $A^{+}$ is satisfying the following conditions

• $A{A}^{+}A=A$
• $A^{+}A{A}^{+}=A^{+}$
• $(A{A}^{+}{)}^{†}=A{A}^{+}$¹
• $(A^{+}A{)}^{†}=A^{+}A$¹

Facts

The pseudoinverse is defined and unique for all matrices whose entries are real and complex numbers.

If the matrix $A$ is a square matrix, then $A^{+}$ is also a square matrix and $A^{+}=A^{-1}$

Footnotes

1. Hermitian Matrix ↩ ↩²

Link to original

Partitioned Matrix

Block Matrix

Definition

A matrix that is interpreted as having been broken into sections called blocks or sub-matrices

Operations

Addition, Subtraction

$\left[\begin{array}{ccc}{\mathbf{A}}_{11}& \dots & {\mathbf{A}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{A}}_{m1}& \dots & {\mathbf{A}}_{mn}\end{array}\right]\pm \left[\begin{array}{ccc}{\mathbf{B}}_{11}& \dots & {\mathbf{B}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{B}}_{m1}& \dots & {\mathbf{B}}_{mn}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{A}}_{11}\pm {\mathbf{B}}_{11}& \dots & {\mathbf{A}}_{1n}\pm {\mathbf{B}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{A}}_{m1}\pm {\mathbf{B}}_{m1}& \dots & {\mathbf{A}}_{mn}\pm {\mathbf{B}}_{mn}\end{array}\right]$

Scalar Multiplication

$c\left[\begin{array}{ccc}{\mathbf{A}}_{11}& \dots & {\mathbf{A}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{A}}_{m1}& \dots & {\mathbf{A}}_{mn}\end{array}\right]=\left[\begin{array}{ccc}c{\mathbf{A}}_{11}& \dots & c{\mathbf{A}}_{1n}\\ ⋮& & ⋮\\ c{\mathbf{A}}_{m1}& \dots & c{\mathbf{A}}_{mn}\end{array}\right]$

Matrix Multiplication

$\left[\begin{array}{ccc}{\mathbf{B}}_{11}& \dots & {\mathbf{B}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{B}}_{k1}& \dots & {\mathbf{B}}_{kn}\end{array}\right]\left[\begin{array}{ccc}{\mathbf{A}}_{11}& \dots & {\mathbf{A}}_{1n}\\ ⋮& & ⋮\\ {\mathbf{A}}_{m1}& \dots & {\mathbf{A}}_{mn}\end{array}\right]=\left[\begin{array}{ccc}({\mathbf{B}}_{11}{\mathbf{A}}_{11}+\cdots +{\mathbf{B}}_{1m}{\mathbf{A}}_{m1})& \dots & ({\mathbf{B}}_{11}{\mathbf{A}}_{1n}+\cdots +{\mathbf{B}}_{1m}{\mathbf{A}}_{mn})\\ ⋮& & ⋮\\ ({\mathbf{B}}_{k1}{\mathbf{A}}_{11}+\cdots +{\mathbf{B}}_{km}{\mathbf{A}}_{m1})& \dots & ({\mathbf{B}}_{k1}{\mathbf{A}}_{1n}+\cdots +{\mathbf{B}}_{km}{\mathbf{A}}_{mn})\end{array}\right]$

Determinant

Let $D$ or $A$ be invertible matrices $\left|\begin{array}{cc}A& B\\ C& D\end{array}\right|=|D|\cdot |A-B{D}^{-1}C|=|A|\cdot |D-C{A}^{-1}B|$

Link to original

Eigenvalues and Eigenvectors

Eigendecomposition

Definition

$\mathbf{Av}=\lambda \mathbf{v}$ where $\mathbf{A}$ is a $n\times n$ matrix and $\mathbf{v}\ne 0$

If $\mathbf{Av}$ is a scalar multiple of non-zero vector $\mathbf{v}$, then $\lambda$ is the eigenvalue and $\mathbf{v}$ is the eigenvector.

Characteristic Polynomial

$|\mathbf{A}-\lambda \mathbf{I}|=0$ The values ${\lambda }_i$ satisfying the characteristic polynomial, are the eigenvalues of the matrix $\mathbf{A}$

Eigenspace

$E=\{\mathbf{v}:(\mathbf{A}-\lambda \mathbf{I})\mathbf{v}=0\}$

The set of all eigenvectors of $A$ corresponding to the same eigenvalue, together with the zero vector.

The Kernel of the matrix $(\mathbf{A}-\lambda \mathbf{I})$

Eigenvector: non-zero vector in the eigenspace

Algebraic Multiplicity

${\mu }_A({\lambda }_i)$

Let ${\lambda }_i$ be an eigenvalue of an $n\times n$ matrix $A$. The algebraic multiplicity ${\mu }_A({\lambda }_i)$ of the eigenvalue is its multiplicity as a root of a Characteristic Polynomial, that is, the largest k such that $(\lambda -{\lambda }_i{)}^k$

Geometric Multiplicity

${\gamma }_A({\lambda }_i)$ Let ${\lambda }_i$ be an eigenvalue of an $n\times n$ matrix $A$. The geometric multiplicity ${\gamma }_A({\lambda }_i)$ of the eigenvalue is the dimension of the Eigenspace associated with the eigenvalue.

Computation

1. Find the solution^[eigenvalues] of the Characteristic Polynomial.
2. Find the solution^[eigenvectors] of the Under-Constrained System $(\mathbf{A}-{\lambda }_i\mathbf{I}){\mathbf{v}}_i=0$ using the found eigenvalue.

Facts

There exists at least one eigenvector corresponding to the eigenvalue $\lambda$

Eigenvectors corresponding to different eigenvalues are always linearly independent.

When $\mathbf{A}$ is a normal or real Symmetric Matrix, the eigendecomposition is called Spectral Decomposition

$\mathrm{tr}(\mathbf{A})={\sum }_{i=1}^n{\lambda }_i$ The Trace of the matrix is equal to the sum of the eigenvalues of the matrix.

Proof Since the matrix $\mathbf{A}-\lambda \mathbf{I}$ only has $\lambda$ term on diagonal, and the calculation of cofactor deletes a row and a column, the coefficient of ${\lambda }^{n-1}$ of $\text{det}(\mathbf{A}-\lambda \mathbf{I})$ is always $-(a_{11}+a_{22}+\cdots +a_{nn})$. Also, since ${\lambda }_i$ are the solution of the Characteristic Polynomial, the expression is factorized as $\text{det}(\mathbf{A}-\lambda \mathbf{I})=(\lambda -{\lambda }_1)(\lambda -{\lambda }_2)\dots (\lambda -{\lambda }_n)$ and the coefficient of ${\lambda }^{n-1}$ become a $-({\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n)$ Therefore, $-(a_{11}+a_{22}+\cdots +a_{nn})=-({\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n)$ and $\text{Trace}(\mathbf{A})={\sum }_{i=1}^n{\lambda }_i$.

$|\mathbf{A}|={\prod }_{i=1}^n{\lambda }_i$ The Determinant of the matrix is equal to the product of the eigenvalues of the matrix.

Not all matrices have $n$ linearly independent eigenvectors.

When $\mathbf{Av}=\lambda \mathbf{v}$ holds,

• the eigenvalues of ${\mathbf{A}}^{-1}$ are $\frac{1}{\lambda }$ and the eigenvectors of the matrix ${\mathbf{A}}^{-1}$ are the same as $\mathbf{A}$.
• the eigenvalues of $A^k$ is ${\lambda }^k$
• the eigenvalues of $c\mathbf{A}$ is $c\lambda$
• the eigenvalues of $\mathbf{A}+c\mathbf{I}$ is $\lambda +c$
• the eigenvalues of $(\mathbf{A}+c\mathbf{I}{)}^{-1}$ is $1/(\lambda +c)$

$1\le {\gamma }_A({\lambda }_i)\le {\mu }_A({\lambda }_i)\le n$

An eigenvalue’s Geometric Multiplicity cannot exceed its Algebraic Multiplicity

$\exists P\text{ s.t. }A=PD{P}^{-1}\iff \forall {\lambda }_i:{\gamma }_A({\lambda }_i)={\mu }_A({\lambda }_i)$

the matrix is diagonalizable $\iff$ the Geometric Multiplicity is equal to the Algebraic Multiplicity for all eigenvalues

Let $A$ be a symmetric matrix. Then, $A$ has $r$ eigenvalues equal to $1$ and the rest zero $\iff A^2=A,\mathrm{rank}(A)=r$.

Link to original

The non-zero eigenvalues of $AB$ are the same as those of $BA$.

Link to original

Quadratic Forms and Positive Definite Matrix

Quadratic Form

Definition

$q(\mathbf{x})={\mathbf{x}}^{⊺}\mathbf{Ax}$ where $A$ is a Symmetric Matrix

A mapping $Q:V\to K$ where $V$ is a Module on Commutative Ring $K$ that has the following properties. $\forall k,l\in K,\forall u,v,w\in V$

• $Q(kv)=k^{2}Q(v)$
• $Q(u+v+w)=Q(u+v)+Q(v+w)+Q(u+w)-Q(u)-Q(v)-Q(w)$
• $Q(ku+lv)=k^{2}Q(u)+l^{2}Q(v)+klQ(u+v)-klQ(u)-klQ(v)$

Matrix Expressions

$a_1{x}_1^2+a_2{x}_2^2+2a_3{x}_1{x}_2\iff \left[\begin{array}{cc}{x}_1& x_2\end{array}\right]\left[\begin{array}{cc}{a}_1& a_3\\ a_3& a_2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$

$a_1{x}_1^2+a_2{x}_2^2+a_3{x}_3^2+2a_4{x}_1{x}_2+2a_5{x}_1{x}_3+2a_6{x}_2{x}_3\iff \left[\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right]\left[\begin{array}{ccc}{a}_1& a_4& a_5\\ a_4& a_2& a_6\\ a_5& a_6& a_3\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$

Facts

Let $\mathbf{x}$ be a Random Vector and $A$ be a symmetric matrix of constants. If $E(\mathbf{x})=\mathbf{\mu }$ and $\mathrm{Var}(\mathbf{x})=\Sigma$, then the expectation of the quadratic form $Q:={\mathbf{x}}^{⊺}A\mathbf{x}$ is $E(Q)={\mathbf{\mu }}^{⊺}A\mathbf{\mu }+\mathrm{tr}(A\Sigma )$

Let $\mathbf{x}=(X_1,\dots ,X_n{)}^{⊺}$ be a Random Vector and $A$ be a symmetric matrix of constants. If $E(\mathbf{x})=\mathbf{\theta }=({\theta }_1,\dots ,{\theta }_n{)}^{⊺}$ and $\mathrm{Var}(X_i)={\mu }_2$, $E[(X_i-{\theta }_i{)}^3]={\mu }_3$, and $E[(X_i-{\theta }_i{)}^4]={\mu }_4$, where $i=1,2,\dots ,n$, then the variance of the quadratic form ${\mathbf{x}}^{⊺}A\mathbf{x}$ is $\mathrm{Var}({\mathbf{x}}^{⊺}A\mathbf{x})=({\mu }_4-3{\mu }_2^2){\mathbf{a}}^{⊺}\mathbf{a}+2{\mu }_2^2\mathrm{tr}(A^2)+4{\mu }_2{\mathbf{\theta }}^{⊺}{A}^2\mathbf{\theta }+4{\mu }_3{\mathbf{\theta }}^{⊺}A\mathbf{a}$ where $\mathbf{a}$ is the column vector of diagonal elements of $A$.

If $X_i\sim N({\theta }_i,{\mu }_2),i=1,\dots ,n$ and $X_i$‘s are independent, then $\mathrm{Var}({\mathbf{x}}^{⊺}A\mathbf{x})=2{\mu }_2^2\mathrm{tr}(A^2)+4{\mu }_2{\mathbf{\theta }}^{⊺}{A}^2\mathbf{\theta }$ If $X_i\sim N(0,1),i=1,\dots ,n$ and $X_i$‘s are independent, then $\mathrm{Var}({\mathbf{x}}^{⊺}A\mathbf{x})=2\mathrm{tr}(A^2)$

Let ${\mathbf{X}}^{⊺}\sim N(0,{\sigma }^2{I}_n)$, $Q:={\mathbf{X}}^{⊺}A\mathbf{X}/{\sigma }^2$, where $A$ is a Symmetric Matrix and $\mathrm{rank}(A)=r\le n$, then the MGF of $Q$ is $M(t)={\prod }_{i=1}^r(1-2t{\lambda }_i{)}^{-1/2}=|I-2tA{|}^{-1/2}$ where ${\lambda }_i$‘s are non-zero eigenvalue of $A$

Let $\mathbf{X}\sim N_n(\mathbf{\mu },\Sigma )$, where $\Sigma$ is Positive-Definite Matrix, then $Q=(\mathbf{X}-\mathbf{\mu }{)}^{⊺}{\Sigma }^{-1}(\mathbf{X}-\mathbf{\mu })\sim {\chi }^2(n)$

Let $\mathbf{X}\sim N(0,{\sigma }^2{I}_n)$, $Q={\mathbf{X}}^{⊺}A\mathbf{X}/{\sigma }^2$, where $A$ is a Symmetric Matrix and $\mathrm{rank}(A)=r\le n$, then $Q\sim {\chi }^2(r)\iff A=A^k$ where $k\in \mathbb{N}$

Let $\mathbf{X}\sim N(0,{\sigma }^2{I}_n)$, $Q_1={\mathbf{X}}^{⊺}A\mathbf{X}/{\sigma }^2,Q_2={\mathbf{X}}^{⊺}B\mathbf{X}/{\sigma }^2$, where $A,B$ are symmetric matrices, then $Q_1,Q_2$ are independent if and only if $AB=O$

Let $Q=Q_1+Q_2+\cdots +Q_{k-1}+Q_k$, where $Q,Q_1,Q_2,\dots ,Q_k$ are quadratic forms in Random Sample from $N(0,{\sigma }^2)$ If $Q/{\sigma }^2\sim {\chi }^2(r),Q_1/{\sigma }^2\sim {\chi }^2(r_1),Q_2/{\sigma }^2\sim {\chi }^2(r_2),\dots ,Q_{k-1}/{\sigma }^2\sim {\chi }^2(r_{k-1})$ and $Q_k$ is non-negative, then

• $Q_1,Q_2,\dots ,Q_k$ are independent
• $Q_k/{\sigma }^2\sim {\chi }^2(r-r_1-r_2-\cdots -r_{k-1})$

Let $\mathbf{X}=(X_1,X_2,\dots ,X_n)\sim N(0,{\sigma }^2{I}_n)$, $\sum _{i=1}^n{X}_i^2=Q_1+Q_2+\cdots +Q_k$, where $Q_j={\mathbf{X}}^{⊺}{A}_j\mathbf{X}$, where $\mathrm{rank}(A_j)=r_j\le n$, then $\forall j=\{1,2,\dots ,k\},\frac{Q_j}{{\sigma }^2}\sim {\chi }^2(r_j)\iff {\sum }_{j=1}^k{r}_j=n$

${\mathbf{x}}^{⊺}A\mathbf{x}-2{\mathbf{x}}^{⊺}\mathbf{b}=(\mathbf{x}-A^{-1}\mathbf{b}{)}^{⊺}A(\mathbf{x}-A^{-1}\mathbf{b})-{\mathbf{b}}^{⊺}{A}^{-1}\mathbf{b}$

Link to original

Positive-Definite Matrix

Definition

$M\succ 0\iff \forall \mathbf{z}\ne 0,{\mathbf{z}}^{⊺}M\mathbf{z}>0$

Matrix $M$, in which ${\mathbf{z}}^{⊺}M\mathbf{z}$ is positive for every non-zero column vector $\mathbf{z}$ is a positive-definite matrix

Facts

Let $M\succ 0$, then

• $A^{-1}\succ 0$
• $\mathrm{rank}(CA{C}^{⊺})=\mathrm{rank}(C)$ where $C$ is a $p\times n$ matrix.
• $\mathrm{rank}(C)=p\Rightarrow CM{C}^{⊺}\succ 0$
• The diagonal elements of $M$ are positive.
• For a Symmetric Matrix $B$, $\exists t,M-tB\succ 0$ for sufficiently small $|t|$
• $\forall \mathbf{b},\underset{\mathbf{h}:\mathbf{h}\ne 0}{\sup }\frac{({\mathbf{h}}^{⊺}\mathbf{b}{)}^2}{{\mathbf{h}}^{⊺}M\mathbf{h}}={\mathbf{b}}^{⊺}{A}^{-1}\mathbf{b}$
• $\exists M^{1/2}\succ 0\text{s.t.}(M^{1/2}{)}^2=M$

$M\succ 0\iff \exists R\text{s.t.}A=R{R}^{⊺}$ where $R$ is a non-singular matrix.

$M\succ 0\iff$ all the leading minor determinants of $M$ are positive.

$\mathrm{rank}(X)=p\Rightarrow X^{⊺}X\succ 0$

Link to original

Projection and Decomposition of Matrix

Projection Matrix

Definition

$P=A(A^{⊺}A{)}^{-1}{A}^{⊺}$

For some vector $\mathbf{b}$, $P\mathbf{b}$ is the projection of $\mathbf{b}$ onto $A$

Facts

The projection matrix $P$ is symmetric Idempotent Matrix

Consider a Symmetric Matrix $P$, then $P$ is idempotent with rank $r$ if and only if $r$ eigenvalues are $1$ and $n-r$ eigenvalues are $0$.

$\mathrm{tr}(P)=\mathrm{rank}(P)$

The projection matrix is Positive Semi-Definite Matrix

If $P_1$ and $P_2$ are projection matrices, and $P_1-P_2$ is Positive Semi-Definite Matrix, then

• $P_1{P}_2=P_2{P}_1=P_2$
• $P_1-P_2$ is a projection matrix.

Link to original

QR Decomposition

Definition

$\mathbf{A}=\mathbf{QR}$ Decomposition of matrix $\mathbf{A}$ into a product $\mathbf{A}=\mathbf{QR}$ of an orthonormal matrix $\mathbf{Q}$ and an upper triangular matrix $\mathbf{R}$.

Computation

Using the Gram Schmidt Orthonormalization

QR decomposition can be computed by Gram Schmidt Orthonormalization. Where $Q=[q_1,\dots ,q_n]$ is the matrix of orthonormal column vectors obtained by the orthonormalization and $\mathbf{R}=\left[\begin{array}{ccccc}{\mathbf{q}}_1^{⊺}{\mathbf{a}}_1& {\mathbf{q}}_1^{⊺}{\mathbf{a}}_2& {\mathbf{q}}_1^{⊺}{\mathbf{a}}_3& \cdots & {\mathbf{q}}_1^{⊺}{\mathbf{a}}_n\\ 0& {\mathbf{q}}_2^{⊺}{\mathbf{a}}_2& {\mathbf{q}}_2^{⊺}{\mathbf{a}}_3& \cdots & {\mathbf{q}}_2^{⊺}{\mathbf{a}}_n\\ 0& 0& {\mathbf{q}}_3^{⊺}{\mathbf{a}}_3& \cdots & {\mathbf{q}}_3^{⊺}{\mathbf{a}}_n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & {\mathbf{q}}_n^{⊺}{\mathbf{a}}_n\end{array}\right]$

Link to original

Cholesky Decomposition

Definition

$A=L{L}^{†}$

Decomposition of a Positive-Definite Matrix into the product of lower triangular matrix and its Conjugate Transpose.

Computation

Let $\mathbf{A}=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{22}& a_{33}\end{array}\right]$ be Positive-Definite Matrix. Then, $\mathbf{A}=\mathbf{L}{\mathbf{L}}^{⊺}=\left[\begin{array}{ccc}{L}_{11}& 0& 0\\ L_{21}& L_{22}& 0\\ L_{31}& L_{32}& L_{33}\end{array}\right]\left[\begin{array}{ccc}{L}_{11}& L_{21}& L_{31}\\ 0& L_{22}& L_{32}\\ 0& 0& L_{33}\end{array}\right]=\left[\begin{array}{ccc}{L}_{11}^2& & \text{(symmetric)}\\ & & \\ L_{21}{L}_{11}& L_{21}^2+L_{22}^2& \\ & & \\ L_{31}{L}_{11}& L_{31}{L}_{21}+L_{32}{L}_{22}& L_{31}^2+L_{32}^2+L_{33}^2\end{array}\right]$ By setting the first-row first-column entry $L_{11}=\sqrt{a_{11}}$, can find other entries using substitution. $L_{21}=a_{21}/L_{11},L_{31}=a_{31}/L_{11},L_{22}=\sqrt{a_{22}-L_{21}^2},L_{32}=(a_{32}=L_{31}{L}_{21})/L_{22},L_{33}=\sqrt{a_{33}-L_{31}^2-L_{32}^2}$ By summarizing, $L_{ii}=\sqrt{a_{ii}-{\sum }_{k=1}^{i-1}{L}_{ik}^2}$, $L_{ij}=\frac{1}{L_{jj}}(a_{ij}-{\sum }_{k=1}^{j-1}{L}_{ik}{L}_{jk})$

Link to original

Spectral Theorem

Definition

$A=A^{†}\iff A=Q\mathbf{\Lambda }{Q}^{†}$ where $Q$ is a Unitary Matrix, and $\mathbf{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$

A matrix $A$ is a Hermitian Matrix if and only if $A$ is Unitary Diagonalizable Matrix.

Facts

Every Hermitian Matrix is diagonalizable, and has real-valued Eigenvalues and orthonormal eigenvector matrix.

For the Hermitian Matrix ${\mathbf{A}}^{†}=\mathbf{A}$, the every eigenvalue is real.

Proof For the eigendecomposition $\mathbf{Ax}=\lambda \mathbf{x}$, ${\mathbf{x}}^{†}\mathbf{Ax}={\mathbf{x}}^{†}\lambda \mathbf{x}=\lambda {\mathbf{x}}^{†}\mathbf{x}=\lambda ||\mathbf{x}||\Rightarrow \lambda =\frac{{\mathbf{x}}^{†}\mathbf{Ax}}{||\mathbf{x}||}$. Since $\mathbf{A}$ is hermitian, ${\mathbf{x}}^{†}\mathbf{Ax}\in \mathbb{R}$ and norm is always real. Therefore, every eigenvalue is real.

For the Hermitian Matrix ${\mathbf{A}}^{†}=\mathbf{A}$, the eigenvectors from different eigenvalues are orthogonal.

Proof Let the eigendecomposition $\mathbf{A}{\mathbf{x}}_1={\lambda }_1{\mathbf{x}}_1,\mathbf{A}{\mathbf{x}}_2={\lambda }_2{\mathbf{x}}_2$ where ${\lambda }_1\ne {\lambda }_2,{\mathbf{x}}_1\ne {\mathbf{x}}_2$. By the property of hermitian matrix, $(\mathbf{A}{\mathbf{x}}_1{)}^{†}{\mathbf{x}}_2={\mathbf{x}}_1^{†}{\mathbf{A}}^{†}{\mathbf{x}}_2={\mathbf{x}}_1^{†}{\lambda }_2{\mathbf{x}}_2={\lambda }_2{\mathbf{x}}_1^{†}{\mathbf{x}}_2$ and $(\mathbf{A}{\mathbf{x}}_1{)}^{†}{\mathbf{x}}_2=({\lambda }_1{\mathbf{x}}_1{)}^{†}{\mathbf{x}}_2={\bar{\lambda }}_1{\mathbf{x}}_1^{†}{\mathbf{x}}_2={\lambda }_1{\mathbf{x}}_1^{†}{\mathbf{x}}_2$ So, ${\lambda }_2{\mathbf{x}}_1^{†}{\mathbf{x}}_2={\lambda }_1{\mathbf{x}}_1^{†}{\mathbf{x}}_2\Rightarrow ({\lambda }_1-{\lambda }_2){\mathbf{x}}_1^{†}{\mathbf{x}}_2=0$. Since ${\lambda }_1\ne {\lambda }_2$, ${\mathbf{x}}_1^{†}{\mathbf{x}}_2=0$. Therefore, ${\mathbf{x}}_1^{†},{\mathbf{x}}_2$ are orthogonal.

Link to original

Singular Value Decomposition

Definition

${\mathbf{A}}_{n\times m}={\mathbf{U}}_{n\times n}{\mathbf{D}}_{n\times m}({\mathbf{V}}_{m\times m}{)}^{⊺}={\sum }_{i=1}^{\min (n,m)}{d}_i{U}_i{V}_i^{⊺}$ An arbitrary matrix ${\mathbf{A}}_{n\times m}$ can be decomposed to $\mathbf{UD}{\mathbf{V}}^{⊺}$.

For $n\ge m$ ${\mathbf{A}}_{n\times m}={\mathbf{U}}_{n\times m}{\mathbf{D}}_{m\times m}({\mathbf{V}}_{m\times m}{)}^{⊺}$ where ${\mathbf{U}}^{⊺}\mathbf{U}={\mathbf{I}}_m$ and ${\mathbf{V}}^{⊺}\mathbf{V}=\mathbf{V}{\mathbf{V}}^{⊺}={\mathbf{I}}_m$ and $\mathbf{D}$ is Diagonal Matrix

For $n\le m$ ${\mathbf{A}}_{n\times m}={\mathbf{U}}_{n\times n}{\mathbf{D}}_{n\times n}({\mathbf{V}}_{m\times n}{)}^{⊺}$ where ${\mathbf{U}}^{⊺}\mathbf{U}=\mathbf{U}{\mathbf{U}}^{⊺}={\mathbf{I}}_n$ and ${\mathbf{V}}^{⊺}\mathbf{V}={\mathbf{I}}_n$ and $\mathbf{D}$ is Diagonal Matrix

Calculation

For matrix ${\mathbf{A}}_{m\times n}=\mathbf{UD}{\mathbf{V}}^{⊺}$

${\mathbf{V}}_{n\times n}$ is the matrix of orthonormal eigenvectors of ${\mathbf{A}}^{⊺}\mathbf{A}$ ${\mathbf{U}}_{m\times m}$ is the matrix of orthonormal eigenvectors of $\mathbf{A}{\mathbf{A}}^{⊺}$ ${\mathbf{D}}_{m\times n}$ is the Diagonal Matrix made of the square roots of the non-zero eigenvalues of ${\mathbf{A}}^{⊺}\mathbf{A}$ and $\mathbf{A}{\mathbf{A}}^{⊺}$ sorted in descending order.

If the eigendecomposition is ${\mathbf{A}}^{⊺}\mathbf{Ax}=\lambda \mathbf{x}$, then the eigenvalues $\lambda \ge 0$ and the eigenvectors $\mathbf{x}\in \mathbb{C}({\mathbf{A}}^{⊺})$ are orthonormal by Spectral Theorem. If the $\mathrm{rank}(\mathbf{A})=r\le n$, then ${\lambda }_1⪈{\lambda }_2⪈\cdots \ge {\lambda }_r,{\lambda }_{r+1}=\dots {\lambda }_n=0$. Where ${\sigma }_i=\sqrt{{\lambda }_i}{}_{i=1,2,\dots ,r}$ are called the singular values.

Now, the Orthonormal Matrix $\mathbf{V}=[{\mathbf{V}}_1,{\mathbf{V}}_2]$ is calculated using the singular values. Where $V_1=[{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_r]$ is the Orthonormal Matrix of eigenvectors corresponding to the non-zero eigenvalues and $V_2=[{\mathbf{v}}_{r+1},{\mathbf{v}}_{r+2},\dots ,{\mathbf{v}}_n]$ is the Orthonormal Matrix of eigenvectors corresponding to zero eigenvalues where each ${\mathbf{v}}_r+i\in \mathbf{N}({\mathbf{A}}^{⊺}\mathbf{A})=\mathbf{N}(\mathbf{A})$, and ${\mathbf{D}}_{m\times n}=\mathrm{diag}({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r,0,\dots ,0)$ is a rectangular Diagonal Matrix.

Since $\mathbf{V}$ is an orthonormal matrix and $\mathbf{D}$ is a Diagonal Matrix, $\mathbf{AV}=\mathbf{UD}\Rightarrow \mathbf{A}[{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_n]=[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_m]\mathbf{D}\Rightarrow \mathbf{A}{\mathbf{v}}_i={\sigma }_i{\mathbf{u}}_i\Rightarrow {\mathbf{u}}_i=\frac{1}{{\sigma }_i}\mathbf{A}{\mathbf{v}}_i$. Now, the Orthonormal Matrix $\mathbf{U}=[{\mathbf{U}}_1,{\mathbf{U}}_2]$ is calculated using the linear system $\frac{1}{{\sigma }_i}\mathbf{A}{\mathbf{v}}_i$ and the null space of $\mathbf{A}$, $\mathbf{N}(\mathbf{A})$ Where ${\mathbf{U}}_1=[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_r]$ is the Orthonormal Matrix of the vectors ${\mathbf{u}}_{\{1,2,\dots ,r\}}\in \mathbf{C}(\mathbf{A})$ obtained from the system equation $\frac{1}{{\sigma }_i}\mathbf{A}{\mathbf{v}}_i$ And ${\mathbf{U}}_2=[{\mathbf{u}}_{r+1},{\mathbf{u}}_{r+2},\dots ,{\mathbf{u}}_m]$ is the Orthonormal Matrix of the vectors ${\mathbf{u}}_{\{r+1,r+2,\dots ,m\}}\in \mathbf{N}({\mathbf{A}}^{⊺})$. which is corresponding to zero eigenvalues

The Orthonormal Matrix $\mathbf{U}=[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_m]$ can also be formed by the eigenvectors of $\mathbf{A}{\mathbf{A}}^{⊺}$ similarly to calculating of $\mathbf{V}$.

Facts

$\mathbf{U},\mathbf{V}$ are the Unitary Matrix

Let $\mathbf{A}$ be a real symmetric Positive Semi-Definite Matrix Then, the Eigendecomposition(Spectral Decomposition) of $\mathbf{A}$ and the singular value decomposition of $\mathbf{A}$ are equal.

$\mathbf{A}=\mathbf{P\Lambda }{\mathbf{P}}^{⊺}=\mathbf{UD}{\mathbf{V}}^{⊺}$ where $\mathbf{P}=\mathbf{U}=\mathbf{V}$ and $\mathbf{\Lambda }=\mathbf{D}$ are non-negative and the shape of the every matrix is $n\times n$

Visualization

every matrix can be decomposed as a

• ${\mathbf{V}}^{⊺}$: rotation and reflection
• $\mathbf{D}$: scaling
• $\mathbf{U}$: rotation and reflection

Link to original

Miscellanea in Matrix

Centering Matrix

Definition

$C=I-\frac{1}{n}{J}_n$ where $J_n=1_n{1}^{⊺}$ is $n\times n$ matrix where all elements are $1$

Summing Vector

$1_n=(1,1,\dots ,1{)}^{⊺}$

Facts

$(n-1)s^2={\mathbf{x}}^{⊺}C\mathbf{x}$ where $s^2$ is a sample variance of $\mathbf{x}$

Link to original

Derivatives of Matrix

Definition

Differentiation by vector

$\frac{\partial }{\partial {\mathbf{x}}^{⊺}}:=\left[\frac{\partial }{\partial x_1},\frac{\partial }{\partial x_2},\dots ,\frac{\partial }{\partial x_n}\right]$ where $\mathbf{x}=[x_1,x_2,\dots ,x_n]$

Differentiation by Matrix

$\frac{\partial }{\partial X}:=\left[\begin{array}{ccc}\frac{\partial }{\partial x_{11}}& \dots & \frac{\partial }{\partial x_{1p}}\\ ⋮& \ddots & ⋮\\ \frac{\partial }{\partial x_{n1}}& \dots & \frac{\partial }{\partial x_{np}}\end{array}\right]$ where $X=\left[\begin{array}{ccc}{x}_{11}& \dots & x_{1p}\\ ⋮& \ddots & ⋮\\ x_{n1}& \dots & x_{np}\end{array}\right]$

Calculation

Let $f$ be a vector function, $\mathbf{x}$ be an $n$-dimensional vector, then $\frac{\partial f}{\partial {\mathbf{x}}^{⊺}}=\left[\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},\dots ,\frac{\partial f}{\partial x_n}\right]$ By the definition of derivative, we hold $\mathbf{df}=f(\mathbf{x}+d\mathbf{x})-f(\mathbf{x})=\frac{\partial f}{\partial x_1}d{x}_1+\frac{\partial f}{\partial x_2}d{x}_2+\cdots +\frac{\partial f}{\partial x_n}d{x}_n=\frac{\partial f}{\partial {\mathbf{x}}^{⊺}}d\mathbf{x}$ where $d\mathbf{x}=[d{x}_1,d{x}_2{]}^{⊺}$

Examples

Let $\mathbf{x},\mathbf{y}$ be $n$-dimensional vectors, then $\frac{\partial }{\partial \mathbf{x}}({\mathbf{x}}^{⊺}\mathbf{y})=\frac{\partial }{\partial \mathbf{x}}({\mathbf{y}}^{⊺}\mathbf{x})=\mathbf{y}$

Let $\mathbf{x}$ be an $n$-dimensional vector and $A$ be an $n\times n$ matrix, then $\frac{\partial }{\partial \mathbf{x}}({\mathbf{x}}^{⊺}A)=A,\frac{\partial }{\partial \mathbf{x}}(A\mathbf{x})=A^{⊺}$

Let $\mathbf{x}$ be an $n$-dimensional vector and $A$ be an $n\times n$ Symmetric Matrix, then the derivative of quadratic form is $\frac{\partial }{\partial \mathbf{x}}({\mathbf{x}}^{⊺}A\mathbf{x})=A\mathbf{x}+A^{⊺}\mathbf{x}=2A\mathbf{x}$

Link to original

Hessian Matrix

Definition

$H_f=\left\{\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\right\}=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \dots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2{x}_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \dots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \dots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]$

A square matrix of second-order partial derivative of a scalar-valued function

Link to original

Vectorization

Definition

Let $A=[{\mathbf{a}}_1,{\mathbf{a}}_2,\dots ,{\mathbf{a}}_n]$ be an $m\times n$ matrix, then $\mathrm{vec}(A)=[a_{11},a_{21},\dots ,a_{m1},a_{12},a_{22},\dots ,a_{m2},\dots ,a_{1n},a_{2n},\dots ,a_{mn}{]}^{⊺}$

A linear transformation which converts the matrix into a vector.

Facts

Let $A,B,C,Z$ be matrices, then

$\mathrm{vec}(ABC)=(C^{⊺}\otimes A)\mathrm{vec}(B)$ $\mathrm{tr}(AB)=\mathrm{vec}(A^{⊺}{)}^{⊺}\mathrm{vec}(B)$ $\mathrm{tr}(A{Z}^{⊺}BZC)=\mathrm{vec}(Z^{⊺}{)}^{⊺}(CA\otimes B^{⊺})\mathrm{vec}(Z)$

Link to original

Kronecker Product

Definition

Let $A$ be an $m\times n$ matrix and $B$ be a $p\times q$ matrix, then the Kronecker product $A\otimes B$ is the $pm\times qn$ block matrix $A\otimes B=\left[\begin{array}{ccc}{a}_{11}B& \cdots & a_{1n}B\\ ⋮& \ddots & ⋮\\ a_{m1}B& \cdots & a_{mn}B\end{array}\right]$

Facts

Let $A,B,X,Y$ be matrices, then

$(A\otimes B{)}^{⊺}=A^{⊺}\otimes B^{⊺}$ $(A\otimes B)(X\otimes Y)=AX\otimes BY$ $(A\otimes B{)}^{-1}=A^{-1}\otimes B^{-1}$, where $A,B$ are square matrices $\mathrm{rank}(A\otimes B)=\mathrm{rank}(A)\mathrm{rank}(B)$ $\mathrm{tr}(A\otimes B)=\mathrm{tr}(A)\mathrm{tr}(B)$ $|A_{p\times p}\otimes B_{m\times m}|=|A{|}^m|B{|}^p$ Eigenvalues of $A\otimes B$ is the product of the eigenvalue of $A$ and the eigenvalue of $B$

Link to original

Norms

Norm

Definition

$\Vert \cdot \Vert :V\to {\mathbb{R}}_0^{+}$ A real-valued function with the following properties

• Positive definiteness: $\forall \mathbf{v}\in V:\Vert \mathbf{v}\Vert =0\iff \mathbf{v}=0$
• Absolute Homogeneity: $\forall \mathbf{v}\in V,k\in K:\Vert k\mathbf{v}\Vert =|k|\Vert \mathbf{v}\Vert$
• Sub additivity or Triange inequality: $\forall \mathbf{v},\mathbf{w}\in V:\Vert \mathbf{v}+\mathbf{w}\Vert \le \Vert \mathbf{v}\Vert +\Vert \mathbf{w}\Vert$ where $K\in \{\mathbb{R},\mathbb{C}\}$ on a vector space $(V,+,\cdot ,K)$

Vector

Real

Let $\mathbf{x}\in {\mathbb{R}}^n$ be an $n$-dimensional vector, then the $L_p$ norm of the $\mathbf{x}$ is defined as $\Vert \mathbf{x}{\Vert }_p:={\left({\sum }_{i=1}^n|x_i{|}^p\right)}^{1/p}$

Complex

Let $\mathbf{x}\in {\mathbb{C}}^n$ be an $n$-dimensional complex vector, then the $L_1$ norm of the $\mathbf{x}$ is defined as $\Vert \mathbf{x}\Vert ={\mathbf{x}}^{†}\mathbf{x}=\sqrt{{\sum }_{i=1}^n{\bar{x}}_i{x}_i}$

Function

$\Vert f\Vert =\sqrt{⟨f,f⟩}=\sqrt{{\int }_a^{b}f(x{)}^{2}dx}$: norm of $f$

Facts

$\Vert \frac{f(x)}{\Vert f\Vert }\Vert =1$

a norm can be induced by a Inner product $d(v,w)=\Vert v-w\Vert$

Link to original

Schatten Norm

Definition

Let $A$ be an $m\times n$ matrix, then the Schatten norm of $A$ is defined as $||A|{|}_p={\left({\sum }_{i=1}^{\min \{m,n\}}{\sigma }_i^p\right)}^{1/p}$ where ${\sigma }_i$ is the $i$-th singularvalue of $A$

Facts

$p=1$: Nuclear Norm

$p=2$: Frobenius Norm $p=\infty$: Spectral Norm

Link to original

Frobenius Norm

Definition

Let $A$ be an $m\times n$ matrix, then the Frobenius norm of $A$ is defined as $||A|{|}_F={\left({\sum }_{i=1}^m{\sum }_{h=1}^n|a_{ij}{|}^2\right)}^{1/2}=\sqrt{\mathrm{tr}(A^{†}A)}$

Link to original

Matrix p-Norm

Definition

Let $A$ be an $m\times n$ matrix, then the matrix p-norm of $A$ is defined as $||A|{|}_p=\underset{\mathbf{x}\ne 0}{\sup }\frac{||A\mathbf{x}|{|}_p}{||x|{|}_p}=\underset{||\mathbf{x}|{|}_p=1}{\max }||A\mathbf{x}|{|}_p$

Facts

When $p=2$, the norm is a Spectral Norm.

Link to original

Spectral Norm

Definition

Let $A$ be an $m\times n$ matrix, then the matrix p-norm of $A$ is defined as $||A|{|}_2=\sqrt{{\lambda }_{\max }(A^{†}A)}={\sigma }_{\max }(A)$

Link to original

Nuclear Norm

Definition

Let $A$ be an $m\times n$ matrix, then the nuclear norm of $A$ is defined as $||A|{|}_{\ast }=\mathrm{tr}(\sqrt{A^{†}A)})={\sum }_{i=1}^{\min \{m,n\}}{\sigma }_i(A)$

Link to original

L-pq Norm

Definition

$L_{pq}$ norm is an entry-wise matrix norm. $||A|{|}_{p,q}={\left(\sum _{i=1}^n{\left(\sum _{j=1}^n|a_{ij}{|}^p\right)}^{q/p}\right)}^{1/q}$ where $p,q\ge 1$

Facts

When $p=q=1$, the norm is the sum of the absolute values of every entry and is called a matrix $L_1$ norm.

Link to original

Rayleigh-Ritz Theorem

Definition

Let $A$ be an $n\times n$ Hermitian Matrix with eigenvalues sorted in descending order ${\lambda }_1,\dots ,{\lambda }_p$, where $p\le n$, then $\underset{\mathbf{x}\ne 0}{\max }\frac{{\mathbf{x}}^{⊺}A\mathbf{x}}{{\mathbf{x}}^{⊺}\mathbf{x}}=\underset{||\mathbf{x}|{|}_2=1}{\max }{\mathbf{x}}^{⊺}A\mathbf{x}={\lambda }_1$ $\underset{\mathbf{x}\ne 0}{\min }\frac{{\mathbf{x}}^{⊺}A\mathbf{x}}{{\mathbf{x}}^{⊺}\mathbf{x}}=\underset{||\mathbf{x}|{|}_2=1}{\min }{\mathbf{x}}^{⊺}A\mathbf{x}={\lambda }_p$

Link to original

Link to original

Elementary Statistical Theory

Random Variable and Random Sampling

Random Variable

Definition

$X:\Omega \to \mathbb{R}\text{s.t.}\forall B\in \mathcal{R},\exists \omega \in \Omega ,\text{s.t.}{X}^{-1}(B)=\{\omega :X(\omega )\in B\}\in \mathcal{F}$

A random variable is a function $X:\Omega \to \mathbb{R}$ whose inverse function $X^{-1}(B)$ is $\mathcal{F}$-measurable for the two measurable spaces $(\Omega ,\mathcal{F})$ and $(\mathbb{R},\mathcal{R})$.

The inverse image of an arbitrary Borel set of Codomain $\mathbb{R}$ is an element of sigma field $\mathcal{F}$.

Notations

Consider a probability space $(\Omega ,\mathcal{F},\mathbb{P})$

• Outcomes: $H,T$
• Set of outcomes (Sample space): $\Omega =\{H,T\}$
• Events: $\varnothing ,\{H\},\{T\},\Omega$
• Set of events (Sigma-Field): $\mathcal{F}=\{\varnothing ,\{H\},\{T\},\{H,T\}\}$
• Probabilities: $\mathbb{P}:\mathcal{F}\to [0,1]$
• Random variable: $X:\Omega \to \mathbb{R}$

For a random variable $X$ on a Probability Space $(\Omega ,\mathcal{F},P)$

• $X\sim {\mu }_X$ if and only if the Distribution of $X$ is ${\mu }_X$
• $X\sim F_X$ if and only if the Distribution Function of $X$ is $F_X$

For a random variable $X$ on Probability Space $({\Omega }_X,{\mathcal{F}}_X,P_X)$ and another random variable $Y$ on Probability Space $({\Omega }_Y,{\mathcal{F}}_Y,P_Y)$

• $X\stackrel{d}{=}Y⟺\forall B\in \mathcal{R}:{\mu }_X(B)={\mu }_Y(B)$
• $X\stackrel{d}{=}Y⟺\forall k\in \mathbb{R}:F_X(k)=F_Y(k)$
• $X\stackrel{d}{=}Y⟺\forall k\in \mathbb{R}:P_X(X\le k)=P_Y(Y\le k)$

Link to original

Expected Value

Definition

$E(X)={\int }_{\Omega }XdP={\int }_{\mathbb{R}}xd{\mu }_X={\int }_{\mathbb{R}}xd{F}_X$ The expected value of the Random Variable $X$ on Probability Space $(\Omega ,\mathcal{F},P)$

Continuous

$E(X)={\int }_{-\infty }^{\infty }x{f}_X(x)dx$ where $f_X(x)$ is a PDF of Random Variable $X$

The expected value of the Random Variable $X$ when $F_X$ satisfies absolute continuous over $\lambda$, ${\mu }_X<<\lambda$

Discrete

$E(X)={\sum }_{x}x{p}_X(x)$ where $p_X(x)$ is a PMF of Random Variable $X$

The expected value of the Random Variable $X$ when $F_X$ satisfies absolute continuous over $\#$, ${\mu }_X<<{\#}_X$ In other words, $F_X$ has at most countable jumps.

Expected Value of a Function

$E(g(x))={\int }_{\mathbb{R}}g(x)d{\mu }_X$

Continuous

$E(g(x))={\int }_{\mathbb{R}}g(x)f_X(x)dx$

Discrete

$E(g(x))=\sum _{x\in S_X}g(x)p_X(x)$

Properties

Linearity

Random Variables

If $\exists k_{i}E[X_i]$, then $E\left[{\sum }_{i=1}^n{k}_i{X}_i\right]={\sum }_{i=1}^n{k}_{i}E(X_i)$

Matrix of Random Variables

Let $W_1,W_2$ be a $m\times n$ matrices of random variables, $A_1,A_2$ be $k\times m$ matrices of constants, and $B$ a $n\times l$ matrix of constant. Then, $E[A_1{W}_1+A_2{W}_2]=A_{1}E[W_1]+A_{2}E[W_2]$

$E[A_1{W}_{1}B]=A_{1}E[W_1]B$

Notations

Expression	Discrete	Continuous
Expression for the event $2^{\Omega }$ and the probability $P$	${\sum }_{\omega \in \Omega }X(\omega )\cdot P(\omega )$	${\int }_{\Omega }X(\omega )\cdot dP(\omega )$
Expression for the measurable space $(\mathbb{R},\mathcal{R})$ and the distribution ${\mu }_X$	${\sum }_{x\in \mathbb{R}}x\cdot f(x)$	${\int }_{\mathbb{R}}f(x)\lambda (x)={\int }_{\mathbb{R}}\frac{d{\mu }_X}{d\lambda }\lambda (x)={\int }_{\mathbb{R}}d{\mu }_X={\int }_{\mathbb{R}}d{F}_X$

Link to original

Multivariate Distribution

Definition

Joint CDF

The Distribution Function $F_{\mathbf{X}}:{\mathbb{R}}^n\to [0,1]$ of a Random Vector $\mathbf{X}=(X_1,\dots ,X_n{)}^{⊺}$ is defined as $F_{\mathbf{X}}(\mathbf{x})=P(\mathbf{X}\le \mathbf{x})=P({\mathbf{X}}_1\le x_1,\dots ,{\mathbf{X}}_n\le x_n)$ where $\mathbf{x}=(x_1,\dots ,x_n{)}^{⊺}$

Joint PDF

$p_{X_1,\dots ,X_n}(x_1,\dots ,x_n)=P(X_1=x_1,\dots ,X_n=x_n)$

The joint probability Distribution Function of $n$ discrete Random Vector $X_1,X_2,\dots ,X_n$

$f_{X_1,\dots ,X_n}(x_1,\dots ,x_n)=\frac{{\partial }^n{F}_{X_1,\dots ,X_n}(x_1,\dots ,x_n)}{\partial x_1\dots \partial x_n}$

The joint probability Distribution Function of $n$ continuous Random Vector $X_1,X_2,\dots ,X_n$

Expected Value of a Multivariate Function

Continuous

E[g(\mathbf{X})] = \idotsint\limits_{x_{n} \dots x_{1}} g(\mathbf{x})f(\mathbf{x})dx_{1} \dots \dots dx_{n}

Discrete

$E[g(\mathbf{X})]=\sum _{x_n}\cdots \sum _{x_1}g(\mathbf{x})f(\mathbf{x})$

Marginal Distribution of a Multivariate Function

Marginal CDF

$F_{X_1}(x_1)=\underset{x_2,\dots ,x_n\to \infty }{\lim }{F}_{\mathbf{X}}(\mathbf{x})$

Marginal PDF

f_{X_{1}}(x_{1}) = E[g(\mathbf{X})] = \idotsint\limits_{x_{n} \dots x_{1}} f(x_{2}, \dots, x_n)dx_{2} \dots \dots dx_{n}

Conditional Distribution of a Multivariate Function

$f_{2,\dots ,n|1}(x_2,\dots ,x_n|x_1)=\frac{f_{\mathbf{X}}(\mathbf{x})}{f_{X_1}(x_1)}$

$f_{1|2,\dots ,n}(x_1|x_2,\dots ,x_n)=\frac{f_{\mathbf{X}}(\mathbf{x})}{f_{X_2,\dots ,X_n}(x_2,\dots ,x_n)}$

Properties

Linearity

If $\exists E[g_1(X_1,X_2)],E[g_2(X_1,X_2)]$, then $E[k_1{g}_1(X_1,X_2)+k_2{g}_2(X_1,X_2)]=k_{1}E[g_1(X_1,X_2)]+k_{2}E[g_2(X_1,X_2)]$

Link to original

Covariance

Definition

$\mathrm{Cov}(X,Y)=E[(X-{\mu }_X)(Y-{\mu }_Y)]$

Properties

Covariance with Itself

$\mathrm{Cov}(X,X)=\mathrm{Var}(X)$

Covariance of Linear Combinations

$\mathrm{Cov}(\sum _{i=1}^n{a}_i{X}_i,\sum _{i=1}^n{b}_i{Y}_i)=\sum _{i=1}^n\sum _{j=1}^m{a}_i{b}_j\mathrm{Cov}(X_i,Y_j)$

$\mathrm{Var}(\sum _{i=1}^n{a}_i{X}_i)=\mathrm{Cov}(\sum _{i=1}^n{a}_i{X}_i,\sum _{i=1}^n{a}_i{X}_i)=\sum _{i=1}^n{a}_i^2\mathrm{Var}(X_i)+2\sum _{i,j:i<j}{a}_i{a}_j\mathrm{Cov}(X_i,X_j)$

Link to original

Transclude of Density-Function

Distribution Function

Definition

$F_X(x):\mathbb{R}\to [0,1]:={\mu }_X((-\infty ,x])=P(X\le x)$

A distribution function $F_X$ is a function for the Random Variable on Probability Space $(\Omega ,\mathcal{F},P)$

Facts

$(\Omega ,\mathcal{F},P)\sim (\mathbb{R},\mathcal{R},{\mu }_X)$

Proof By definition, $F_X(x)={\mu }_X((-\infty ,x])$ So, defining $\forall x\in \mathbb{R},F_X(x)$ is the Equivalence Relation to defining $\mathcal{A}=\{(-\infty ,x]:x\in \mathbb{R}\},{\mu }_X:\mathcal{A}\to [0,1]$

Since $\mathcal{A}$ is a Pi-System, ${\mu }_X:\sigma (\mathcal{R})=\mathcal{R}\to [0,1]$ is uniquely determined by ${\mu }_X:\mathcal{A}\to [0,1]$ by Extension from Pi-System

Now, $(\mathbb{R},\mathcal{R},{\mu }_X)$ is a Measurable Space induced by $X$. Therefore, defining ${\mu }_X$ on $(\mathbb{R},\mathcal{R})$ is equivalent to the defining $P$ on $(\Omega ,\mathcal{F})$

Distribution function(CDF) $F:\mathbb{R}\to [0,1]$ has the following properties

• Monotonic increasing: $\forall a,b,\in \mathbb{R},a<b\Rightarrow F(a)\le F(b)$
• $\underset{x\to -\infty }{\lim }F(x)=0$
• $\underset{x\to \infty }{\lim }F(x)=1$
• Right-continuous: $\underset{x\to x_0^{+}}{\lim }F(x)=F(x_0)$

If a function $F:\mathbb{R}\to [0,1]$ satisfies the following properties, then $F$ is a distribution function(CDF) of some Random Variable $X$

• Monotonic increasing: $\forall a,b,\in \mathbb{R},a<b\Rightarrow F(a)\le F(b)$
• $\underset{x\to -\infty }{\lim }F(x)=0$
• $\underset{x\to \infty }{\lim }F(x)=1$
• Right-continuous: $\underset{x\to x_0^{+}}{\lim }F(x)=F(x_0)$

$\forall a,b\in \mathbb{R},a<b\Rightarrow (P(a<X\le b)=F_X(b)-F_X(a))$

Link to original

Bernoulli Distribution

Definition

$X\sim \mathrm{Ber}(p)=B(1,p)$ where $p\in [0,1]$ is a probability of success

number of success in a single trial with success probability $p$

Bernoulli Process

The i.i.d. Random Vector $X_1,\dots ,X_n$ with Bernoulli distribution

Properties

PDF

$f(x)=p^x(1-p{)}^{1-x},x\in \{0,1\}$

Mean

$E(X)=p$

Variance

$\mathrm{Var}(X)=p(1-p)$

Link to original

Binomial Distribution

Definition

$X\sim B(n,p)$ where $n$ is the length of bernoulli process, and $p\in [0,1]$ is a probability of success

The number of successes in length $n$ bernoulli process with success probability $p$

Properties

PDF

$f(x)=\left(\genfrac{}{}{0ex}{}{n}{x}\right)p^x(1-p{)}^{n-x},x\in \{0,1,\dots ,n\}$

Mean

$E(X)=np$

Variance

$\mathrm{Var}(X)=np(1-p)$

MGF

$M_X(t)=[(1-p)+p{e}^t{]}^n$

Facts

Let $X_i\sim B(n_i,p)$ be independent random variables following binomial distribution, then $Y=\sum _{i=1}^n{X}_i\sim B(\sum _{i=1}^n{n}_i,p)$

Link to original

Multinomial Distribution

Definition

$X\sim \mathrm{Mult}(n,p_1,p_2,\dots ,p_{k-1})$ where $n$ is the number of trials, $p_i$ is a probability of category $i$

Distribution that describes the probability of observing a specific combination of outcomes

Properties

PDF

$f(x_1,x_2,\dots ,x_{k-1})=n!{\prod }_{i=1}^k\frac{p_i^{x_i}}{x_i!}$ where $x_k=n-\sum _{i=1}^{k-1}{x}_i$, $p_k=1-\sum _{i=1}^{k-1}{p}_i$, and $0\le \sum _{i=1}^{k-1}{x}_i\le n$

MGF

$M(t_1,t_2,\dots ,t_{k-1})=({\sum }_{i=1}^{k-1}{p}_i{e}^{t_i}+p_k{)}^n$

Marginal PDF

Each one-variable marginal pdf is Binomial Distribution, each two-variables marginal pdf is Trinomial Distribution, and so on.

Link to original

Poission Distribution

Definition

$X\sim \mathrm{Pois}(\lambda )$ where $\lambda$ is the average number of occurrences in a fixed interval of time

The number of occurrences in a fixed interval of time with mean $\lambda$

Properties

PDF

$f(x)=\frac{e^{-\lambda }{\lambda }^x}{x!}\in [0,1]$ where $x\in {\mathbb{N}}_0$

Mean

$E(X)=\lambda$

Variance

$\mathrm{Var}(X)=\lambda$

MGF

$M(t)=\exp [m(e^t-1)]$

Summation

Let $X_i\sim \mathrm{Pois}({\lambda }_i)$, and $X_i$‘s are independent, then $\sum _{i=1}^n{X}_i\sim \mathrm{Pois}\left(\sum _{i=1}^n{\lambda }_i\right)$

Link to original

Normal Distribution

Definition

$X\sim N(\mu ,{\sigma }^2)$ where $\mu \in \mathbb{R}$ is the location parameter(mean), and ${\sigma }^2\in {\mathbb{R}}_0^{+}$ is the scale parameter(variance)

Standard Normal Distribution

$X\sim N(0,1)$ $f(x)=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{x^2}{2}\right)\in \mathbb{R}$

Properties

PDF

$f(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{1}{2}(\frac{x-\mu }{\sigma }{)}^2\right)\in \mathbb{R}$

Mean

$E(X)=\mu$

Variance

$\mathrm{Var}(X)={\sigma }^2$

MGF

$M_X(t)=\exp \left(\mu t+\frac{{\sigma }^2{t}^2}{2}\right)$

Higher Order Moments

$E(Z^k)=\{\begin{array}{l}E(Z^{2k})=\frac{(2k)!}{2^{k}k!}\\ E(Z^{2k+1})=0\end{array}$

$E(X^k)=\sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right){\sigma }^{j}E(Z^j){\mu }^{k-j}$

Sum of Normally Distributed Random Variables

Let $X_i\sim N({\mu }_i,{\sigma }_i^2)$ be independent random variables following normal distribution, then $\sum _{i=1}^n{a}_i{X}_i\sim N\left(\sum _{i=1}^n{a}_i{\mu }_i,{\sum }_{i=1}^n{a}_i^2{\sigma }_i^2\right)$

Relationship with Chi-squared Distribution

Let $Z\sim N(0,1)$ be a standard normal distribution, then $Z^2\sim {\chi }^2(1)$

Facts

Let $X\sim N(0,{\sigma }^2)$, then $\mathrm{Var}(X^2)=2{\sigma }^4$

Link to original

Chi-squared Distribution

Definition

$X\sim {\chi }^2(k)=\Gamma \left(\frac{k}{2},2\right)$ where $k\in \mathbb{N}$ is the degrees of freedom

squared sum of independent standard normal distributions

Properties

PDF

$f(x)=\frac{1}{\Gamma (\frac{k}{2})2^{\frac{k}{2}}}{x}^{\frac{k}{2}-1}{e}^{-\frac{x}{2}}\in {\mathbb{R}}_0^{+}$

Mean

$E(X)=k$

Variance

$\mathrm{Var}(X)=2k$

MGF

$M(t)=(1-2t{)}^{-k/2}$

Additivity

Let $X_i\sim {\chi }^2(k_i)$‘s are independent chi-squared distributions $\sum _{i=1}^n{X}_i\sim {\chi }^2\left(\sum _{i=1}^n{k}_i\right)$

Facts

Let $Q_1\sim {\chi }^2(r_1)$ and $Q_2\sim {\chi }^2(r_2)$, and $Q=Q_1-Q_2$ is independent of $Q_2$, then $Q\sim {\chi }^2(r_1-r_2)$

Let $Q=Q_1+Q_2+\cdots +Q_{k-1}+Q_k$, where $Q,Q_1,Q_2,\dots ,Q_k$ are quadratic forms in $\mathbf{x}$, where each element of the $\mathbf{x}$ is a Random Sample from $N(\mu ,{\sigma }^2)$ If $Q/{\sigma }^2\sim {\chi }^2(r),Q_1/{\sigma }^2\sim {\chi }^2(r_1),Q_2/{\sigma }^2\sim {\chi }^2(r_2),\dots ,Q_{k-1}/{\sigma }^2\sim {\chi }^2(r_{k-1})$, then

• $Q_1,Q_2,\dots ,Q_k$ are independent
• $Q_k/{\sigma }^2\sim {\chi }^2(r-r_1-r_2-\cdots -r_{k-1})$

Link to original

Student's t-Distribution

Definition

Let $W\sim N(0,1)$ be a standard normal distribution, $V\sim {\chi }^2(r)$ be a Chi-squared Distribution, and $W,V$ be independent, then $\frac{W}{\sqrt{V/r}}\sim t(r)$ where $r\in {\mathbb{R}}^{+}$ is the degrees of freedom

Properties

PDF

$f_r(t)=\frac{\Gamma (\frac{r+1}{2})}{\sqrt{\pi r}\Gamma (\frac{r}{2})}{\left(1+\frac{t^2}{r}\right)}^{-(r+1)/2}\in \mathbb{R}$

Mean

$E(X)=0$

Variance

$\mathrm{Var}(X)=\{\begin{array}{ll}\frac{r}{r-2}& r>2\\ \infty & 1<r\le 1\end{array}$

Link to original

F-Distribution

Definition

Let $U\sim {\chi }^2(r_1),V\sim {\chi }^2(r_2)$ be independent random variables following Chi-squared distributions, then $\frac{U/r_1}{V/r_2}\sim F(r_1,r_2)$ where $r_1,r_2\in \mathbb{N}$ are the degrees of freedoms

Properties

PDF

{\Gamma(\frac{r_{1} + r_{2}}{2}) (\frac{r_{1}}{r_{2}})^{\frac{r_{1}}{2}} x^{\frac{r_{1}}{2}-1}} {\Gamma(\frac{r_{1}}{2}) \Gamma(\frac{r_{2}}{2}) (\frac{r_{1}}{r_{2}}x+1)^{(r_{1} + r_{2})/2}}$$ ## Mean $$E(X) = \frac{r_{2}}{r_{2} - 2},\quad r_{2}>2$$ ## Variance $$\operatorname{Var}(X) = \frac{2(r_{1}+r_{2}-2)}{r_{1} (r_{2}-2)^{2}(r_{2}-4)},\quad r_{2}>4$$Link to original

Random Vector

Definition

$\mathbf{X}:(\Omega ,\mathcal{F})\to ({\mathbb{R}}^d,{\mathcal{R}}^d),\mathbf{X}=(X_1,X_2,\dots ,X_d)$

Column vector whose components are random variables $X:(\Omega ,\mathcal{F})\to (\mathbb{R},\mathcal{R})$.

$\mathbf{X}(\omega )=(X_1,X_2,\dots ,X_d)(\omega )=(X_1(\omega ),X_2(\omega ),\dots ,X_d(\omega ))=(x_1,x_2,\dots ,x_d)$

Link to original

Multivariate Normal Distribution

Definition

$\mathbf{X}\sim N_n(\mathbf{\mu },\Sigma )$ where $n$ is the number of dimensions, $\mathbf{\mu }$ is the vector of location parameters, and $\mathbf{\Sigma }$ is the vector of scale parameters

Standard Multivariate Normal Distribution

$\mathbf{X}\sim N_n(0,I_n)$

PDF

$f(\mathbf{z})=(2\pi {)}^{-\frac{n}{2}}\exp \left(-\frac{1}{2}{\mathbf{z}}^{⊺}\mathbf{z}\right)$

MGF

$M_{\mathbf{z}}(\mathbf{t})=\exp \left(\frac{1}{2}{\mathbf{t}}^{⊺}\mathbf{t}\right)$

Properties

PDF

$f(\mathbf{x})=(2\pi {)}^{-\frac{n}{2}}|\Sigma {|}^{-\frac{1}{2}}\exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{⊺}{\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu })\right)$

Mean

$E(\mathbf{X})=\mu$

Variance

$\mathrm{Var}(\mathbf{X})=\Sigma$

MGF

$M_{\mathbf{x}}(\mathbf{t})=\exp \left({\mathbf{t}}^{⊺}\mathbf{\mu }+\frac{1}{2}{\mathbf{t}}^{⊺}\Sigma \mathbf{t}\right)$

Affine Transformation

Let $\mathbf{X}\sim N_n(\mathbf{\mu },\Sigma )$ be a Random Variable following multivariate normal distribution, $A$ be a $m\times n$ matrix, and $\mathbf{b}$ be a $m$ dimensional vector, then $A\mathbf{X}+\mathbf{b}\sim N_m(A\mathbf{\mu }+\mathbf{b},A\Sigma A^{⊺})$

Relationship with Chi-squared Distribution

Suppose $\mathbf{X}\sim N_n(\mathbf{\mu },\Sigma )$ be a Random Variable following multivariate normal distribution, then ${\Sigma }^{-\frac{1}{2}}(\mathbf{X}-\mathbf{\mu })\sim N_n(0,I_n)$ $(\mathbf{X}-\mathbf{\mu }{)}^{⊺}{\Sigma }^{-1}(\mathbf{X}-\mathbf{\mu })\sim {\chi }^2(n)$

Facts

Let $\mathbf{X}\sim N_n(\mathbf{\mu },\Sigma )$, $C$ be an $m\times n$, and $\mathbf{d}$ be a $m$-dimensional vector, then $C\mathbf{X}+\mathbf{D}\sim N_m(C\mu +\mathbf{d},C\Sigma C^{⊺})$

Let $\mathbf{Y}\sim N_n(\mathbf{\mu },\Sigma )$, $\mathbf{Y}=(Y_1,Y_2{)}^{⊺}$, $\mathbf{\mu }=({\mu }_1,{\mu }_2{)}^{⊺}$, $\mathbf{\Sigma }=\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right)$, where $Y_1$ is $p\times 1$ and $Y_2$ is $(n-p)\times 1$ vectors. Then, $Y_1$ and $Y_2$ are independent $\iff {\Sigma }_{12}=O_{p\times (n-p)}$

Let $\mathbf{Y}\sim N_n(\mathbf{\mu },\Sigma )$, $\mathbf{u}=A\mathbf{Y}$, and $\mathbf{v}=B\mathbf{Y}$, where $A$ is $m\times n$, $B$ is $l\times n$ matrices. $\mathbf{u}$ and $\mathbf{v}$ are independent $\iff \mathrm{Cov}(\mathbf{u},\mathbf{v})=A\Sigma B^{⊺}=O_{m\times l}$

Link to original

Statistical Estimation and Testing

Bias of an Estimator

Definition

Let $X_1,X_2,\dots ,X_n$ be a Random Sample from $f(\mathbf{x}|\mathbf{\theta })$ where $\mathbf{\theta }\in \Omega$ is a parameter, and $\hat{\mathbf{\theta }}=T(X_1,\dots ,X_n)$ be a Statistic. $\mathrm{Bias}(\hat{\mathbf{\theta }})=E((\hat{\mathbf{\theta }}))-\mathbf{\theta }$

$\forall \mathbf{\theta }\in \Omega ,E(\hat{\mathbf{\theta }})=\mathbf{\theta }\Rightarrow \hat{\mathbf{\theta }}$ is unbiased estimator An estimator is unbiased if its bias is equal to zero for all values of parameter $\mathbf{\theta }$.

Link to original

Consistency

Definition

${\hat{\mathbf{\theta }}}_n\stackrel{p}{\to }\mathbf{\theta }$

A Statistic ${\hat{\mathbf{\theta }}}_n$ is called consistent estimator of $\mathbf{\theta }$ if $\hat{\mathbf{\theta }}$ converges in probability to $\mathbf{\theta }$

Link to original

Minimum Variance Unbiased Estimator

Definition

An estimator $Y=u(X_1,X_2,\dots ,X_n)$ satisfying the following is the minimum variance unbiased estimator (MVUE) for $\theta$ $(E(Y)=\theta )\wedge (\forall T,(E(T)=\theta )\Rightarrow \mathrm{Var}(Y)\le \mathrm{Var}(T))$ where $T$ is an unbiased estimator

An Unbiased Estimator that has lower variance than any other unbiased estimator for the parameter $\theta$

Facts

A minimum variance unbiased estimator does not always exist.

If some unbiased estimator’s variance is equal to the Rao-Cramer Lower Bound for all $\theta$, then it is a minimum variance unbiased estimator.

Link to original

Least Squares Estimator

Definition

Simple Linear Regression Case

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The least square estimator is the estimator that minimizes $\sum _{i=1}^n{ϵ}_i^2=\sum _{i=1}^n(Y_i-{\beta }_0+{\beta }_1{X}_i{)}^2$

The least square estimator of the model is ${\hat{\beta }}_0=\bar{Y}-{\hat{\beta }}_1\bar{X},{\hat{\beta }}_1=\frac{S_{XY}}{S_{XX}}$ where $S_{XY}=\sum _{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})$ and $S_{XX}=\sum _{i=1}^n(X_i-\bar{X}{)}^2$

Estimation of ${\sigma }^2$

$s^2=\frac{1}{n-2}{\sum }_{i=1}^n{e}_i^2$ where $e_i=(Y_i-{\hat{Y}}_i)$‘s are residuals

Multiple Linear Regression Case

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ The least square estimator is the estimator that minimizes ${ϵ}^{⊺}ϵ=(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^{⊺}(\mathbf{y}-\mathbf{X}\mathbf{\beta })$

The least square estimator of the model is $\hat{\mathbf{\beta }}=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$

Fitted response vector is expressed as $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}=\mathbf{X}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$ where $H:=\mathbf{X}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}$ is called the hat matrix.

Estimation of ${\sigma }^2$

$s^2=\frac{{\mathbf{e}}^{⊺}\mathbf{e}}{n-p}$ where $\mathbf{e}=(e_1,e_2,\dots ,e_n)$, $e_i=(Y_i-{\hat{Y}}_i)$‘s are residuals, and $p$ is the number of the explanatory variables

Facts

The least square estimator ${\hat{\beta }}_1$ is a linear combination of $Y_i$‘s Let $w_i=\frac{(X-\bar{X})}{S_{XX}}$, then ${\hat{\beta }}_1=\sum _{i=1}^n{w}_i{Y}_i$

$E({\hat{\beta }}_1)={\beta }_1$ $E({\hat{\beta }}_0)={\beta }_0$ $\mathrm{Var}({\hat{\beta }}_1)=\frac{{\sigma }^2}{S_{XX}}$ $\mathrm{Var}({\hat{\beta }}_0)={\sigma }^2(\frac{1}{n}+\frac{{\bar{X}}^2}{S_{XX}})$

$\sum _{i=1}^n{e}_i=0$ $\sum _{i=1}^n{X}_i{e}_i=0,{\mathbf{X}}^{⊺}\mathbf{ϵ}=0$ $\sum _{i=1}^n{\hat{Y}}_i{e}_i=0,{\hat{\mathbf{y}}}^{⊺}ϵ=0$ The fitted line always go through $(\bar{X},\bar{Y})$

$E(\hat{\mathbf{\beta }})=\mathbf{\beta }$ $\mathrm{Cov}(\hat{\beta })=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\sigma }^2$ $E(s^2)={\sigma }^2$

Let $\mathbf{y}=(Y_1,\dots ,Y_n{)}^{⊺}$, where $Y_i$‘s are independent, $E(\mathbf{Y})=\mathbf{X}\mathbf{\beta }$, and ${\sigma }^2,{\mu }_3,{\mu }_4$ be the 2nd, 3rd, and 4th Central Moment of $Y_i$ respectively. Then, $s^2$ is the unique non-negative quadratic Unbiased Estimator of ${\sigma }^2$ with minimum variance when the excess kurtosis is $({\mu }_4-3{\sigma }^4)/{\sigma }^4=0$ or when the diagonal elements of the hat matrix $X(X^{⊺}X{)}^{-1}{X}^{⊺}$ are equal.

Link to original

Maximum Likelihood Estimation

Definition

MLE is the method of estimating the parameters of an assumed Distribution

Let $X_1,X_2,\dots ,X_n$ be Random Sample with PDF $f(x|\theta )$, where $\theta \in \Omega$, then the MLE ${\hat{\theta }}_{\text{MLE}}$ of $\theta$ is estimated as ${\hat{\theta }}_{\text{MLE}}=\underset{\theta }{\mathrm{argmax}}L(\theta |\mathbf{x})$

Regularity Conditions

• R0: The pdfs are distinct, i.e. $\theta \ne {\theta }^{\prime }\Rightarrow f(x_i|\theta )\ne f(x_i|{\theta }^{\prime })$
• R1: The pdfs have same supports $\forall \theta$
• R2: The true value ${\theta }_0$ is an interior point in $\Omega$
• R3: The pdf $f(x|\theta )$ is twice differentiable with respect to $\theta \in \Omega$
• R4: $\frac{\partial }{\partial {\theta }^2}\int f(x|\theta )dx=\int \frac{\partial }{\partial {\theta }^2}f(x|\theta )dx$
• R5: The pdf $f(x|\theta )$ is three times differentiable with respect to $\theta \in \Omega$, $\forall \theta \in \Omega ,\left|\frac{{\partial }^3}{\partial {\theta }^3}\ln f(x|\theta )\right|\le M(x)$, and $\exists c\in \mathbb{R},\exists M(x),\forall |\theta -{\theta }_0|<c,\forall$ interior point $x,E_{{\theta }_0}[M(X)]<\infty$

Properties

Functional Invariance

If $\hat{\theta }$ is the MLE for $\theta$, then $g(\hat{\theta })$ is the MLE of $g(\theta )$

Consistency

Under R0 ~ R2 Regularity Conditions, let ${\theta }_0$ be a true parameter, $f(x|\theta )$ is differentiable with respect to $\theta \in \Omega$, then $\frac{\partial }{\partial \theta }L(\theta )=0$ has a solution ${\hat{\theta }}_n$ such that ${\hat{\theta }}_n\stackrel{P}{\to }{\theta }_0$

Asymptotic Normality

Under the R0 ~ R5 Regularity Conditions, let $X_1,X_2,\dots ,X_n$ be Random Sample with PDF $f(x|\theta )$, where $\theta \in \Omega$, ${\hat{\theta }}_n$ be a consistent Sequence of solutions of MLE equation $\frac{\partial l(\theta )}{\partial \theta }=0$, and $0<I({\theta }_0)<\infty$, then $\sqrt{n}({\hat{\theta }}_n-{\theta }_0)\stackrel{D}{\to }N\left(0,\frac{1}{I({\theta }_0)}\right)$ where $I({\theta }_0)$ is the Fisher Information.

By the asymptotic normality, the MLE estimator is asymptotically efficient under R0 ~ R5 Regularity Conditions

Asymptotic Confidence Interval

By the asymptotic normality of MLE, $\sqrt{nI(\hat{\theta })}(\hat{\theta }-\theta )\stackrel{D}{\to }N(0,1)$ Thus, $100(1-\alpha )\%$ confidence interval of for $\theta$ is $\left(\hat{\theta }-z_{\alpha /2}\frac{1}{\sqrt{nI(\hat{\theta })}},\hat{\theta }+z_{\alpha /2}\frac{1}{\sqrt{nI(\hat{\theta })}}\right)$

Delta method for MLE Estimator

Under the R0 ~ R5 Regularity Conditions, let $g(x)$ be a continuous function and $g^{\prime }({\theta }_0)\ne 0$, then $\sqrt{n}(g({\hat{\theta }}_n)-g({\theta }_0))\stackrel{D}{\to }N\left(0,\frac{g^{\prime }({\theta }_0{)}^2}{I({\theta }_0)}\right)$

Facts

Under R0 and R1 regularity conditions, let ${\theta }_0$ be a true parameter, then $\forall \theta \ne {\theta }_0,\underset{n\to \infty }{\lim }{P}_{{\theta }_0}[L({\theta }_0)>L(\theta )]=1$

Link to original

Taylor Approximation

Taylor Series

Definition

Taylor Series for One Variable

Let $f:\mathbb{R}\to \mathbb{R}$ be an infinitely differentiable function at the point $c\in \mathbb{R}$, then $f(x)=\sum _{n=0}^{\infty }\frac{f^{(n)}(c)}{n!}(x-c{)}^n=f(a)+f^{\prime }(a)(x-a)+\frac{f^{\prime \prime }(a)}{2!}(x-a{)}^2+\cdots +\frac{f^{(m)}(a)}{m!}(x-a{)}^m+\dots$

Taylor Series for Two Variables

Let $f:{\mathbb{R}}^2\to \mathbb{R}$ be an infinitely differentiable function at the point $(t_0,x_0)\in {\mathbb{R}}^2$, then

&+ \frac{\partial f}{\partial t}(t_{0}, x_{0})(t-t_{0}) + \frac{\partial f}{\partial x}(t_{0}, x_{0})(x-x_{0})\\ &+ \frac{1}{2}\frac{\partial^{2} f}{\partial t^{2}}(t_{0}, x_{0})(t-t_{0})^{2} + \frac{1}{2}\frac{\partial^{2} f}{\partial x^{2}}(t_{0}, x_{0})(x-x_{0})^{2} + \frac{\partial^{2} f}{\partial t \partial x}(t_{0}, x_{0})(x-x_{0})(t-t_{0})\\ &+ \dots \end{aligned}$$ ## Approximation Let $f: \mathbb{R} \to \mathbb{R}$ be a $k$-the differentiable function at the point $c \in \mathbb{R}$ and $|x-c| \to 0$, then $$f(x) = \sum\limits_{n=0}^{k} \cfrac{f^{(n)}(c)}{n!}(x-c)^{n} + o(|x-c|^{k})$$ where $o$ is [[Big O Notation|little oh]] notation Let $f: \mathbb{R} \to \mathbb{R}$ be a $k+1$-the differentiable function at any point and $|f^{(k+1)}(x)| \leq M \leq \infty$, then by the [[Rolle's Theorem]] $$f(x) = \sum\limits_{n=0}^{k} \cfrac{f^{(n)}(c)}{n!}(x-c)^{n} + \frac{1}{(k+1)!}f^{(k+1)}(\xi)(x-c)^{k+1}$$ where $\xi \in (x, c)$ ## Maclaurin Series $$f(x) = \sum\limits_{n=0}^{\infty} \cfrac{f^{(n)}(0)}{n!}x^{n} = f(0)+f'(0)x + \cfrac{f''(0)}{2!}x^{2}+ \dots +\cfrac{f^{(m)}(0)}{m!}x^m + \dots$$ Taylor Series with $c=0$ ## Derivation from the [[Fundamental Theorem of Calculus]] Let $f$ be an infinitely [[Differentiability|differentiable]] function. By applying [[Fundamental Theorem of Calculus|FTC]] $n$ times, we can expande the function $f(x)$ as follows: $$\begin{aligned} f(x) &= \int\limits_{c}^{x} f'(x_{1}) dx_{1} + f(c)\\ &= \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} f''(x_{2}) dx_{2} + f'(c) dx_{1} + f(c)\\ &= \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} f''(x_{2}) dx_{2} dx_{1} + \int\limits_{c}^{x} f'(c) dx_{1} + f(c)\\ &= \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\int\limits_{c}^{x_{2}} f'''(x_{3}) dx_{3} dx_{2} dx_{1} + \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} f''(c) dx_{2} dx_{1} + \int\limits_{c}^{x} f'(c) dx_{1} + f(c)\\ &= \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{n-1}} f^{(n)}(x_{n}) dx_{n} dx_{n-1} \dots dx_{1} + \sum_{i=0}^{n} \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{i-1}} f^{(i)}(c) dx_{i} dx_{i-1} \dots dx_{1} \end{aligned}$$ The integrals inside the summation, can be simplified $$\begin{aligned} \int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{i-1}} f^{(i)}(c) dx_{i} dx_{i-1} \dots dx_{1} &= f^{(i)}(c)\int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{i-1}} 1 dx_{i} dx_{i-1} \dots dx_{1} \\ &= f^{(i)}(c)\int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{i-2}} (x_{i-1} - c) dx_{i-1} dx_{i-2} \dots dx_{1} \\ &= f^{(i)}(c)\int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{i-3}} \frac{1}{2}((x_{i-2} - c)^2 - (c - c)^2) dx_{i-2} dx_{i-3} \dots dx_{1} \\ &\dots \\ &= f^{(i)}(c)\frac{(x - c)^i}{i!} \end{aligned}$$ Therefore, we have $$f(x) = \underbrace{\int\limits_{c}^{x} \!\int\limits_{c}^{x_{1}} \!\dots \!\int\limits_{c}^{x_{n-1}} f^{(n)}(x_{n}) dx_{n} dx_{n-1} \dots dx_{1}}_{\text{Error term}} + \underbrace{\sum_{i=0}^{n} f^{(i)}(c)\frac{(x - c)^i}{i!}}_{\text{Polynomial terms}}$$Link to original

Delta Method

Definition

Univariate Delta Method

Let $\{X_n\}$ be a sequence of random variables satisfying $\sqrt{n}(X_n-\theta )\stackrel{D}{\to }N(0,{\sigma }^2)$, $g$ be a differentiable function at $x=\theta$, and $g^{\prime }(\theta )\ne 0$, then $\sqrt{n}(g(X_n)-g(\theta ))\stackrel{D}{\to }{g}^{\prime }(\theta )N(0,{\sigma }^2)=N(0,{\sigma }^2{g}^{\prime }(\theta {)}^2)$

Proof

By Taylor Series approximation $g(X_n)=g(\theta )+g^{\prime }(\theta )(X_n-\theta )+o_P(|X_n-\theta |)$ $\Rightarrow \sqrt{n}(g(X_n)-g(\theta ))=g^{\prime }(\theta )\sqrt{n}(X_n-\theta )+o_P(\sqrt{n}|X_n-\theta |)$

Where $\sqrt{n}(X_n-\theta )\stackrel{D}{\to }N(0,{\sigma }^2)$ by the assumption. By the continuous mapping theorem, $\sqrt{n}|X_n-\theta |$, which is a function of $X_n$, also converges to random variable, so it is Boundedness in Probability $\sqrt{n}|X_n-\theta |=O_P(1)$ by the property of converging random vector

Therefore, $\sqrt{n}(g(X_n)-g(\theta ))=g^{\prime }(\theta )N(0,{\sigma }^2)+o_P(O_P(1))=g^{\prime }(\theta )N(0,{\sigma }^2)=N(0,{\sigma }^2{g}^{\prime }(\theta {)}^2)$ by the property of sequence of random variables bounded in probability

Examples

Estimation of the sample variance of Bernoulli Distribution

Let $X_1,X_2,\dots ,X_n\sim Ber(p)$ $\sqrt{n}({\bar{X}}_n-p)\to N(0,p(1-p))$ by CLT

$\sqrt{n}(g({\bar{X}}_n)-g(p))\to g^{\prime }(p)N(0,p(1-p))$ by Delta method

Let $g(x)=x(1-x)$, then $\sqrt{n}({\bar{X}}_n(1-{\bar{X}}_n)-p(1-p))\to N(0,p(1-p)(1-2p{)}^2)$

Therefore, the sample mean and variance follow such distributions ${\bar{X}}_n\to N\left(p,\frac{p(1-p)}{n}\right)=N\left(\mu ,\frac{\sigma }{n}\right)$ ${\bar{X}}_n(1-{\bar{X}}_n)=\hat{\sigma }\to N\left(p(1-p),\frac{p(1-p)(1-2p{)}^2}{n}\right)=N\left(\sigma ,\frac{\sigma (1-2p{)}^2}{n}\right)$

Visualization

• x-axis: $\mu =p$
• y_axis: ${\sigma }^2=p(1-p)$
• y1: $y=p(1-p)$^[variance by mean]
• y2: ${\bar{X}}_n$^[sample mean]
• y3, y4: ${\bar{X}}_n(1-{\bar{X}}_n)$^[sample variance that calculated by the sample mean]
• y5: first order approximated line at $(p,p(1-p))$

If sample size $n\to \infty$, variance of sample mean $V({\bar{X}}_n)\to 0$. So, sample variance $\hat{\sigma }$ can be well approximated by first-order approximation.

Link to original

Newton's Method

Definition

An iterative algorithm for finding the roots of a differentiable function, which are solution to the equation $f(x)=0$

Algorithm

Find the next point such that the Taylor series of the given point is 0 Taylor first approximation: $f(x)\simeq f(x_n)+f^{\prime }(x_n)(x-x_n)$ The point such that the Taylor series is 0: $x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$

multivariate version: ${\mathbf{x}}_{n+1}=f({\mathbf{x}}_n)-\nabla f({\mathbf{x}}_n{)}^{-1}f({\mathbf{x}}_n)$

In convex optimization,

Find the minimum point^[its derivative is 0] of Taylor quadratic approximation. Taylor quadratic approximation: $f(x)\simeq f(x_n)+f^{\prime }(x_n)(x-x_n)+\frac{f^{\prime \prime }(x_n)}{2}(x-x_n{)}^2$ The derivative of the quadratic approximation: $f^{\prime }(x_n)+f^{\prime \prime }(x_n)(x-x_n)$ The minimum point of the quadratic approximation^[the point such that the derivative of the quadratic approximation is 0]: $x_{n+1}=x_n-\frac{f^{\prime }(x_n)}{f^{\prime \prime }(x_n)}$ multivariate version: ${\mathbf{x}}_{n+1}={\mathbf{x}}_n-{\nabla }^{2}f({\mathbf{x}}_n{)}^{-1}\nabla f({\mathbf{x}}_n)$

Examples

Solution of a linear system

Solve $A\mathbf{x}=\mathbf{b}$ with an MSE loss

The cost function is $f(\mathbf{x})=(A\mathbf{x}-\mathbf{b}{)}^{⊺}(A\mathbf{x}-\mathbf{b})$ and its gradient and hessian are $\nabla f(x)=-2(\mathbf{b}-A\mathbf{x}{)}^{⊺}$, ${\nabla }^{2}f(x)=2A^{⊺}A$

Then, solution is ${\mathbf{x}}_{n+1}={\mathbf{x}}_n+(A^{⊺}A{)}^{-1}{A}^{⊺}(\mathbf{b}-A{\mathbf{x}}_n)$ If $(A^{⊺}A)$ is invertible, ${\mathbf{x}}_{n+1}=(A^{⊺}A{)}^{-1}{A}^{⊺}\mathbf{b}$ is a Least Square solution.

Link to original

Multiple Testing

Family-Wise Error Rate

Definition

Consider a result of Multiple Testing

fact \ test result	do not reject $H_0$	reject $H_0$	total
$H_0$ is true	$U$	$V$	$m_0$
$H_1$ is true	$T$	$S$	$m_1$
total	$m-R$	$R$	$m$
where $R=V+S$ is the number of rejected null hypothesis, and $m$ is the total number of hypotheses tested.

$FWER=P(V\ge 1|m_0>0)$

The family-wise error rate (FWER) is the probability of making at least one type 1 error in the family Therefore, by assuring $FWER=P(V\ge 1|m_0>0)\le \alpha$, the probability of making one or more type 1 errors in the family is controlled at level $\alpha$

Link to original

Bonferroni's Method

Definition

Consider $m$ Multiple Testing problem $H_{0i}:{\mu }_{1i}={\mu }_{2i},i=1,2,\dots ,m$, and let $E_i$ be an event of type 1 error for $H_{0i}$. Then, by the Bonferroni Inequality, $FWER=P(V\ge 1)=P(\bigcup _{i=1}^m{E}_i)\le \sum _{i=1}^{m}P(E_i)$ Therefore, to control Family-Wise Error Rate $P(V\ge 1)$ at $\alpha$, we need to control $P(E_i)\le \frac{\alpha }{m}$

The Bonferroni method is too conservative and consequently, it is not a powerful test.

Link to original

Sidak Method

Definition

Consider $m$ Multiple Testing problem $H_{0i}:{\mu }_{1i}={\mu }_{2i},i=1,2,\dots ,m$, and let $E_i$ be an event of type 1 error for $H_{0i}$ By the definition of the type 1 error in $m$ multiple testing, $FWER=P(V\ge 1)=P(\bigcup _{i=1}^m{E}_i)\le 1-(1-\alpha {)}^m$ Therefore, to control Family-Wise Error Rate $P(V\ge 1)$ at $\alpha$, Sidak method uses ${\alpha }^{\ast }=1-(1-\alpha {)}^{1/m}$ instead of the $\alpha$

Link to original

False Discovery Rate

Definition

$FDR=E(V/R)$ where $V$ is the number of false discoveries and $S$ is the number of discoveries

The false discovery rate (FDR) is the expected proportion of falsely rejected $H_0$ we control FDR at some $q$, $FDR\le q$

Facts

If Family-Wise Error Rate is controlled, then FDR also be controlled.

Let $p_1\le p_2\le \cdots \le p_m$ be the $p$-values of $m$ hypotheses, and $q$ is given. If $j=\underset{1\le 1\le m}{\max }\left\{i|p_i\le \frac{i}{m}q\right\}$, then FDR rejects hypotheses corresponding to $p_1\le p_2\le \cdots \le p_j$

Link to original

Simple Linear Regression Model

Simple Linear Regression Model

Simple Linear Regression

Definition

$Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i,i=1,2,\dots ,n$ where ${ϵ}_i$‘s are i.i.d. error terms, with $E({ϵ}_i)=0$ and $\mathrm{Var}({ϵ}_i)={\sigma }^2$

Link to original

Estimation of Regression Coefficients

Least Squares Estimator

Definition

Simple Linear Regression Case

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The least square estimator is the estimator that minimizes $\sum _{i=1}^n{ϵ}_i^2=\sum _{i=1}^n(Y_i-{\beta }_0+{\beta }_1{X}_i{)}^2$

The least square estimator of the model is ${\hat{\beta }}_0=\bar{Y}-{\hat{\beta }}_1\bar{X},{\hat{\beta }}_1=\frac{S_{XY}}{S_{XX}}$ where $S_{XY}=\sum _{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})$ and $S_{XX}=\sum _{i=1}^n(X_i-\bar{X}{)}^2$

Estimation of ${\sigma }^2$

$s^2=\frac{1}{n-2}{\sum }_{i=1}^n{e}_i^2$ where $e_i=(Y_i-{\hat{Y}}_i)$‘s are residuals

Multiple Linear Regression Case

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ The least square estimator is the estimator that minimizes ${ϵ}^{⊺}ϵ=(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^{⊺}(\mathbf{y}-\mathbf{X}\mathbf{\beta })$

The least square estimator of the model is $\hat{\mathbf{\beta }}=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$

Fitted response vector is expressed as $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}=\mathbf{X}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$ where $H:=\mathbf{X}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}$ is called the hat matrix.

Estimation of ${\sigma }^2$

$s^2=\frac{{\mathbf{e}}^{⊺}\mathbf{e}}{n-p}$ where $\mathbf{e}=(e_1,e_2,\dots ,e_n)$, $e_i=(Y_i-{\hat{Y}}_i)$‘s are residuals, and $p$ is the number of the explanatory variables

Facts

The least square estimator ${\hat{\beta }}_1$ is a linear combination of $Y_i$‘s Let $w_i=\frac{(X-\bar{X})}{S_{XX}}$, then ${\hat{\beta }}_1=\sum _{i=1}^n{w}_i{Y}_i$

$E({\hat{\beta }}_1)={\beta }_1$ $E({\hat{\beta }}_0)={\beta }_0$ $\mathrm{Var}({\hat{\beta }}_1)=\frac{{\sigma }^2}{S_{XX}}$ $\mathrm{Var}({\hat{\beta }}_0)={\sigma }^2(\frac{1}{n}+\frac{{\bar{X}}^2}{S_{XX}})$

$\sum _{i=1}^n{e}_i=0$ $\sum _{i=1}^n{X}_i{e}_i=0,{\mathbf{X}}^{⊺}\mathbf{ϵ}=0$ $\sum _{i=1}^n{\hat{Y}}_i{e}_i=0,{\hat{\mathbf{y}}}^{⊺}ϵ=0$ The fitted line always go through $(\bar{X},\bar{Y})$

$E(\hat{\mathbf{\beta }})=\mathbf{\beta }$ $\mathrm{Cov}(\hat{\beta })=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\sigma }^2$ $E(s^2)={\sigma }^2$

Let $\mathbf{y}=(Y_1,\dots ,Y_n{)}^{⊺}$, where $Y_i$‘s are independent, $E(\mathbf{Y})=\mathbf{X}\mathbf{\beta }$, and ${\sigma }^2,{\mu }_3,{\mu }_4$ be the 2nd, 3rd, and 4th Central Moment of $Y_i$ respectively. Then, $s^2$ is the unique non-negative quadratic Unbiased Estimator of ${\sigma }^2$ with minimum variance when the excess kurtosis is $({\mu }_4-3{\sigma }^4)/{\sigma }^4=0$ or when the diagonal elements of the hat matrix $X(X^{⊺}X{)}^{-1}{X}^{⊺}$ are equal.

Link to original

Least Absolute Deviation Estimator

Definition

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The least absolute deviation estimator ${\hat{\beta }}_0,{\hat{\beta }}_1$ are the estimator that minimize $\sum _{i=1}^n|{ϵ}_i|=\sum _{i=1}^n|Y_i-{\beta }_0+{\beta }_1{X}_i|$

Link to original

Gauss-Markov Theorem

Definition

Simple Linear Regression Case

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ Let $E({ϵ}_i)=0$ and $\forall i\ne j,\mathrm{Cov}(e_i,e_j)=0$, then the Least Squares Estimator ${\hat{\beta }}_0,{\hat{\beta }}_1$ has the minimum variance among the all unbiased linear estimators, i.e. BLUE

Multiple Linear Regression Case

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ Let $E(\mathbf{ϵ})=0,\mathrm{Cov}(\mathbf{ϵ})=I_n{\sigma }^2$, then the Least Squares Estimator $\hat{\mathbf{\beta }}$ has the minimum variance among all the unbiased and linear estimators, i.e. BLUE

Link to original

Maximum Likelihood Estimation

Definition

MLE is the method of estimating the parameters of an assumed Distribution

Let $X_1,X_2,\dots ,X_n$ be Random Sample with PDF $f(x|\theta )$, where $\theta \in \Omega$, then the MLE ${\hat{\theta }}_{\text{MLE}}$ of $\theta$ is estimated as ${\hat{\theta }}_{\text{MLE}}=\underset{\theta }{\mathrm{argmax}}L(\theta |\mathbf{x})$

Regularity Conditions

• R0: The pdfs are distinct, i.e. $\theta \ne {\theta }^{\prime }\Rightarrow f(x_i|\theta )\ne f(x_i|{\theta }^{\prime })$
• R1: The pdfs have same supports $\forall \theta$
• R2: The true value ${\theta }_0$ is an interior point in $\Omega$
• R3: The pdf $f(x|\theta )$ is twice differentiable with respect to $\theta \in \Omega$
• R4: $\frac{\partial }{\partial {\theta }^2}\int f(x|\theta )dx=\int \frac{\partial }{\partial {\theta }^2}f(x|\theta )dx$
• R5: The pdf $f(x|\theta )$ is three times differentiable with respect to $\theta \in \Omega$, $\forall \theta \in \Omega ,\left|\frac{{\partial }^3}{\partial {\theta }^3}\ln f(x|\theta )\right|\le M(x)$, and $\exists c\in \mathbb{R},\exists M(x),\forall |\theta -{\theta }_0|<c,\forall$ interior point $x,E_{{\theta }_0}[M(X)]<\infty$

Properties

Functional Invariance

If $\hat{\theta }$ is the MLE for $\theta$, then $g(\hat{\theta })$ is the MLE of $g(\theta )$

Consistency

Under R0 ~ R2 Regularity Conditions, let ${\theta }_0$ be a true parameter, $f(x|\theta )$ is differentiable with respect to $\theta \in \Omega$, then $\frac{\partial }{\partial \theta }L(\theta )=0$ has a solution ${\hat{\theta }}_n$ such that ${\hat{\theta }}_n\stackrel{P}{\to }{\theta }_0$

Asymptotic Normality

Under the R0 ~ R5 Regularity Conditions, let $X_1,X_2,\dots ,X_n$ be Random Sample with PDF $f(x|\theta )$, where $\theta \in \Omega$, ${\hat{\theta }}_n$ be a consistent Sequence of solutions of MLE equation $\frac{\partial l(\theta )}{\partial \theta }=0$, and $0<I({\theta }_0)<\infty$, then $\sqrt{n}({\hat{\theta }}_n-{\theta }_0)\stackrel{D}{\to }N\left(0,\frac{1}{I({\theta }_0)}\right)$ where $I({\theta }_0)$ is the Fisher Information.

By the asymptotic normality, the MLE estimator is asymptotically efficient under R0 ~ R5 Regularity Conditions

Asymptotic Confidence Interval

By the asymptotic normality of MLE, $\sqrt{nI(\hat{\theta })}(\hat{\theta }-\theta )\stackrel{D}{\to }N(0,1)$ Thus, $100(1-\alpha )\%$ confidence interval of for $\theta$ is $\left(\hat{\theta }-z_{\alpha /2}\frac{1}{\sqrt{nI(\hat{\theta })}},\hat{\theta }+z_{\alpha /2}\frac{1}{\sqrt{nI(\hat{\theta })}}\right)$

Delta method for MLE Estimator

Under the R0 ~ R5 Regularity Conditions, let $g(x)$ be a continuous function and $g^{\prime }({\theta }_0)\ne 0$, then $\sqrt{n}(g({\hat{\theta }}_n)-g({\theta }_0))\stackrel{D}{\to }N\left(0,\frac{g^{\prime }({\theta }_0{)}^2}{I({\theta }_0)}\right)$

Facts

Under R0 and R1 regularity conditions, let ${\theta }_0$ be a true parameter, then $\forall \theta \ne {\theta }_0,\underset{n\to \infty }{\lim }{P}_{{\theta }_0}[L({\theta }_0)>L(\theta )]=1$

Link to original

Goodness-of-Fit of the Regression Line

Decomposition of Sum of Squares

Definition

Simple Linear Regression Case

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The total sum of squares $SST={\sum }_{i=1}^n(Y_i-\bar{Y}{)}^2$ is expressed by the sum of the unexplained variation $SSE$ and the explained variation $SSR$ $SST=SSE+SSR$ where

• $SST={\sum }_{i=1}^n(Y_i-\bar{Y}{)}^2$ is the total sum of squares and has the degree of freedom $(n-1)$,
• $SSE={\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2$ is the error sum of squares and has the degree of freedom $(n-2)$, and
• $SSR={\sum }_{i=1}^n({\hat{Y}}_i-\bar{Y}{)}^2$ is the regression sum of squares and has the degree of freedom $1$

Multiple Linear Regression Case

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$

• $SST={\sum }_{i=1}^n(Y_i-\bar{Y}{)}^2$ has the degree of freedom $(n-1)$,
• $SSE={\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2$ has the degree of freedom $(n-p)$, and
• $SSR={\sum }_{i=1}^n({\hat{Y}}_i-\bar{Y}{)}^2$ has the degree of freedom $p-1$

Facts

$SST/{\sigma }^2\sim {\chi }^2\left((n-1),\frac{1}{2{\sigma }^2}{\mathbf{\beta }}^{⊺}{\mathbf{X}}^{⊺}(I-\frac{1}{n}J)\mathbf{X}\mathbf{\beta }\right)$ $SSE/{\sigma }^2\sim {\chi }^2(n-p)$ $SSR/{\sigma }^2\sim {\chi }^2\left((p-1),\frac{1}{2{\sigma }^2}{\mathbf{\beta }}^{⊺}{\mathbf{X}}^{⊺}(I-\frac{1}{n}J)\mathbf{X}\mathbf{\beta }\right)$

Link to original

Coefficient of Determination

Definition

The coefficient of determination of the linear regression model is defined as $R^2=\frac{SSR}{SST}=\frac{{\sum }_{i=1}^n({\hat{Y}}_i-\bar{Y}{)}^2}{{\sum }_{i=1}^n(Y_i-\bar{Y}{)}^2}$

Facts

The coefficient of determination is a squared value of Pearson Correlation Coefficient

Link to original

Analysis of Variance

One-way ANOVA

One-way ANOVA is used to analyze the significance of differences of means between groups. Let an $i$-th response of the $j$-th group be $X_{ij}={\mu }_j+{ϵ}_{ij}$ where ${ϵ}_{ij}\sim N(0,{\sigma }^2),\forall i=1,2,\dots ,n_j,\forall j=1,2,\dots ,c$

We want to test the null hypothesis $H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_c$ i.e. there is no treatment effect. $F=\frac{\sum _{j=1}^c{n}_j({\bar{X}}_{.j}-{\bar{X}}_{..}{)}^2/(c-1)}{\sum _{j=1}^c\sum _{i=1}^{n_j}(X_{ij}-{\bar{X}}_{.j}{)}^2/(N-c)}=\frac{SSB/(c-1)}{SSW/(N-c)}\sim F(c-1,N-c)$ where $N=\sum _{j=1}^c{n}_j$ is the total number of observations, ${\bar{X}}_{.j}$ represents the mean of the $j$-th group and ${\bar{X}}_{..}$ represents the overall mean, $\frac{1}{{\sigma }^2}\sum _{j=1}^c{n}_j({\bar{X}}_{.j}-{\bar{X}}_{..}{)}^2\sim {\chi }^2(c-1)$ of the numerator indicates between-group variance (SSB) and $\frac{1}{{\sigma }^2}\sum _{j=1}^c\sum _{i=1}^{n_j}(X_{ij}-{\bar{X}}_{.j}{)}^2\sim {\chi }^2(N-c)$ of the denominator indicates within-group variance (SSW)

The Likelihood Ratio Test rejects $H_0$ if $F>F_{\alpha }(c-1,N-c)$

Link to original

Statistical Inference

Confidence Interval for the Mean in Simple Linear Regression

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The $100(1-\alpha )\%$ confidence region for the mean response $E(Y)$, when $x$ is given, is defined as $({\hat{\beta }}_0+{\hat{\beta }}_{1}x)\pm t_{\alpha /2}(n-2)\times s\sqrt{\frac{1}{n}+\frac{(x-\bar{X}{)}^2}{S_{XX}}}$ where $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-2}$

Link to original

Prediction interval for a New Response in Linear Regression

Definition

Consider a Simple Linear Regression model $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ The prediction interval for a new response $Y(x)$, when $x$ is given, is defined as $\hat{Y}(x)\pm t_{\alpha /2}((n-2))\times s\sqrt{1+\frac{1}{n}+\frac{(x-\bar{X}{)}^2}{S_{XX}}}$ where $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-2}$

Facts

The prediction interval is always wider than the Confidence Interval.

Link to original

Multiple Linear Regression Model

Multiple Linear Regression Model

Multiple Linear Regression

Definition

$Y_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_{p-1}{X}_{i,p-1}+{ϵ}_i,i=1,2,\dots ,n$ where ${ϵ}_i$‘s are i.i.d. error terms, with $E({ϵ}_i)=0$ and $\mathrm{Var}({ϵ}_i)={\sigma }^2$

Matrix Notations

Let $\mathbf{\beta }=({\beta }_0,{\beta }_1,{\beta }_2,\dots ,{\beta }_{p-1}{)}^{⊺}$ and ${\mathbf{x}}_i=(1,X_{i1},X_{i2},\dots ,X_{i,p-1}{)}^{⊺}$, then the regression model can be express the model as $Y_i={\mathbf{x}}_i^{⊺}\mathbf{\beta }+{ϵ}_i,i=1,2,\dots ,n$

Let $\mathbf{y}=[Y_1,Y_2,\dots ,Y_n{]}^{⊺}$, $\mathbf{X}=[{\mathbf{x}}_1^{⊺},{\mathbf{x}}_2^{⊺},\dots ,{\mathbf{x}}_n^{⊺}{]}^{⊺}$, $\mathbf{ϵ}=[{ϵ}_1,{ϵ}_2,\dots ,{ϵ}_n{]}^{⊺}$, then the regression model can be express as $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ where $E[\mathbf{ϵ}]=0$, and $\mathrm{Cov}[\mathbf{ϵ}]={\sigma }^2{I}_n$

Link to original

Statistical Inference

Distribution of Regression Coefficient

Definition

Distribution of Regression Coefficient

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The Least Squares Estimator $\hat{\mathbf{\beta }}=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$ also follows Normal Distribution. $\hat{\mathbf{\beta }}\sim N_p(\mathbf{\beta },({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\sigma }^2)$ where $p-1$ is the number of explanatory variables.

Marginal Distribution of Regression Coefficient

The marginal distribution of the Multivariate Normal Distribution is univariate Normal Distribution ${\hat{\beta }}_j=N({\beta }_j,c_{j+1}{\sigma }^2)$ where $c_{j+1}=\mathrm{diag}[({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{]}_{j+1}$

Link to original

Confidence Interval for Regression Coefficient

Definition

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The $100(1-\alpha )\%$ confidence interval for ${\beta }_j$ is defined as ${\hat{\beta }}_j\pm t_{\alpha /2}(n-p)SE({\hat{\beta }}_j)$ where $SE({\hat{\beta }}_j)$ is the Standard Error for Regression Coefficient ${\beta }_j$, and $p-1$ is the number of explanatory variables.

Link to original

Standard Error for Regression Coefficient

Definition

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The standard error for regression coefficient ${\beta }_j$ is calculated as $SE({\hat{\beta }}_j)=\sqrt{c_{j+1}}s$ where $c_{j+1}=\mathrm{diag}[({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{]}_{j+1}$, $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-p}$, and $p-1$ is the number of explanatory variables.

Link to original

Hypothesis testing for Regression Coefficient

Definition

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The null hypothesis $H_0:{\beta }_j={\beta }_{j,0}$ can be tested with the test statistic $t_j=\frac{{\hat{\beta }}_j-{\beta }_{j,0}}{SE({\hat{\beta }}_j)}\sim t(n-p)$ where $SE({\hat{\beta }}_j)$ is the Standard Error for Regression Coefficient ${\beta }_j$.

Link to original

Joint Confidence Region for Regression Coefficient

Definition

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The $100(1-\alpha )\%$ join confidence region for $\mathbf{\beta }$ is defined as $\{\mathbf{\beta }|(\mathbf{\beta }-\hat{\mathbf{\beta }}{)}^{⊺}{\mathbf{X}}^{⊺}\mathbf{X}(\mathbf{\beta }-\hat{\mathbf{\beta }})\le p{s}^2{F}_{\alpha }(p,n-p)\}$ where $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-p}$ and $p-1$ is the number of explanatory variables.

The joint confidence region is ellipsoidal shape whose center is $\hat{\mathbf{\beta }}$ by the inequality.

Link to original

Simultaneous Confidence Interval for Regression Coefficient

Definition

Simultaneous confidence interval is used when computing confidence intervals for $g$ parameters ${\beta }_1,{\beta }_2,\dots ,{\beta }_g$ simultaneously.

Joint confidence region provides the accurate elliptical area as the confidence region, but the calculations and interpretations are complex. To obtain a rectangular-shaped confidence interval, the simultaneous confidence interval is used. When multiplying the confidence intervals of each coefficient, the resulting confidence region becomes smaller than the desired area. Therefore, to obtain the “at least $100(1-\alpha )\%$” confidence region, correction method is used. Conservative methods like the Bonferroni’s Method, provide a confidence region much larger than our desired $100(1-\alpha )\%$, which satisfies the condition of being at least $100(1-\alpha )\%$ but reduces the power of the test. To address this issue, other methods have been devised.

Bonferroni’s Method

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The $100(1-\alpha )\%$ Bonferroni confidence interval for ${\beta }_j$ is defined as ${\hat{\beta }}_j\pm t_{\alpha /2g}(n-p)SE({\hat{\beta }}_j)$ where $SE({\hat{\beta }}_j)$ is the Standard Error for Regression Coefficient ${\beta }_j$ and $p-1$ is the number of explanatory variables.

Scheffe Method

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The $100(1-\alpha )\%$ Scheffe confidence interval for ${\beta }_j$ is defined as ${\hat{\beta }}_j\pm \sqrt{g{F}_{\alpha }(g,n-p)}SE({\hat{\beta }}_j)$ where $SE({\hat{\beta }}_j)$ is the Standard Error for Regression Coefficient ${\beta }_j$ and $p-1$ is the number of explanatory variables.

Maximum Modulus Method

Assume that $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\mathbf{I}}_n{\sigma }^2)$ The $100(1-\alpha )\%$ Maximum Modulus confidence interval for ${\beta }_j$ is defined as ${\hat{\beta }}_j\pm u_{\alpha }(g,n-p)SE({\hat{\beta }}_j)$ where $u$ is a Random Variable representing the maximum of the absolute of $g$ independent random variables following $t(n-p)$, $SE({\hat{\beta }}_j)$ is the Standard Error for Regression Coefficient ${\hat{\beta }}_j$, and $p-1$ is the number of explanatory variables.

Facts

Maximum Modulus method < Bonferroni’s method < Scheffe’s method Here, the second inequality holds only when the degree of freedom of comparisons is relatively small compared to the number of groups being compared.

Link to original

Confidence Interval for the Mean in Multiple Linear Regression

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ The $100(1-\alpha )\%$ confidence region for the mean response $E(Y)$, when ${\mathbf{x}}_0$ is given, is defined as ${\mathbf{x}}_0^{⊺}\mathbf{\beta }\pm t_{\alpha /2}(n-p)\times s\sqrt{{\mathbf{x}}_0^{⊺}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{x}}_0}$ where $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-p}$ and $p-1$ is the number of explanatory variables

Link to original

Partial F-Test

Definition

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ and the hypothesis $H_0:{\beta }_r={\beta }_{r+1}=\cdots ={\beta }_{p-1}=0$

Define the full model as $Y_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_{p-1}{X}_{i,p-1}+{ϵ}_i,i=1,2,\dots ,n$ and under $H_0$, define the reduced model as $Y_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_{r-1}{X}_{i,r-1}+{ϵ}_i,i=1,2,\dots ,n$

We reject $H_0$ if $SSR(F)-SSR(R)$ is large, where $SSR(F)$ and $SSR(R)$ are the regression sum of squares of the full and reduced model, respectively. By the equality $SST=SSE(F)+SSR(F)=SSE(R)+SSR(R)$, we can test the $SSE(R)-SSE(F)$ instead of the $SSR(F)-SSR(R)$

$\frac{SSE(R)-SSE(F)}{{\sigma }^2}\sim {\chi }^2(p-r)$ and $\frac{SSE(F)}{{\sigma }^2}\sim {\chi }^2(n-p)\sim {\chi }^2(n-p)$. Therefore, $F=\frac{(SSE(R)-SSE(F))/({\text{df}}_R-{\text{df}}_F)}{SSE(F)/{\text{df}}_F}\sim F(p-r,n-p)$ where ${\text{df}}_R=n-r$ and ${\text{df}}_F=n-p$

Hence, we reject $H_0$ if $F>F_{\alpha }(p-r,n-p)$ and we call it partial F-test

Link to original

The General Linear Test

Definition

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ The general linear test can be written as $H_0:\mathbf{C}\mathbf{\beta }=\mathbf{m},H_1:\mathbf{C}\mathbf{\beta }\ne \mathbf{m}$ where $\mathbf{C}$ is $q\times p$ matrix with rank $q$ and $\mathbf{m}$ is $q\times 1$ vector.

e.g. for the hypothesis $H_0:{\beta }_1-{\beta }_2=0,{\beta }_3+2{\beta }_4=3$, $C=\left[\begin{array}{ccccc}0& 1& -1& 0& 0\\ 0& 0& 0& 1& 2\end{array}\right],\mathbf{m}=[0,3{]}^{⊺}$

${\text{df}}_R=n-(p-q)$ and ${\text{df}}_F=(n-p)$. Therefore, $F=\frac{SSE(R)-SSE(F)}{q\cdot MSE(F)}\sim F(q,n-p)$ where $MSE(F)=SSE(F)/{\text{df}}_F$

Hence, we reject $H_0$ if $F>F_{\alpha }(q,n-p)$

Method of Lagrange Multiplier

By the method of Lagrange multiplier, $SSE(R)-SSE(F)=(\mathbf{C}\hat{\mathbf{\beta }}-\mathbf{m}{)}^{⊺}[\mathbf{C}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{C}}^{⊺}{]}^{-1}(\mathbf{C}\hat{\mathbf{\beta }}-\mathbf{m})$ Therefore, we can calculate the F-statistic without fitting the reduced model $F=\frac{(\mathbf{C}\hat{\mathbf{\beta }}-\mathbf{m}{)}^{⊺}[\mathbf{C}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{C}}^{⊺}{]}^{-1}(\mathbf{C}\hat{\mathbf{\beta }}-\mathbf{m})}{q\cdot MSE(F)}\sim F(q,n-p)$

Link to original

Lack of Fit Test

Lack of Fit Test

Definition

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ We want to test the linearity of the model $H_0:E(Y_i)={\mathbf{X}}_i^{⊺}\mathbf{\beta },i=1,2,\dots ,n$

Assume that $Y_{il}(i=1,2,\dots ,m;l=1,2,\dots ,n_i)$ is the $l$-th replication at $x_i$, then The $SSE$ is decomposed to the sum of the pure error sum of squares $SSPE$ and the lack-of-fit sum of squares $SSLF$. $\underset{SSE}{\underbrace{\sum _{i=1}^m\sum _{l=1}^{n_i}(Y_{il}-{\hat{Y}}_i{)}^2}}=\underset{SSPE}{\underbrace{\sum _{i=1}^m\sum _{l=1}^{n_i}(Y_{il}-{\bar{Y}}_i{)}^2}}+\underset{SSLF}{\underbrace{\sum _{i=1}^m\sum _{l=1}^{n_i}({\bar{Y}}_i-{\hat{Y}}_i{)}^2}}$

The degree of freedoms are ${\text{df}}_{SSE}=n-p$, ${\text{df}}_{SSPE}=n-m$, and ${\text{df}}_{SSLF}=m-p$, where $n:=\sum _{i=1}^m{n}_i$

Also, the mean squares are obtained by dividing the sum of squares by its degree of freedom $M{S}_{PE}=\frac{{\sum }_{i=1}^m{\sum }_{l=1}^{n_i}(Y_{il}-{\bar{Y}}_i{)}^2}{n-m},M{S}_{LF}=\frac{{\sum }_{i=1}^m{\sum }_{l=1}^{n_i}({\bar{Y}}_i-{\hat{Y}}_i{)}^2}{m-p}$

Under $H_0$, the $M{S}_{LF}$ is an unbiased estimator of ${\sigma }^2$, and $M{S}_{PE}$ is always an unbiased estimator of ${\sigma }^2$, regardless of $H_0$. Therefore, we can test the hypothesis by comparing ratios of both estimators.

$\frac{S{S}_{LF}}{{\sigma }^2}\sim {\chi }^2(n-p)$ and $\frac{S{S}_{PE}}{{\sigma }^2}={\chi }^2(n-m)$. Therefore, $F_{L}F=\frac{M{S}_{LF}}{M{S}_{PE}}\sim F(m-p,n-m)$

Hence, we reject $H_0$ if $F>F_{\alpha }(m-p,n-m)$

Facts

The lack of fit test uses data’s replication, Therefore, there should be more than one data (replication) for some $x_i$‘s.

Link to original

Miscellanea

Multicollinearity

Definition

If two or more covariates are highly correlated, in a multiple regression model, the model has a multicollinearity The multicollinearity makes the estimation of coefficient very unstable. So, the estimated coefficients are not reliable.

Link to original

Regression Diagnostics

Residual

Studentized Residual

Definition

A studentized residual is a technique used to detect outliers

Internally Studentized Residual

$r_i=\frac{e_i}{s\sqrt{1-h_{ii}}}$ where $s=\frac{1}{n-p}\sum _{i=1}^n{\widehat{ϵ}}_i^2$

$\frac{r_i^2}{n-p}\sim \mathrm{Beta}(\frac{1}{2},\frac{n-p-1}{2})$

Externally Studentized Residual

$r_i^{\ast }=\frac{e_i}{s_{(i)}\sqrt{1-h_{ii}}}\sim (n-p-1)$ where $s_{(i)}^2=\frac{1}{n-p}\sum _{j=1}^n{\widehat{ϵ}}_j^2\cdot 1(j\ne i)$ is the unbiased estimator of ${\sigma }^2$ based on $n-1$ observations after deleting the $i$-th observation

Facts

$s_{(i)}^2=s^2\frac{n-p-r_i^2}{n-p-1}$ $r_i^{\ast }=r_i\sqrt{\frac{n-p-1}{n-p-r_i^2}}$

Link to original

Leverage (Statistics)

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $h_{ii}={\mathbf{x}}_i^{⊺}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{x}}_i$

$h_{ii}$ is the $i$-th leverage, the distance between ${\mathbf{x}}_i$ and $\bar{\mathbf{x}}$

A measure of how far away the independent variable value of an observation is from those of the other observations It is used to detect outlier

Facts

$\sum _{i=1}^n{h}_{ii}=p$, where $p-1$ is the number of explanatory variables $\frac{1}{n}\le h_{ii}\le 1$ $h_{ij}^2\le h_{ii}{h}_{jj}$

Link to original

Influence Measures

Cook's Distance

Definition

Single Observation

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $C_i=\frac{1}{p}(\hat{\mathbf{\beta }}-{\hat{\mathbf{\beta }}}_{(i)}{)}^{⊺}\mathrm{Cov}(\hat{\mathbf{\beta }}{)}^{-1}(\hat{\mathbf{\beta }}-{\hat{\mathbf{\beta }}}_{(i)})=\frac{1}{p{\sigma }^2}\frac{e_i^2{h}_{ii}}{(1-h_{ii}{)}^2}$ where $p-1$ is the number of explanatory variables, and $h_{ii}$ is the Leverage (Statistics)

$C_i$ is the Cook’s distance. The more influential the data point, the larger the Cook’s distance

Set of Observations

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $C_K=\frac{1}{p}(\hat{\mathbf{\beta }}-{\hat{\mathbf{\beta }}}_{(K)}{)}^{⊺}\mathrm{Cov}(\hat{\mathbf{\beta }}{)}^{-1}(\hat{\mathbf{\beta }}-{\hat{\mathbf{\beta }}}_{(K)})=\frac{1}{p{\sigma }^2}{\mathbf{e}}_K^{⊺}(I-{\mathbf{H}}_K{)}^{-1}{\mathbf{H}}_K(I-{\mathbf{H}}_K{)}^{-1}{\mathbf{e}}_K$ where $p-1$ is the number of explanatory variables, and $\mathbf{H}$ is the matrix of leverages

Link to original

Andrews-Pregibon Statistic

Definition

Single Observation

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $A{P}_i=(1-h_{ii})\frac{1-e_i^2}{(1-h_{ii}){\mathbf{e}}^{⊺}\mathbf{e}}=1-h_{ii}^{\ast }$ where $h_{ii}^{\ast }={\mathbf{x}}_i^{\ast ⊺}({\mathbf{X}}^{\ast ⊺}{\mathbf{X}}^{\ast }{)}^{-1}{\mathbf{x}}_i^{\ast }$ is the Leverage (Statistics) of the matrix ${\mathbf{X}}^{\ast }$

$A{P}_i$ is the Andrews-Pregibon statistic. The more influential the data point, the smaller the Andrews-Pregibon statistic

Set of Observations

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $A{P}_K=\det [I-{\mathbf{H}}_K^{\ast }]$ where ${\mathbf{H}}_K^{\ast }={\mathbf{X}}_K^{\ast }({\mathbf{X}}^{\ast ⊺}{\mathbf{X}}^{\ast }{)}^{-1}{\mathbf{X}}_K^{\ast ⊺}$ is the matrix of leverages

Link to original

DFBETAS

Definition

The DFBEATS consider the influence of the $i$-th observation on ${\hat{\mathbf{\beta }}}_j$

Single Observation

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $DFBETA{S}_{j,i}=\frac{{\hat{\mathbf{\beta }}}_j-{\hat{\mathbf{\beta }}}_{j(i)}}{s_{(i)}\sqrt{({\mathbf{X}}^{⊺}\mathbf{X}{)}_{jj}^{-1}}}=\frac{c_{ji}}{\sqrt{{\sum }_{l=1}^n{c}_{jl}^2}}\frac{e_i}{s_{(i)}(1-h_{ii})}$ where $s_{(i)}$ is the externally studentized residual and $c_{ji}$ is the element of the matrix $C:=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}$

Set of Observations

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $DFBETA{S}_{j,K}=\frac{{\hat{\mathbf{\beta }}}_j-{\hat{\mathbf{\beta }}}_{j(K)}}{s_{(K)}\sqrt{({\mathbf{X}}^{⊺}\mathbf{X}{)}_{jj}^{-1}}}=\frac{{\mathbf{c}}_{jk}^{⊺}(I-{\mathbf{H}}_K{)}^{-1}{\mathbf{e}}_K}{s_{(K)}\sqrt{{\sum }_{l=1}^n{c}_{jl}^2}}$ where $\mathbf{H}$ is the matrix of leverages and $c_{ji}$ is the element of the matrix $C:=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}$

Link to original

DFFITS

Definition

The DFFITS consider the influence of the $i$-th observation on $\widehat{y_i}$

Single Observation

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $DFFIT{S}_i=\frac{|{\mathbf{x}}_i^{⊺}(\hat{\mathbf{\beta }}-{\hat{\mathbf{\beta }}}_{(i)})|}{s_{(i)}\sqrt{h_{ii}}}=\frac{1}{s_{(i)}}\frac{|e_i|\sqrt{h_{ii}}}{1-h_{ii}}=\frac{|r_i^{\ast }|\sqrt{h_{ii}}}{\sqrt{1-h_{ii}}}$ where $s_{(i)}$ is the externally studentized residual

Set of Observations

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $DFFIT{S}_K=\frac{1}{s_{(K)}}\frac{{\mathbf{x}}_i^{⊺}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}_K^{⊺}(I-{\mathbf{H}}_K{)}^{-1}{\mathbf{e}}_K}{\sqrt{h_{ii}}}$ where $s_{(K)}$ is the externally studentized residual and $\mathbf{H}$ is the matrix of leverages

Link to original

COVRATIO

Definition

The COVRATIO consider the influence of the $i$-th observation on $\mathrm{Cov}(\hat{\mathbf{\beta }})$

Single Observation

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $COVRATI{O}_i=\frac{\det (s_{(i)}^2({\mathbf{X}}_{(i)}^{⊺}{\mathbf{X}}_{(i)}{)}^{-1})}{\det (s^2({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1})}$

Set of Observations

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $COVRATI{O}_K=\frac{\det (s_{(K)}^2({\mathbf{X}}_{(K)}^{⊺}{\mathbf{X}}_{(K)}{)}^{-1})}{\det (s^2({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1})}$

Link to original

Multicollinearity

Variance Inflation Factor

Definition

Consider a Multiple Linear Regression model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ $VI{F}_k=\frac{1}{1-R_k^2},k=1,2,\dots ,p-1$ where $p-1$ is the number of explanatory variables.

The variance inflation factor (VIF) measures how much the variance of an estimated regression coefficient is increased because of Multicollinearity

Link to original

Condition Number

Definition

Consider a matrix $X$ with singular values $d_1\ge d_2\ge \cdots \ge d_p\ge 0$ $\kappa (X)=\frac{d_1}{d_p}$

$\kappa (X)$ is the conditional number of $X$

Facts

In linear regression, if the conditional value of the matrix $X^{⊺}X$ is large, where $X$ is a design matrix, the model has a Multicollinearity problem.

Link to original

Autocorrelation and Durbin-Watson Test

Durbin-Watson Test

Definition

Test used to detect the presence of autocorrelation at lag $1$ in the residual

Assume $e_i$ is the residual structure of $e_t=\rho e_{t-1}+{\nu }_t$ $DW=\frac{{\sum }_{t=2}^T(e_t-e_{t-1}{)}^2}{{\sum }_{t=1}^T{e}_t^2}\simeq 2(1-\hat{\rho })$ where $T$ is the number of observations and $\hat{\rho }=\frac{\sum _{t=2}^T{e}_t{e}_{t-1}}{\sum _{t=2}^T{e}_{t-1}^2}$

$DW=2$ indicates no autocorrelation

Link to original

Selection of Regression Models

Adjusted R_squared Value

Definition

$R_a^2=1-\frac{SS{E}_p/(n-p)}{SST/(n-1)}=1-\frac{s_p^2}{SST/(n-1)}$ where $p-1$ is the number of explanatory variables.

Link to original

Mallow's Cp

Definition

$C_p=\frac{SS{E}_p}{s^2}-(n-2p)=\overline{\text{err}}+\frac{p}{n}{\hat{\sigma }}^2$ where $p-1$ is the number of explanatory variables, $s^2$ is a sample variance under the full model, and $\overline{\text{err}}=\frac{RSS}{n}$

Mallow’s $C_p$ is used to assess the fit of a regression model. A small value of $C_p$ means that the model is relatively precise.

Link to original

PRESS statistic

Definition

$PRES{S}_p=\sum _{i=1}^n(y_i-{\hat{y}}_{i(i)}{)}^2=\sum _{i=1}^n{\left(\frac{e_i}{1-h_{ii}}\right)}^2$ where $y_{i(i)}$ is an estimated value of $y_i$ calculated by the $n-1$ data excluding $y_i$

The predicted residual error sum of squares (PRESS) is a form of cross validation used in regression analysis

Link to original

Model Validation

Cross Validation

Definition

Partition a set of data into $K$ sets $I_j,j=1,2,\dots ,K$ and denote a function $\kappa :\{1,2,\dots ,n\}\to \{1,2,\dots ,K\}$ such that $\kappa (i)=j\iff i\in I_j$. For the dataset $\{({\mathbf{X}}_1,Y_1),({\mathbf{X}}_2,Y_2),\dots ,({\mathbf{X}}_n,Y_n)\}$, let ${\hat{f}}_{-j}$ be the estimator based on observations except $I_j$, then the cross validation estimation for the prediction error is defined as ${\widehat{PE}}_{CV}(\hat{f})=\frac{1}{n}\sum _{i=1}^{n}L(Y_i,{\hat{f}}_{-\kappa (i)}({\mathbf{X}}_i))$ where $L$ is a loss function

Facts

If the data is partitioned into $k$ group with equal size, then it is called a $k$-fold cross validation.

Link to original

Miscellanea in Model Selection

Jensen's Inequality

Definition

Let $f$ be a Convex Function on an interval $I$, $X$ be a Random Variable with support $S\subset I$, and $E(X)<\infty$, then $f(E[X])\le E[f(X)]$

Facts

$\frac{1}{\frac{1}{n}\sum \frac{1}{a_i}}\le (\prod _{i=1}^n{a}_i{)}^{1/n}\le \frac{1}{n}\sum _{i=1}^n{a}_i$

By Jensen’s Inequality, the following relation is satisfied arithmetic mean $\le$ geometric mean $\le$ harmonic mean

Link to original

Kullback-Leibler Divergence

Definition

Assume that two Probability Distributions $P$ and $Q$ are given. Then the Kullback-Leibler divergence between $P$ and $Q$ is defined as $D_{KL}(P||Q)=E_{x\sim P}\left[\ln \left(\frac{p(x)}{q(x)}\right)\right]={\int }_{x}p(x)\ln \frac{p(x)}{q(x)}dx$

Kullback-Leibler divergence measures how different two distributions are.

It also can be expressed as a difference between the cross entropy $H(p,q)$ (difference between distributions $p$ and $q$)and entropy $H(p)$ (inherent uncertainty of $p$). $H(P,Q)-H(P)={\int }_{x}p(x)\ln \frac{1}{q(x)}dx+{\int }_{x}p(x)\ln p(x)dx={\int }_{x}p(x)\ln \frac{p(x)}{q(x)}dx=D_{KL}(P\Vert Q)$

Facts

Let $\{P_n{\}}_{n\in \mathbb{N}}$ be a sequence of distributions. Then, $D_{KL}(P_n||P)\to 0⟹D_{JS}(P_n,P)\to 0⟺\delta (P_n,P)\to 0⟹W(P_n,P)\to 0⟺P_n\stackrel{D}{\to }P$ The convergence of the KL-Divergence to zero implies that the JS-Divergence also converges to zero. The convergence of the JS-Divergence to zero is equivalent to the convergence of the Total Variation Distance to zero. The convergence of the Total Variation Distance to zero implies that the Wasserstein Distance also converges to zero. The convergence of the Wasserstein Distance to zero is equivalent to the Convergence in Distribution of the sequence.

Link to original

Akaike Information Criterion

Definition

$AIC=-\frac{2}{n}\sum _{i=1}^n\ln f(X_i\hat{\mathbf{\theta }})+\frac{2p}{n}$ where $\hat{\mathbf{\theta }}$ is the MLE of the postulated model’s parameter $\mathbf{\theta }$

Link to original

Bayesian Information Criterion

Definition

$BIC=\frac{SS{E}_p}{n{\sigma }^2}+\frac{p\ln n}{n}$

Facts

BIC more penalize than AIC when $p$ is large

Link to original

Transformation of the Linear Regression Model

The Use of Dummy Variables

Dummy Variable

Definition

$D=\{\begin{array}{ll}1& X\in A\\ 0& \text{otherwisse}\end{array}$

A dummy variable is one that takes a binary variable. It is commonly used in regression analysis to represent categorical variables that have more than two categories.

Link to original

One-way ANOVA with a Regression Model

$Y={\beta }_0+{\beta }_1{D}_1+{\beta }_2{D}_2+\cdots +{\beta }_{c-1}{D}_{c-1}+ϵ$ where $D_i$ are dummy variables represent each category, and $c$ is the number of categories

In the setting, a corner-point constraint is used.

The null hypothesis $H_0:{\beta }_1={\beta }_2=\cdots ={\beta }_{c-1}=0$ i.e. there is no treatment effect, yields same test statistic as the null hypothesis of the one-way ANOVA with equal replication $H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_c$

Link to original

Two-way ANOVA with a Regression Model

$Y={\beta }_0+\sum _{i=1}^{a-1}{\alpha }_i{D}_i+\sum _{j=1}^{b-1}{\beta }_j{E}_j+\sum _{i=1}^{a-1}\sum _{j=1}^{b-1}{\gamma }_{ij}{D}_i{E}_j+ϵ$ where $D_i$ and $E_j$ are dummy variables representing categories of the two factors, $a$ is the number of categories for the first factor, and $b$ is the number of categories for the second factor

In this setting, a corner-point constraint is used for both factors.

The null hypotheses can be tested by Deviance

Link to original

Polynomial Regression

Polynomial Regression

Definition

Consider explanatory variables $X_i$ and a response variable $Y$ The $k$-variables 2nd order polynomial regression model is defined as $Y={\beta }_0+\sum _{i=1}^k{\beta }_i{X}_i+\underset{\text{quadratic terms}}{\underbrace{\sum _{i=1}^k{\beta }_{ii}{X}_i^2}}+\underset{\text{interaction terms}}{\underbrace{\sum _{i<j}^k{\beta }_{ij}{X}_i{X}_j}}+ϵ$

The single variable $k$-th order polynomial regression model is defined as $Y={\beta }_0+\sum _{i=1}^k{\beta }_i{X}^iϵ$

Facts

The column vectors of the design matrix $\mathbf{X}$ in the polynomial regression are highly correlated, so proper transformation is needed. Among them, orthogonal polynomial is often used.

Link to original

Response Surface Analysis

Definition

The goal of response surface analysis is finding the optimal condition for the response. It employs Polynomial Regression model to model response surface.

Consider explanatory variables $X_i$ and a response variable $Y$. The second-order response surface model is defined as $Y={\beta }_0+{\mathbf{x}}^{⊺}\mathbf{\beta }+{\mathbf{x}}^{⊺}\mathbf{Bx}+ϵ$ where $\mathbf{x}=(X_1,X_2,\dots ,X_k{)}^{⊺}$, $\mathbf{\beta }=({\beta }_1,{\beta }_2,\dots ,{\beta }_k{)}^{⊺}$, $\mathbf{B}=\left[\begin{array}{cccc}{\beta }_{11}& \frac{{\beta }_{12}}{2}& \dots & \frac{{\beta }_{1k}}{2}\\ \frac{{\beta }_{12}}{2}& {\beta }_{22}& \dots & \frac{{\beta }_{2k}}{2}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\beta }_{1k}}{2}& \frac{{\beta }_{2k}}{2}& \dots & {\beta }_{kk}\end{array}\right]$ Let the fitted model is $\hat{Y}={\hat{\beta }}_0+{\mathbf{x}}^{⊺}\hat{\mathbf{\beta }}+{\mathbf{x}}^{⊺}\hat{\mathbf{B}}\mathbf{x}$. Then, the stationary point is obtained by ${\mathbf{x}}_s=-\frac{1}{2}{\hat{\mathbf{B}}}^{-1}\hat{\mathbf{\beta }}$ and the fitted value at the stationary point is ${\hat{Y}}_s={\hat{\beta }}_0+{\mathbf{x}}_s^{⊺}\hat{\mathbf{\beta }}+{\mathbf{x}}_s^{⊺}\hat{\mathbf{B}}{\mathbf{x}}_s$. If $\hat{\mathbf{B}}$ is Positive-Definite Matrix, then ${\hat{Y}}_s$ is minimum If $\hat{\mathbf{B}}$ is Negative-Definite Matrix, then ${\hat{Y}}_s$ is maximum If $\hat{\mathbf{B}}$ is neither Positive-Definite Matrix nor Negative-Definite Matrix, then ${\hat{Y}}_s$ is saddle point.

Link to original

Weighted Least Squares Method

Weighted Least Squares

Definition

Weighted least squares (WLS) is a generalization of Least Square to cope with Heteroskedasticity.

We assume that $\mathrm{Cov}(\mathbf{ϵ})={\sigma }^2{\mathbf{W}}^{-1}$ where $\mathbf{W}=\mathrm{diag}(w_1,w_2,\dots ,w_n)$ The weight matrix is positive definite, so there exists a non-singular matrix $\mathbf{P}:=W^{-1/2}$. Consider a transformed model $\mathbf{Py}=\mathbf{PX}\mathbf{\beta }+\mathbf{P}\mathbf{ϵ}$ Now, the new error term $\mathbf{P}\mathbf{ϵ}$ is i.i.d. $\mathrm{Cov}(\mathbf{P}\mathbf{ϵ})={\sigma }^2\mathbf{I}$ And the weighted least squares estimator is obtained by ${\hat{\mathbf{\beta }}}_{WLS}=({\mathbf{X}}^{⊺}\mathbf{WX}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{Wy}$

Link to original

Box-Cox Transformation Model

Box-Cox Transformation Model

Definition

Box-Cox transformation is useful when dealing with heteroskedasticity by stabilizing variance.

Box-Cox Transformation

$w_i=\{\begin{array}{ll}\frac{y_i^{\lambda }-1}{\lambda }& \lambda \ne 0\\ \ln y_i& \lambda =0\end{array}$

Box-Cox Transformation Model

The model using the Box-Cox Transformed variable as a response variable is called a Box-Cox transformation model $\mathbf{w}(\lambda )=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ where $\mathbf{ϵ}\sim N(0,{\sigma }^2\mathbf{I})$

Link to original

Robust Regression

M-estimator

Definition

Consider an error (loss) function $\rho (\cdot )$, The M-estimator is obtained by ${\hat{\mathbf{\beta }}}_M=\underset{\mathbf{\beta }}{\mathrm{argmin}}\sum _{i=1}^n\rho \left(\frac{{ϵ}_i}{s}\right)=\underset{\mathbf{\beta }}{\mathrm{argmin}}\sum _{i=1}^n\rho \left(\frac{y_i-{\mathbf{x}}_i^{⊺}\mathbf{\beta }}{s}\right)$ where $s$ is a robust estimator, and given by $s=\frac{\mathrm{med}(|e_i-\mathrm{med}(e_i)|)}{0.6745}$

If $\rho (u)=\frac{1}{2}{u}^2$, then it is LSE, and if $\rho (u)=|u|$, then it is $L_1$-norm regression estimator. The derivative of $\rho$, denoted $\psi =\frac{\partial \rho }{\partial \mathbf{\beta }}$, is called the influence function. Since $\psi$ is a non-liner in $\beta$, we van not get explicit solution. Instead, we use an iterative method, called IRLS (Iterative Reweighted Least Squares).

Iterative Reweighted Least Squares

1. Initialize ${\mathbf{\beta }}_0$, often obtained from Ordinary Least Squares. ${\hat{\mathbf{\beta }}}_0=({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$
2. Calculate the weight $w_{i,j}=\frac{\psi [(y_i-{\mathbf{x}}_i^{⊺}{\hat{\mathbf{\beta }}}_j)/s]}{(y_i-{\mathbf{x}}_i^{⊺}{\hat{\mathbf{\beta }}}_j)/s}$ If $y_i-{\mathbf{x}}_i^{⊺}\mathbf{\beta }=0$, then $w_i=1$
3. Update $\mathbf{\beta }$ using Weighted Least Squares ${\hat{\mathbf{\beta }}}_{j+1}=({\mathbf{X}}^{⊺}{\mathbf{W}}_j\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}{\mathbf{W}}_j\mathbf{y}$ where ${\mathbf{W}}_j=\mathrm{diag}(w_{1,j},w_{2,j},\dots ,w_{n,j})$
4. Iterate step 1 to 3 until the estimator converges.

Link to original

Inverse Regression

Inverse Regression

Definition

The problem predicting $x_0$ for a given response $y_0$ is called inverse regression (calibration, discrimination).

Consider a simple linear regression model $Y={\beta }_0+{\beta }_{1}X+ϵ$ The point estimation of $x_0$ for a given $y_0$ is obtained by ${\hat{x}}_0=\frac{y_0-{\hat{\beta }}_0}{{\hat{\beta }}_1}$ Also, the $100(1-\alpha )\%$ confidence interval for $x_0$ is obtained by the second-order inequality $d^2\left[{\hat{\beta }}_1^2-\frac{t^2{s}^2}{{\sum }_{i=1}^n(x_1-\bar{x}{)}^2}\right]-2d{\hat{\beta }}_1(y_0-\bar{y})+\left[(y_0-\bar{y}{)}^2-t^2{s}^2\left(1+\frac{1}{n}\right)\right]\le 0$ where $t=t_{\alpha /2}(n-2)$ and $s^2=\frac{{\sum }_{i=1}^n(Y_i-\hat{Y}{)}^2}{n-2}$. Let $d_1,d_2$ be solutions for the inequality, and we have a confidence interval $\bar{X}+d_1\le X_0\le \bar{X}+d_2$

Link to original

Biased Estimation

James-Stein Shrinkage Method

James-Stein Estimator

Definition

The James–Stein estimator is a biased estimator of the mean of correlated Multivariate Normal Distribution.

Let $\mathbf{z}\sim N_p(\mathbf{\mu },\mathbf{I})$. The James–Stein estimator is defined as ${\hat{\mathbf{\mu }}}_{JS}=\left(1-\frac{p-2}{{\mathbf{z}}^{⊺}\mathbf{z}}\right)\mathbf{z}$

The Least Squares Estimator ${\hat{\mathbf{\mu }}}_{LS}=\mathbf{z}$ is MLE and UMVUE. However, in terms of Mean Squared Error, James-Stein estimator ${\hat{\mathbf{\mu }}}_{JS}$ is better than the least squares estimator ${\hat{\mathbf{\mu }}}_{LS}$ for $p\ge 3$ cases.

Link to original

Ridge Regression

Ridge Regression

Definition

$\begin{array}{rl}{\hat{\beta }}_{\text{ridge}}& =({\mathbf{X}}^{⊺}\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}\\ & =\underset{\beta }{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\beta {)}^{⊺}(\mathbf{y}-\mathbf{X}\beta )+\lambda {\beta }^{⊺}\beta \\ & =\underset{\beta }{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\beta {)}^{⊺}(\mathbf{y}-\mathbf{X}\beta ),\text{subject to}{\beta }^{⊺}\beta \le c\end{array}$ where $\lambda \ge 0$ is a complexity parameter that controls the amount of shrinkage.

Ridge regression is particularly useful to mitigate the problem of Multicollinearity in linear regression

Facts

${\hat{\mathbf{y}}}_{\text{ridge}}=\mathbf{X}{\hat{\beta }}_{\text{ridge}}=\mathbf{X}({\mathbf{X}}^{⊺}\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}=\sum _{i=1}^p{\mathbf{u}}_j\frac{d_j^2}{d_j^2+\lambda }{\mathbf{u}}_j^{⊺}\mathbf{y}$ where $\mathbf{X}=\mathbf{UD}{\mathbf{V}}^{⊺}$ by Singular Value Decomposition

Link to original

Principal Component Regression

Principal Component Analysis

Definition

PCA is a linear dimensionality reduction technique. The correlated variables are linearly transformed onto a new coordinate system such that the directions capturing the largest variance in the data.

Population Version

Given a random vector $\mathbf{x}$, we find a $\mathbf{\alpha }$ such that $\mathrm{Var}({\mathbf{\alpha }}^{⊺}\mathbf{x})$ is maximized: $\underset{\mathbf{\alpha }}{\mathrm{argmax}}\mathrm{Var}({\mathbf{\alpha }}^{⊺}\mathbf{x})\text{s.t.}{\mathbf{\alpha }}^{⊺}\mathbf{\alpha }=1$ Equivalently, by the Method of Lagrange Multipliers with ${\mathbf{\alpha }}^{⊺}\mathbf{\alpha }=1$, $\underset{\mathbf{\alpha }}{\mathrm{argmax}}{\mathbf{\alpha }}^{⊺}\Sigma \mathbf{\alpha }-\lambda ({\mathbf{\alpha }}^{⊺}\mathbf{\alpha }-1)$ By differentiation, the $\mathbf{\alpha }$ is given by the eigen value problem $\Sigma \mathbf{\alpha }=\lambda \mathbf{\alpha }$ Thus the $\mathbf{\alpha }$ maximizing the variance of ${\mathbf{\alpha }}^{⊺}\mathbf{x}$ is the eigenvector corresponding to the largest Eigenvalue.

Sample Version

Given a data matrix $\mathbf{X}$, by Singular Value Decomposition, A matrix $\mathbf{X}$ can be factorized as $\mathbf{X}=\mathbf{UD}{\mathbf{V}}^{⊺}$. By algebra, $\mathbf{XV}=\mathbf{UD}=:\mathbf{Z}\Rightarrow \mathbf{X}{\mathbf{v}}_i=d_i{\mathbf{u}}_i=:{\mathbf{z}}_i$, where we call ${\mathbf{z}}_i$ the $i$-th principal component.

Facts

Since $\mathrm{Var}({\mathbf{z}}_i)=\mathrm{Var}(\mathbf{X}{\mathbf{v}}_i)=\frac{d_i^2}{n}$ and $d_1\ge d_2\ge \cdots \ge d_p\ge 0$ $\mathrm{Var}({\mathbf{z}}_1)\ge \mathrm{Var}({\mathbf{z}}_2)\ge \cdots \ge \mathrm{Var}({\mathbf{z}}_p)\ge 0$

Link to original

Principal Component Regression

Definition

${\hat{\mathbf{y}}}_{\text{PCR}}=\hat{\beta }\mathbf{Z}$ where ${\mathbf{Z}}_{n\times M}(M\le p)$ is a matrix whose columns are ${\mathbf{z}}_i:=\mathbf{X}{\mathbf{v}}_i=d_i{\mathbf{u}}_i$. In PCR, instead of regressing the dependent variables on explanatory variables directly, principal components(PC) of the explanatory variables are used as regressors. One typically uses only first a few PCs for regression, making PCR a kind of regularized procedure.

Facts

Since ${\mathbf{z}}_m$ are linear combination of the original ${\mathbf{x}}_j$, we can express the solution in therms of coefficients of the ${\mathbf{x}}_j$.

Link to original

PLS: Partial Least Squares

Partial Least Squares Regression

Definition

${\hat{\phi }}_m=\underset{\alpha }{\mathrm{argmax}}{\mathrm{Corr}}^2(\mathbf{y},\mathbf{X}a)\mathrm{Var}(\mathbf{X}\alpha )\text{subject to}||\alpha ||=1,{\alpha }^{⊺}\mathbf{S}{\hat{\phi }}_l=0,l=1,\dots ,m-1$ where $\mathbf{S}=\frac{1}{n}{\mathbf{X}}^{⊺}\mathbf{X}$ is the sample covariance matrix.

Unlike PCR maximizing $\mathrm{Var}(\mathbf{X}a)$ only, PLS finds directions maximizing both $\mathrm{Var}(\mathbf{X}a)$ and ${\mathrm{Corr}}^2(\mathbf{y},\mathbf{X}a)$

Link to original

LASSO

Lasso Regression

Definition

$\begin{array}{rl}{\hat{\beta }}_{\text{lasso}}& =\underset{\beta }{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\beta {)}^{⊺}(\mathbf{y}-\mathbf{X}\beta )+\lambda ||\beta ||\\ & =\underset{\beta }{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\beta {)}^{⊺}(\mathbf{y}-\mathbf{X}\beta ),\text{subject to}||\beta ||\le c\end{array}$

Lasso model assume that the coefficients of the model are sparse.

Link to original

Bayes Estimator and Biased Estimators

Bayes Theorem

Definition

Discrete Case

$P(C_j|C)=\frac{P(C|C_j)P(C_j)}{P(C)}=\frac{P(C|C_j)P(C_j)}{\sum _{i=1}^{k}P(C|C_i)P(C_i)},j=1,\dots ,k$ where $\sum _{i=1}^{k}P(C_i)P(C|C_i)=P(C)$ by Law of Total Probability

$P(C_j)$ is called a prior probability, and $P(C_j|C)$ is called a posterior probability

Continuous Case

$f(\theta |x;\gamma )=\frac{f(x|\theta )f(\theta ;\gamma )}{f(x)}=\frac{f(x|\theta )f(\theta ;\gamma )}{\int f(x|\theta )f(\theta ;\gamma )d\theta }\propto f(x|\theta )f(\theta ;\gamma )$ where $\theta$ is not a constant, but an unknown parameter follows a certain distribution with a parameter $\gamma$.

$f(\theta ;\gamma )$ is called a prior probability, $f(x|\theta )$ is called a likelihood, $f(x)$ is called an evidence or marginal likelihood, and $f(\theta |x;\gamma )$ is called a posterior probability

Parameter-Centric Notation

$p(\theta |x)=\frac{p(x|\theta )\pi (\theta )}{p(x)}$

Examples

Consider a random variable follows Binomial Distribution $X\sim B(n,\theta )$ and a prior distribution follows Beta Distribution $\pi (\theta ;\gamma )\sim \mathrm{Beta}(\alpha ,\beta )$ where $\gamma =(\alpha ,\beta )$.

The PDFs are defined as $p(x|\theta )=(\genfrac{}{}{0ex}{}{n}{x}){\theta }^x(1-\theta {)}^{n-x}$ and $\pi (\theta ;\gamma )=\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1}$. Then, by Bayes theorem, $p(\theta |x)\propto p(x|\theta )\pi (\theta ;\gamma )\propto {\theta }^{x+\alpha -1}(1-{\theta }^{n-x+\beta -1})\sim \mathrm{Beta}(x+\alpha ,n-x+\beta )$ Under Squared Error Loss, the Bayes Estimator is a mean of the posterior distribution. ${\hat{\delta }}_{\text{Bayes}}=\frac{x+\alpha }{n+\alpha +\beta }$

Link to original

Loss Function

Definition

Let $\theta \in \Omega$ be a parameter, $X:=u(X_1,X_2,\dots ,X_n)$ be a Statistic for the parameter $\theta$, and $\delta (x)$ be a Decision Function

A loss function is a non-negative function defined as $L(\theta ,\delta (x))$

It indicates the difference or discrepancy between $\theta$ and $\delta (x)$

Examples

• Absolute Error Loss
• Squared Error Loss
• Sum of Squared Errors Loss
• Cross-Entropy Loss
• Goal Post Error Loss
• Huber Loss
• Binary Loss
• Triplet Loss
• Pairwise Loss

Link to original

Risk Function

Definition

$R(\theta ,\delta )=E[L(\theta ,\delta (x))]={\int }_{X}L(\theta ,\delta (x))f(x|\theta )dx$

Risk function is an expectation of Loss Function

Link to original

Bayes Risk

Definition

r(\theta, \delta) &= \int_{\Theta} R(\theta, \delta) \pi(\theta) d\theta = E_\theta[R(\theta, \delta)] = E_\theta[E_{x}[L(\theta, \delta(x))]]\\ &= \int_{\Theta} \int_{X} L(\theta, \delta(x))p(x|\theta) \pi(\theta) dx d\theta = \int_{X} \int_{\Theta} L(\theta, \delta(x))p(\theta|x)p(x) d\theta dx\\ &= \int_{X}E_\theta[L(\theta, \delta(x))|X=x]p(x)dx = \int_{X}\rho(x, \pi)p(x)dx \end{aligned}$$ where $L(\theta, \delta)$ is [[Loss Function]], $R(\theta, \delta)$ is [[Risk Function]], and $\rho(x, \pi)$ is a [[Posterior Risk]].Link to original

Bayes Estimator

Definition

${\hat{\delta }}_{\text{Bayes}}=\underset{\delta }{\mathrm{argmin}}r(\theta ,\delta )=\underset{\delta }{\mathrm{argmin}}\rho (\pi ,\delta )$

Estimator that minimizes the Bayes Risk $r(\theta ,\delta )$ or Posterior Risk $\rho (\pi ,\delta )$

Facts

Under Squared Error Loss, the Bayes estimator is a posterior mean, and a posterior mode under Absolute Error Loss.

Consider a regression model $\mathbf{y}\sim N_n(\mathbf{X}\mathbf{\beta },{\sigma }^2\mathbf{I})$ and the prior distribution for $\mathbf{\beta }$, $\mathbf{\beta }\sim N_p(\mathbf{m},{\sigma }^2\mathbf{V})$. Then, the Bayes estimator under the Squared Error Loss is obtained as ${\hat{\mathbf{\beta }}}_{\text{Bayes}}=({\mathbf{X}}^{⊺}\mathbf{X}+{\mathbf{V}}^{-1}{)}^{-1}({\mathbf{V}}^{-1}\mathbf{m}+{\mathbf{X}}^{⊺}\mathbf{y})$ If $\mathbf{V}={\lambda }^{-1}{\mathbf{I}}_p$ for some $\lambda >0$ and $\mathbf{m}=0$, then the Bayes estimator is the same as ridge estimator ${\hat{\mathbf{\beta }}}_{\text{Bayes}}={\hat{\mathbf{\beta }}}_{\text{ridge}}=({\mathbf{X}}^{⊺}\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{y}$ If $\mathbf{V}={\lambda }^{-1}({\mathbf{X}}^{⊺}\mathbf{X}{)}^{-1}$ and $\mathbf{m}=0$, then the Bayes estimator is the James-Stein regression estimator ${\hat{\mathbf{\beta }}}_{\text{Bayes}}=(1+\lambda {)}^{-1}\hat{\mathbf{\beta }}$

Link to original

Empirical Bayes Estimator

Definition

Empirical Bayes estimator is a Bayes Estimator whose prior distribution is estimated from the data.

The parameter of the prior distribution $\hat{\gamma }$ is estimated by Maximum Likelihood Estimation. $\hat{\gamma }=\underset{\gamma }{\mathrm{argmax}}p(x|\gamma )=\underset{\gamma }{\mathrm{argmax}}\int p(x|\theta )\pi (\theta ;\gamma )d\theta$ And a posterior distribution is calculated with the $\hat{\gamma }$. $p(\theta |x;\hat{\gamma })\propto p(x|\theta )\pi (\theta ;\hat{\gamma })$ Estimator is obtained using the posterior distribution.

Examples

Consider data following a Poission Distribution $X_i\sim \mathrm{Pois}({\lambda }_i)$ and a prior distribution follows a Gamma Distribution $\pi (\lambda ;\gamma )\sim \Gamma (\alpha ,\beta )$ where $\gamma =(\alpha ,\beta )$ with $\alpha$ known, $\beta$ unknown.

Then the marginal likelihood is defined as

p(x_{i}|\beta) &= \int \operatorname{Pois}(x|\lambda_{i})\Gamma(\lambda_{i};\alpha, \beta) d\lambda_{i} = \int \left[ \frac{e^{-\lambda_{i}}\lambda^{x_{i}}}{x_{i}!} \right]\left[ \frac{\beta^{\alpha}\lambda_{i}^{\alpha-1}e^{-\beta \lambda_{i}}}{\Gamma(\alpha)} \right]d\lambda_{i}\\ &= \binom{x_{i}+\alpha-1}{\alpha-1}\left( \frac{\beta}{\beta+1} \right)^{\alpha} \left( \frac{1}{\beta+1} \right)^{x_{i}} \sim NB\left( \alpha, \frac{\beta}{\beta+1} \right) \end{aligned}$$ And the [[Maximum Likelihood Estimation|MLE]] of $\beta$, $\hat{\beta}_\text{MLE}$ is $$\hat{\beta}_\text{MLE} = \underset{\beta}{\operatorname{argmax}} \prod_{i=1}^{n}p(x_{i}|\beta) = \frac{\alpha}{\bar{X}}$$ The posterior distribution with $\hat{\beta}_\text{MLE}$ is defined as $$p(\lambda_{i}|x; \hat{\beta}_{\text{MLE}}) \propto p(x|\lambda_{i})\pi(\lambda_{i};\alpha, \hat{\beta}) \sim \Gamma(x_{i}+\alpha, 1+\hat{\beta})$$ Under [[Squared Error Loss]], the [[Bayes Estimator]] is a mean of the posterior distribution $\Gamma(x_{i}+\alpha, 1+\hat{\beta})$. $$\hat{\delta}_{\text{Bayes}} = \frac{\bar{X}(X_{i}+\alpha)}{\bar{X} + \alpha}$$ # Facts > Assume that $\mathbf{z} \sim N_{p}(\boldsymbol{\mu}, \mathbf{I})$ and the prior distribution for $\boldsymbol{\mu}$ is $\boldsymbol{\mu} \sim N_{p}(\mathbf{0}, \sigma^{2}\mathbf{I})$. > Then, the empirical Bayes estimator under the [[Squared Error Loss]] is [[James-Stein Estimator]] $\hat{\boldsymbol{\mu}}_{JS} = \left( 1 - \cfrac{p-2}{\mathbf{z}^{\intercal}\mathbf{z}} \right)\mathbf{z}$Link to original

Generalized Linear Model

Exponential Family

Exponential Family

Definition

General Representation

$p(x|\mathbf{\theta })-\frac{1}{Z(\mathbf{\theta })}h(\mathbf{x})\exp (\mathbf{t}(\mathbf{x})\cdot \mathbf{\theta })$ Where:

• $\mathbf{\theta }$ is the vector of natural parameters.
• $\mathbf{t}(\mathbf{x})$ is the vector of sufficient statistics.
• $h(\mathbf{x})$ is the non-negative volume of $x$
• $Z(\mathbf{\theta })=\int h(\mathbf{x})\exp (\mathbf{t}(\mathbf{x})\cdot \mathbf{\theta })\nu (d\mathbf{x})$ is the normalizer $\nu (d\mathbf{x})$ refers to the measure of $\mathbf{x}$.

Canonical Form with Dispersion Parameter

$f(x|\theta ,\varphi )=\exp \left[\frac{x\theta -b(\theta )}{a(\varphi )}+c(x,\varphi )\right]$ where $a(\cdot )$, $b(\cdot )$, and $c(\cdot )$ are known functions, $V(\mu ):=b^{\prime \prime }(\theta )$ is called a variance function, and $\varphi$ is called a dispersion parameter. This form is mainly used for GLMs.

Examples

Distributions	Normal	Poisson	Binomial	Gamma	Inverse Gaussian
Notation	$N(\mu ,{\sigma }^2)$	$\mathrm{Pois}(\lambda )$	$B(n,\pi )$	$\Gamma (\alpha ,\beta )$	$IG(\mu ,\lambda )$
$a(\varphi )$	$\varphi$	$1$	$1/n$	$\varphi$	$\varphi$
$b(\theta )$	${\theta }^2/2$	$e^{\theta }$	$\ln (1+e^{\theta })$	$-\ln (-\theta )$	$-(-2\theta {)}^{1/2}$
$c(x,\varphi )$	$-\frac{1}{2}[\ln (2\pi {\sigma }^2)+\frac{x^2}{{\sigma }^2}]$	$-\ln (x!)$	$\ln \left(\genfrac{}{}{0ex}{}{n}{x}\right)$	$\alpha \ln (\alpha x)-\ln \Gamma (\alpha )-\alpha \ln (x)$	$-\frac{1}{2}[\ln (2\pi \lambda x^3)+\frac{1}{\lambda x}]$
$\mu (\theta )=b^{\prime }(\theta )$	$\theta$	$e^{\theta }$	$e^{\theta }/(1+e^{\theta })$	$-1/\theta$	$(-2\theta {)}^{-1/2}$
$V(\mu )$	$1$	$\mu$	$\mu (1-\mu )$	${\mu }^2$	${\mu }^3$
$\varphi$	${\sigma }^2$	$1$	$1/n$	$1/\alpha$	$\lambda$
Natural link	Identity	log	logit	inverse	${\mu }^{-2}$

Facts

For an exponential family, $E(X)=\mu =b^{\prime }(\theta )$ and $\mathrm{Var}(X)=b^{\prime \prime }(\theta )a(\varphi )$ holds by Bartlett Identities

Link to original

Bartlett Identities

Definition

First Bartlett Identity

$0=E\left[\frac{\partial \ln f(X|\theta )}{\partial \theta }\right]=E[s(\theta |x)]$ where $f(X|\theta )$ is a Likelihood Function and $s(\theta |x)$ is a Score Function

Second Bartlett Identity

$0=E\left[\frac{{\partial }^2\ln f(X|\theta )}{\partial {\theta }^2}\right]+E\left[{\left(\frac{\partial \ln f(X|\theta )}{\partial \theta }\right)}^2\right]$ where $f(X|\theta )$ is a Likelihood Function

Link to original

Score Function

Definition

$s(\theta |x):=\frac{\partial \ln L(\theta |x)}{\partial \theta }$

The gradient of the log-likelihood function with respect to the parameter vector. The score indicates the steepness of the log-likelihood function

Facts

The score will vanish at a local Extremum

Link to original

Fisher Information

Definition

Fisher Information

I(\theta) &:= E\left[ \left( \frac{\partial \ln f(X|\theta)}{\partial \theta} \right)^{2} \right] = \int_{\mathbb{R}} \left( \frac{\partial \ln f(X|\theta)}{\partial \theta} \right)^{2} p(x, \theta)dx\\ &= -E\left[ \frac{\partial^{2} \ln f(X|\theta)}{\partial \theta^{2}} \right] = -\int_{\mathbb{R}} \frac{\partial^{2} \ln f(X|\theta)}{\partial \theta^{2}} p(x, \theta)dx\\ \end{aligned}$$ by the [[Bartlett Identities#second-bartlett-identity|second Bartlett identity]] $$I(\theta) = \operatorname{Var}\left( \frac{\partial \ln f(X|\theta)}{\partial \theta} \right) = \operatorname{Var}(s(\theta|x))$$ where $s(\theta|x)$ is a [[Score Function]] ## Fisher Information Matrix ![[Pasted image 20231224171415.png|800]] Let $\mathbf{X}$ be a [[Random Vector]] with [[Density Function|PDF]] $f(x|\boldsymbol{\theta})$, where $\boldsymbol{\theta} \in \Omega \subset R^{p}$, then the **Fisher information matrix** for on $\boldsymbol{\theta}$ is a $p \times p$ matrix defined as $$I(\boldsymbol{\theta}) := \operatorname{Cov}\left( \cfrac{\partial}{\partial \boldsymbol{\theta}} \ln f(x|\boldsymbol{\theta}) \right) = E\left[ \left( \cfrac{\partial}{\partial \boldsymbol{\theta}} \ln f(x|\boldsymbol{\theta}) \right) \left( \cfrac{\partial}{\partial \boldsymbol{\theta}} \ln f(x|\boldsymbol{\theta}) \right)^\intercal \right] = -E\left[ \cfrac{\partial^{2}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^\intercal} \ln f(x|\boldsymbol{\theta}) \right]$$ and $jk$-th element of $I(\boldsymbol{\theta})$, $I_{jk} = - E\left[ \cfrac{\partial^{2}}{\partial \theta_{j} \partial\theta_{k}} \ln f(x|\boldsymbol{\theta}) \right]$ # Properties ## Chain Rule The information in length $n$ [[Random Sample]] $X_{1}, X_{2}, \dots, X_{n}$ is $n$ times the information in a single sample $X_{i}$ $I_\mathbf{X}(\theta) = n I_{X_{1}}(\theta)$ # Facts > In a location model, information is not dependent on a location parameter. > $$I(\theta) = \int_{-\infty}^{\infty}\left( \frac{f'(z)}{f(z)} \right)^{2} f(z)dz$$Link to original

Construction of GLMs

Generalized Linear Model

Definition

A generalized linear model (GLM) is a generalization of linear regression model.

The GLM consists of three elements.

1. Response variable: The i.i.d. response variables $Y_i$ follow Exponential Family.
2. Linear predictor: $\mathbf{\eta }:={\mathbf{X}}^{⊺}\mathbf{\beta }$ is called a linear predictor
3. Link function: There exists a link function satisfies $g(\mu )=\mathbf{\eta }$, where $g$ is assumed to be monotone and differentiable.

Link Functions

Link functions for binomial response.

• Logit: $\eta =\ln \left(\frac{\mu }{1-\mu }\right)$
• Probit: $\eta ={\Phi }^{-1}(\mu )$ where $\Phi$ is a CDF of Normal Distribution
• Complementary log-log: $\eta =\ln [-\ln (1-\mu )]$

For an Exponential Family, the link function satisfying $\theta =\eta$ is called natural link.

Facts

The parameter $\mathbf{\beta }$ is estimated by MLE. Since the score function is non-linear, it can’t have explicit solution, so we use Newton–Raphson method or Fisher’s Scoring Method.

Link to original

Estimation of Regression Coefficients

Newton's Method

Definition

An iterative algorithm for finding the roots of a differentiable function, which are solution to the equation $f(x)=0$

Algorithm

Find the next point such that the Taylor series of the given point is 0 Taylor first approximation: $f(x)\simeq f(x_n)+f^{\prime }(x_n)(x-x_n)$ The point such that the Taylor series is 0: $x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$

multivariate version: ${\mathbf{x}}_{n+1}=f({\mathbf{x}}_n)-\nabla f({\mathbf{x}}_n{)}^{-1}f({\mathbf{x}}_n)$

In convex optimization,

Find the minimum point^[its derivative is 0] of Taylor quadratic approximation. Taylor quadratic approximation: $f(x)\simeq f(x_n)+f^{\prime }(x_n)(x-x_n)+\frac{f^{\prime \prime }(x_n)}{2}(x-x_n{)}^2$ The derivative of the quadratic approximation: $f^{\prime }(x_n)+f^{\prime \prime }(x_n)(x-x_n)$ The minimum point of the quadratic approximation^[the point such that the derivative of the quadratic approximation is 0]: $x_{n+1}=x_n-\frac{f^{\prime }(x_n)}{f^{\prime \prime }(x_n)}$ multivariate version: ${\mathbf{x}}_{n+1}={\mathbf{x}}_n-{\nabla }^{2}f({\mathbf{x}}_n{)}^{-1}\nabla f({\mathbf{x}}_n)$

Examples

Solution of a linear system

Solve $A\mathbf{x}=\mathbf{b}$ with an MSE loss

The cost function is $f(\mathbf{x})=(A\mathbf{x}-\mathbf{b}{)}^{⊺}(A\mathbf{x}-\mathbf{b})$ and its gradient and hessian are $\nabla f(x)=-2(\mathbf{b}-A\mathbf{x}{)}^{⊺}$, ${\nabla }^{2}f(x)=2A^{⊺}A$

Then, solution is ${\mathbf{x}}_{n+1}={\mathbf{x}}_n+(A^{⊺}A{)}^{-1}{A}^{⊺}(\mathbf{b}-A{\mathbf{x}}_n)$ If $(A^{⊺}A)$ is invertible, ${\mathbf{x}}_{n+1}=(A^{⊺}A{)}^{-1}{A}^{⊺}\mathbf{b}$ is a Least Square solution.

Link to original

Fisher's Scoring Method

Definition

Fisher’s scoring method is a variation of Newton–Raphson method that uses Fisher Information instead of the Hessian Matrix

${\hat{\theta }}_{n+1}={\hat{\theta }}_n-\frac{l^{\prime }({\hat{\theta }}_n)}{I({\hat{\theta }}_n)}$ where $I(\theta )=E[l^{\prime \prime }(\theta )]$ is the Fisher Information.

Link to original

Goodness-of-Fit Measures for GLMs

Deviance

Definition

Deviance is a goodness-of-fit statistic for Generalized Linear Model.

Let ${\hat{\beta }}_{\text{max}}$ be an estimator of $\beta$ under maximal model $(n=p)$. The deviance is defined as $D=2\ln \lambda =2[l{\hat{\beta }}_{\text{max}};y)-l\hat{\beta };y)]\stackrel{a}{\sim }{\chi }^2(n-p)$ where $l$ is a log-likelihood, $p$ is the number of parameter of the current model, and $n$ is the sample size.

Consider a null hypothesis $H_0$: the current model is not good, can be tested with the deviation. We reject $H_0$ if $D>{\chi }_{\alpha }^2(n-p)$

Examples

Distribution	Deviance
Normal	$\sum (y_i-{\hat{\mu }}_i{)}^2/{\sigma }^2$
Poisson	$2\sum [y_i\log (y_i/{\hat{\mu }}_i)-(y_i-{\hat{\mu }}_i)]$
Binomial	$2\sum [y_i\log (y_i/{\hat{\mu }}_i)+(m_i-y_i)\log [(m_i-y_i)/(m_i-{\hat{\mu }}_i)]]$
Gamma	$2\sum [-\log (y_i/{\hat{\mu }}_i)+(y_i-{\hat{\mu }}_i)/{\hat{\mu }}_i]$
Inverse Gaussian	$\sum (y_i-{\hat{\mu }}_i{)}^2/({\hat{\mu }}_i^2{y}_i)$

Link to original

Pearson's Chi-squared Statistic

Definition

The Pearson’s ${\chi }^2$ statistic is goodness-of-fit measure for GLM defined as $X^2=\sum _{i=1}^n\frac{(y_i-{\hat{\mu }}_i{)}^2}{V({\hat{\mu }}_i)}\stackrel{a}{\sim }{\chi }^2(n-p)$ where $\hat{\mu }=\mu (\theta )=b^{\prime }(\theta )$, $V(\hat{\mu })=b^{\prime \prime }(\theta )$ is the variance function, $p$ is the number of parameter of the current model, and $n$ is the sample size.

Facts

Under the Gaussian distribution, the Deviance and the Pearson’s ${\chi }^2$ statistic is the same and follow Chi-squared Distribution $D=X^2\sim {\chi }^2(n-p)$

Link to original

Testing and Residuals

Goodness-of-Fit Test with Deviance

Definition

Assume that a current model with $p$-parameters $\mathbf{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_q,\dots ,{\beta }_{p-1}{)}^{⊺}$ and consider a null hypothesis $H_0:{\beta }_q=\cdots ={\beta }_{p-1}=0$. Then, the test Statistic is defined as $\Delta D=D_0-D_1=2[l({\hat{\mathbf{\beta }}}_1;\mathbf{y})-l({\hat{\mathbf{\beta }}}_0;\mathbf{y})]\stackrel{a}{\sim }{\chi }^2(p-q)$ where $D$ are the Deviance of the generalized linear models.

We reject $H_0$ if $\Delta D>{\chi }_{\alpha }^2(p-q)$

Link to original

Pearson Residual

Definition

The Pearson residual is defined as $r_{P_i}=\frac{y_i-{\hat{\mu }}_i}{\sqrt{V({\hat{\mu }}_i)}}$ where $\hat{\mu }=\mu (\theta )=b^{\prime }(\theta )$, $V(\hat{\mu })=b^{\prime \prime }(\theta )$ is the variance function.

Facts

$X^2=\sum _{i=1}^n{r}_{P_i}^2$ where $X^2$ is the Pearson’s Chi-squared Statistic

Link to original

Anscombe Residual

Definition

Anscombe Transformation

$A(x)=\int V^{-1/3}(x)dx$ where $V$ is the variance function.

Anscombe Residual

The Anscombe residual is the Anscombe Transformed Pearson Residual. The transformation make the residual follows a normal distribution. $r_{A_i}=\frac{A(y_i)-A({\hat{\mu }}_i)}{\sqrt{V(A({\hat{\mu }}_i))}}=\frac{A(y_i)-A({\hat{\mu }}_i)}{A^{\prime }({\hat{\mu }}_i)\sqrt{V({\hat{\mu }}_i)}}$ where $\hat{\mu }=\mu (\theta )=b^{\prime }(\theta )$, $V(\hat{\mu })=b^{\prime \prime }(\theta )$ is the variance function.

Link to original

Deviance Residual

Definition

The deviance residual is defined as $r_{D_i}=\mathrm{sign}(y_i-{\hat{\mu }}_i)\sqrt{d_i}=\mathrm{sign}(y_i-{\hat{\mu }}_i)\sqrt{2[l{\hat{\beta }}_{\text{max}};y_i)-l\hat{\beta };y_i)]}$ where $\hat{\mu }=\mu (\theta )=b^{\prime }(\theta )$, $V(\hat{\mu })=b^{\prime \prime }(\theta )$ is the variance function, and the Deviance $D=\sum _{i=1}^n{d}_i$

Link to original

ANOVA Models

Constraints for Dummy Variable

Definition

When using Dummy Variable, the design matrix $X=[1,2,D_1,\dots ,D_c]$ where $c$ is the number of groups, is not a full rank. The problem can be solved by assigning some constraints.

Sum-to-Zero Constraint

Sum-to-zero constraint uses row-wise demeaned dummy variables. $\sum _{i=1}^c{D}_j=0$

With sum-to-zero constraint, each coefficient indicates discrepancy from overall mean, and the intercept term indicates overall mean.

Corner-Point Constraint

Corner-point constraint omit a base group from the model. $X=[1,D_2,D_3,\dots ,D_c]$

With corner-point constraint, each coefficient indicates discrepancy from base group, and the intercept term indicates mean of base group.

Link to original

Deviance Test for One-way ANOVA

The null hypothesis $H_0:D_1=D_2=\cdots =D_{c-1}=0$ i.e. there is no treatment effect, can be tested with the Deviance. If ${\sigma }^2$ is known, the test statistic is defined as $\Delta D=D_0-D_1=\frac{1}{{\sigma }^2}\left(\frac{1}{n_c}\sum _{j=1}^c{y}_{j.}^2-\frac{1}{N}{Y}_{..}^2\right)\sim {\chi }^2(c-1)$ And reject the $H_0$ if $\Delta D>{\chi }_{\alpha }^2(c-1)$

If ${\sigma }^2$ is unknown, the test statistic is defined as $F=\frac{D_0-D_1}{(N-1)-(N-c)}/\frac{D_1}{N-c}\sim F(c-1,N-c)$ And reject the $H_0$ if $F>F_{\alpha }(c-1,N-c)$

Link to original

Deviance Test for Two-way ANOVA

For a two-way ANOVA with factors $A$ and $B$, we have three null hypotheses to test:

• $H_{0A}:{\alpha }_1={\alpha }_2={\alpha }_{a-1}=0$ i.e. there is no treatment effect of factor $A$
• $H_{0B}:{\beta }_1={\beta }_2={\beta }_{b-1}=0$ i.e. there is no treatment effect of factor $B$
• $H_{0AB}:{\gamma }_{ij}=0,\forall i,j$ i.e. there is no interaction effect between factor $A$ and $B$ They can be tested with the Deviance.

If ${\sigma }^2$ is known, the test statistic for $H_{0A}$ is defined as $\Delta D_A=D_0-D_1=\frac{1}{{\sigma }^2}\left(\frac{1}{bn}\sum _{i=1}^a{y}_{i..}^2-\frac{1}{N}{Y}_{...}^2\right)\sim {\chi }^2(a-1)$ And reject the $H_{0A}$ if $\Delta D_A>{\chi }^2(a-1)$

If ${\sigma }^2$ is unknown, the test statistic for $H_{0A}$ is defined as $F_A=\frac{D_0-D_1}{a-1}/\frac{D_1}{N-ab}\sim F(a-1,N-ab)$ And reject the $H_{0A}$ if $F_A>F(a-1,N-ab)$

If ${\sigma }^2$ is known, the test statistic for $H_{0B}$ is defined as $\Delta D_B=D_0-D_1=\frac{1}{{\sigma }^2}\left(\frac{1}{an}\sum _{j=1}^b{y}_{.j.}^2-\frac{1}{N}{Y}_{...}^2\right)\sim {\chi }^2(b-1)$ And reject the $H_{0B}$ if $\Delta D_B>{\chi }^2(b-1)$

If ${\sigma }^2$ is unknown, the test statistic for $H_{0B}$ is defined as $F_B=\frac{D_0-D_1}{b-1}/\frac{D_1}{N-ab}\sim F(b-1,N-ab)$ And reject the $H_{0B}$ if $F_B>F(b-1,N-ab)$

If ${\sigma }^2$ is known, the test statistic for $H_{0AB}$ is defined as $\Delta D_{AB}=D_0-D_1=\frac{1}{{\sigma }^2}\left(\frac{1}{n}\sum _{i=1}^a\sum _{j=1}^b{y}_{ij.}^2-\frac{1}{bn}\sum _{i=1}^a{Y}_{i..}^2-\frac{1}{an}\sum _{j=1}^b{Y}_{.j.}^2+\frac{1}{N}{Y}_{...}^2\right)\sim {\chi }^2((a-1)(b-1))$ And reject the $H_{0AB}$ if $\Delta D_{AB}>{\chi }^2((a-1)(b-1))$

If ${\sigma }^2$ is unknown, the test statistic for $H_{0AB}$ is defined as $F_{AB}=\frac{D_0-D_1}{(a-1)(b-1)}/\frac{D_1}{N-ab}\sim F((a-1)(b-1),N-ab)$ And reject the $H_{0AB}$ if $F_{AB}>F((a-1)(b-1),N-ab)$

Link to original

Logistic Regression Model

Logistic Regression

Definition

Logistic regression is a Generalized Linear Model with Logit link function. Logistic regression models the log-odds of an event as a linear combination of independent variables. $\mathrm{logit}(\pi )=\ln \left(\frac{\pi }{1-\pi }\right)={\mathbf{x}}^{⊺}\mathbf{\beta }$ The probability is calculated as $\pi (\mathbf{x})=P(Y=1|X=\mathbf{x})=\frac{\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })}{1+\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })}=\frac{1}{1+\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta })}$

Link to original

Weighted Least Squares for Binomial Distribution

Definition

Let $Y_i\sim B(m_i,{\pi }_i),i=1,\dots ,n$ and random variables $p_i=\frac{y_i}{m_i}$ and consider a variance stabilizing transformation $\psi (\cdot )$. Then, by Taylor expansion and Delta Method $E[\psi (p_i)]\approx \psi ({\pi }_i)$ and $\mathrm{Var}(\psi (p_i))\approx {\psi }^{\prime }({\pi }_i{)}^2\frac{{\pi }_i(1-{\pi }_i)}{m_i}$ hold. Now, the Pearson $X^2$ is defined as $X^2=\sum _{i=1}^n\frac{[\psi (p_i)-\psi ({\pi }_i){]}^2}{{\psi }^{\prime }({\pi }_i{)}^2{\pi }_i(1-{\pi }_i)/m_i}$ By minimizing $X^2$ we can find the WLS estimator $\hat{\pi }$. It can be seen as a Weighted Least Squares with $y_i=\psi (p_i)$, ${\hat{y}}_i=\psi ({\pi }_i)$, and $w_i=[{\psi }^{\prime }({\pi }_i{)}^2{\pi }_i(1-{\pi }_i)/m_i{]}^{-1}$

Widely Used Variance Stabilizing Transformations

• Logit function: $\psi (x)=\ln \left(\frac{x}{1-x}\right)$
• The arcsin squared-root function: $\psi (x)=\arcsin (\sqrt{x})$
• Empirical logistic function: $\psi (x)=\ln \left(\frac{x+0.5}{1-x+0.5}\right)$

Link to original

Multinomial Logistic Regression

Definition

Multinomial logistic regression is used for polychotomous (multi-class) data.

Suppose the number of classes is $K$ $\ln \frac{P(Y=k|X=x)}{P(Y=K|X=x)}={\mathbf{\beta }}_k^{⊺}x,k=1,\dots ,K-1$ The probability is calculated as ${\pi }_k=P(Y=k|X=\mathbf{x})=\{\begin{array}{ll}\frac{\exp ({\mathbf{\beta }}_k^{⊺}\mathbf{x})}{1+{\sum }_{l=1}^{K-1}\exp ({\mathbf{\beta }}_l^{⊺}\mathbf{x})}& k=1,\dots ,K-1\\ \frac{1}{1+{\sum }_{l=1}^{K-1}\exp ({\mathbf{\beta }}_l^{⊺}\mathbf{x})}& k=K\end{array}$ $\underset{1\le k\le K}{\mathrm{argmax}}{\hat{f}}_k(\mathbf{x})=\underset{1\le k\le K}{\mathrm{argmax}}\hat{P}(Y=k,X=\mathbf{x})$

Link to original

Proportional Odds Model

Definition

The proportional odds model is a Generalized Linear Model used for modeling the ordinal response.

Consider the response $Y$ having $k$ possible categories. $\mathrm{logit}({\gamma }_j)=\ln \left(\frac{{\gamma }_j}{1-{\gamma }_j}\right)={\alpha }_j-{\mathbf{\beta }}^{⊺}\mathbf{x},j=1,\dots ,k-1$ And the cumulative probability of a response is ${\gamma }_j(\mathbf{x})=P(Y\le j|\mathbf{x})$.

Link to original

Log-Linear Model

Contingency Table

Definition

$J\times K$ contingency table

	$B_1$	$B_2$	$\dots$	$B_K$	Total
$A_1$	$Y_{11}$	$Y_{12}$	$\dots$	$Y_{1K}$	$Y_{1.}$
$A_2$	$Y_{21}$	$Y_{22}$	$\dots$	$Y_{2K}$	$Y_{2.}$
$⋮$	$⋮$	$⋮$	$\ddots$	$⋮$	$⋮$
$A_J$	$Y_{J1}$	$Y_{J2}$	$\dots$	$Y_{JK}$	$Y_{J.}$
Total	$Y_{.1}$	$Y_{.2}$	$\dots$	$Y_{.K}$	$n=Y_{..}$

Contingency table is a table displays the multivariate frequency distribution of the variables.

Examples

Poisson and Multinomial Distribution Cases

2-Dimensional Case

Consider a $J\times K$ contingency table.

No Margin Fixed

If $Y_{jk}$‘s are independent and $Y_{jk}\sim \mathrm{Pois}({\lambda }_{jk})$, then the joint PDF becomes $f(\mathbf{y};\mathbf{\lambda })={\prod }_{j=1}^J{\prod }_{k=1}^K\frac{{\lambda }_{jk}^{y_{jk}}{e}^{-{\lambda }_{jk}}}{y_{jk}!}$

Fixed Total

Given the total $n=Y_{..}={\sum }_{j=1}^J{\sum }_{k=1}^K{Y}_{jk}$, the conditional distribution is Multinomial Distribution. $f(\mathbf{y}|Y_{..}=n)=n!{\prod }_{j=1}^J{\prod }_{k=1}^K\frac{{\theta }_{jk}^{y_{jk}}}{y_{jk}!}\sim \mathrm{Mult}(n,{\theta }_{11},{\theta }_{12},\dots ,{\theta }_{JK})$ where ${\theta }_{jk}^{y_{jk}}=\frac{{\lambda }_{jk}}{{\sum }_{j=1}^J{\sum }_{k=1}^K{\lambda }_{jk}}$ and ${\sum }_{j=1}^J{\sum }_{k=1}^K{\theta }_{jk}=1$

One Margin Fixed

Given one fixed margin (in this case, row) $n_j=Y_{j.}$, the conditional distribution of frequency of $j$-th row is Multinomial Distribution. $f(y_{j1},y_{j2},\dots ,y_{jK}|y_{j.})=y_{j.}!{\prod }_{k=1}^K\frac{{\theta }_{jk}^{y_{jk}}}{y_{jk}!}\sim \mathrm{Mult}(n,{\theta }_{j1},{\theta }_{j2},\dots ,{\theta }_{jK})$ And if each row is independent, the conditional distribution is product-Multinomial Distribution (joint multinomial distribution). $f(\mathbf{y}|y_{j.}=n_j,j=1,\dots ,J)={\prod }_{j=1}^J\left[y_{j.}!{\prod }_{k=1}^K\frac{{\theta }_{jk}^{y_{jk}}}{y_{jk}!}\right]$

3-Dimensional Case

Consider a $J\times K\times L$ contingency table.

No Margin Fixed

If $Y_{jkl}$‘s are independent and $Y_{jkl}\sim \mathrm{Pois}({\lambda }_{jkl})$, then the joint PDF becomes $f(\mathbf{y};\mathbf{\lambda })={\prod }_{j=1}^J{\prod }_{k=1}^K{\prod }_{l=1}^L\frac{{\lambda }_{jkl}^{y_{jkl}}{e}^{-{\lambda }_{jkl}}}{y_{jkl}!}$

Fixed Total

Given the total $n=Y_{...}={\sum }_{j=1}^J{\sum }_{k=1}^K{\sum }_{k=1}^L{Y}_{jkl}$, the conditional distribution is Multinomial Distribution. $f(\mathbf{y}|Y_{...}=n)=n!{\prod }_{j=1}^J{\prod }_{k=1}^K{\prod }_{l=1}^L\frac{{\theta }_{jkl}^{y_{jkl}}}{y_{jkl}!}\sim \mathrm{Mult}(n,{\theta }_{111},{\theta }_{112},\dots ,{\theta }_{JKL})$ where ${\theta }_{jkl}^{y_{jkl}}=\frac{{\lambda }_{jkl}}{{\sum }_{j=1}^J{\sum }_{k=1}^K{\sum }_{l=1}^L{\lambda }_{jkl}}$ and ${\sum }_{j=1}^J{\sum }_{k=1}^K{\sum }_{l=1}^L{\theta }_{jkl}=1$

One Margin Fixed

Given one fixed margin (in this case, $L$) $n_l=Y_{..l}$, and each $l$-th rows are independent, the conditional distribution is product-Multinomial Distribution (joint multinomial distribution). $f(\mathbf{y}|y_{..l}=n_l,l=1,\dots ,L)={\prod }_{l=1}^L\left[y_{..l}!{\prod }_{j=1}^J{\prod }_{k=1}^K\frac{{\theta }_{jkl}^{y_{jkl}}}{y_{jkl}!}\right]$

Two Margin Fixed

If two margins are fixed (in this case, $J$ and $L$) $n_{jl}=Y_{j.l}$, and each $j$-th and $l$-th rows are independent, the conditional distribution is product-Multinomial Distribution (joint multinomial distribution). $f(\mathbf{y}|y_{j.l}=n_{jl},j=1,\dots ,J;l=1,\dots ,L)={\prod }_{j=1}^J{\prod }_{l=1}^L\left[y_{j.l}!{\prod }_{k=1}^K\frac{{\theta }_{jkl}^{y_{jkl}}}{y_{jkl}!}\right]$

Link to original

Poisson Regression

Definition

Poisson Regression

Log-linear model is a Generalized Linear Model with log link function. Log-linear model models the log expected counts of an event as a linear combination of independent variables. $\ln (\mu )=\ln [E(y)]={\mathbf{x}}^{⊺}\mathbf{\beta }$

Log-Linear Model for Contingency Table

Consider a $J\times K$ Contingency Table, then the model is defined as ${\eta }_{jk}=\ln ({\mu }_{jk})=\ln [E(y_{jk})]=\mu +{\alpha }_j+{\beta }_k+(\alpha \beta {)}_{jk}$ where $\mu$ is the overall effect, ${\alpha }_j$ and ${\beta }_k$ are main effect, and $(\alpha \beta {)}_{jk}$ are interaction effect.

The Deviance of the model is defined as $D=2\left[\sum _{j=1}^J\sum _{k=1}^K\left(y_{jk}\ln \left(\frac{y_{jk}}{{\hat{\mu }}_{jk}}\right)-(y_{jk}-{\widehat{mu}}_{jk})\right)\right]$ and the Pearson’s Chi-squared Statistic is defined as $X^2=\sum _{j=1}^J\sum _{k=1}^K\frac{(y_{jk}-{\hat{\mu }}_{jk}{)}^2}{{\hat{\mu }}_{jk}}$

Link to original

Overdispersion

Overdispersion

Definition

Overdispersion is the presence of greater observed variance than what would be expected under a given statistical model. It can happen due to heterogeneity or lack of independence between trials, and clustering or grouping in the data.

Overdispersion on Binomial Distribution

Assume that there are $k$ clusters, the number of observations of each cluster is $m$. Let the Random Variable following Binomial Distribution $X_i\sim B(m,{\pi }_i)$ be the number of successes out of $m$ observation, where ${\pi }_i$ is also Random Variable with mean $E({\pi }_i)=\pi$ and variance $\mathrm{Var}({\pi }_i)={\sigma }^2\pi (1-\pi )$. Also, let $Y=\sum _{i=1}^k{X}_i\sim B(n,\pi )$ be the total number of successes in the all clusters, where $m:=mk$. Then, the mean and variance of $Y$ is defined as $E(Y)=n\pi$ $Var(Y)=n\pi (1-\pi )[1+(m-1){\sigma }^2]$ Here, $\varphi =1+(m-1){\sigma }^2$ is called a dispersion parameter. If $\varphi >1$, then overdispersion occurred.

The Lexis Ratio is also used for detecting overdispersion $Q=\frac{{\sum }_{i=1}^k({\hat{\pi }}_i-\hat{\pi }{)}^2/(k-1)}{\hat{\pi }(1-\hat{\pi })/n}$ where ${\hat{\pi }}_i=\frac{X_i}{m}$ and $\hat{\pi }=\frac{Y}{n}$

Overdispersion on Poisson Distribution

If $\mathrm{Var}(Y)>E(Y)$ in Poission Distribution setting, then the overdispersion might occur. Let i.i.d. $Z_i$ be the number of occurrences at the $i$-th cluster, and the number of clusters $N$ is also a Random Variable. Also, let $Y=\sum _{i=1}^N{Z}_i$ be the total number of occurrences in $N$ clusters, where $N$ follows a Poission Distribution and independent to $Z$. Then, the mean and variance of $Y$ is defined as $E(Y)=E(N)E(Z)$ $\mathrm{Var}(Y)=E(N)E(Z^2)$ If $E(Z^2)\gg E(Z)$, then the overdispersion occurs.

Link to original

Nonlinear Regression Model

Nonlinear Regression Model

Nonlinear Regression Model

Definition

Let $Y$ be a response and $\mathbf{X}=(X_1,X_2,\dots ,X_k{)}^{⊺}$ be $k$-dimensional explanatory variable. The non-linear regression model can be written as $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ where $f({\mathbf{X}}_i,\mathbf{\theta })$ is a known regression function and $\mathbf{\theta }=({\theta }_1,{\theta }_2,\dots ,{\theta }_p{)}^{⊺}$ is a $p$-dimensional parameter vector.

The parameter vector $\mathbf{\theta }$ can not be obtained analytically. So it estimated numerically with such methods Newton–Raphson method, Gauss-Newton Method

Facts

If the function $f$ is linear, then the model is Multiple Linear Regression.

Link to original

Estimation of Parameter

Gauss-Newton Method

Definition

Gauss-Newton method is an iterative algorithm used to solve used to solve non-linear least squares problems. This method approximates Hessian Matrix using the Jacobian Matrix.

Algorithm

Consider a non-linear regression model with $k$-dimensional explanatory variable and $p$-dimensional parameter vector $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ The first order Taylor expansion gives the linear approximation of the model. $\mathbf{y}=\mathbf{\eta }({\mathbf{\theta }}_0)+{\mathbf{J}}_0(\mathbf{\theta }-{\mathbf{\theta }}_0)+ϵ$ where ${\mathbf{\theta }}_0=({\theta }_{1,0},{\theta }_{2,0},\dots ,{\theta }_{p,0}{)}^{⊺}$ is the initial vector for $\mathbf{\theta }$ given by domain knowledge, $\mathbf{\eta }({\mathbf{\theta }}_0)=[f({\mathbf{X}}_1,{\mathbf{\theta }}_0),f({\mathbf{X}}_2,{\mathbf{\theta }}_0),\dots ,f({\mathbf{X}}_n,{\mathbf{\theta }}_0){]}^{⊺}$ is the estimation vector with ${\mathbf{\theta }}_0$, and ${\mathbf{J}}_0=\left(\frac{\partial f({\mathbf{X}}_i,\mathbf{\theta })}{\partial {\theta }_j}\right){|}_{\mathbf{\theta }={\mathbf{\theta }}_0}$ is the $n\times p$ Jacobian Matrix of the function $f$ at $\mathbf{\theta }={\mathbf{\theta }}_0$.

Now, the approximated model has a Multiple Linear Regression. Since $(\mathbf{y}-\mathbf{\eta }({\mathbf{\theta }}_0))$ and ${\mathbf{J}}_0$ are constants given ${\mathbf{\theta }}_0$, we can find LSE of $(\mathbf{\theta }-{\mathbf{\theta }}_0)$. $({\hat{\mathbf{\theta }}}_1-{\hat{\mathbf{\theta }}}_0)=({\mathbf{J}}_0^{⊺}{\mathbf{J}}_0{)}^{-1}{\mathbf{J}}_0^{⊺}(\mathbf{y}-\mathbf{\eta }({\theta }_0))$ And the updating formula is defined as ${\hat{\mathbf{\theta }}}_{n+1}={\hat{\mathbf{\theta }}}_n-({\mathbf{J}}_n^{⊺}{\mathbf{J}}_n{)}^{-1}{\mathbf{J}}_n^{⊺}(\mathbf{y}-\mathbf{\eta }({\hat{\mathbf{\theta }}}_n))$ The updating is repeated iteratively until it converges. ${\hat{\mathbf{\theta }}}_i$ converges to the $\mathbf{\theta }$ minimizing Sum of Squared Errors Loss with a proper initial value ${\mathbf{\theta }}_0$.

Link to original

Inference in the Nonlinear Model

Distribution of Parameter of Nonlinear Regression

Definition

Suppose a Nonlinear Regression Model $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ If the error term ${ϵ}_i$ follows a Normal Distribution, LSE $\hat{\mathbf{\theta }}$ is the same as MLE.

For a large sample size, the distribution of the estimated parameter $\hat{\mathbf{\theta }}$ follows a Normal Distribution $\hat{\mathbf{\theta }}\stackrel{D}{\to }N(\mathbf{\theta },({\hat{\mathbf{J}}}^{⊺}\hat{\mathbf{J}}{)}^{-1}{\sigma }^2)$ where $\hat{\mathbf{J}}=\left(\frac{\partial f({\mathbf{X}}_i,\mathbf{\theta })}{\partial {\theta }_j}\right){|}_{\mathbf{\theta }=\hat{\mathbf{\theta }}}$ is the Jacobian Matrix of $f$ and ${\sigma }^2$ is the variance of the error term.

The variance of the error term can be estimated by $s^2=\frac{{\sum }_{i=1}^n(Y_i-{\hat{Y}}_i{)}^2}{n-p}=\frac{{\sum }_{i=1}^n(Y_i-f({\mathbf{X}}_i,\hat{\mathbf{\theta }}){)}^2}{n-p}$

Link to original

Standard Error for Parameter of Nonlinear Regression

Definition

Suppose a Nonlinear Regression Model $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ The standard error for parameter of nonlinear regression ${\theta }_j$ is calculated as $SE({\hat{\theta }}_j)=\sqrt{c_j}s$ where $c_j=\mathrm{diag}[({\hat{\mathbf{J}}}^{⊺}\hat{\mathbf{J}}{)}^{-1}{]}_j$, $s^2=\frac{{\sum }_{i=1}^n(Y_i-f({\mathbf{X}}_i,\hat{\mathbf{\theta }}){)}^2}{n-p}$, and $\hat{\mathbf{J}}=\left(\frac{\partial f({\mathbf{X}}_i,\mathbf{\theta })}{\partial {\theta }_j}\right){|}_{\mathbf{\theta }=\hat{\mathbf{\theta }}}$is the Jacobian Matrix of $f$.

Link to original

Joint Confidence Region for Nonlinear Regression

Definition

Suppose a Nonlinear Regression Model $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ The $100(1-\alpha )\%$ join confidence region for $\mathbf{\theta }$ is obtained as $\{\mathbf{\theta }|(\mathbf{\theta }-\hat{\mathbf{\theta }}{)}^{⊺}{\hat{\mathbf{J}}}^{⊺}\hat{\mathbf{J}}(\mathbf{\theta }-\hat{\mathbf{\theta }})\le p{s}^2{F}_{\alpha }(p,n-p)\}$ where $s^2=\frac{{\sum }_{i=1}^n(Y_i-f({\mathbf{X}}_i,\hat{\mathbf{\theta }}){)}^2}{n-p}$, $p-1$ is the number of explanatory variables, and $\hat{\mathbf{J}}=\left(\frac{\partial f({\mathbf{X}}_i,\mathbf{\theta })}{\partial {\theta }_j}\right){|}_{\mathbf{\theta }=\hat{\mathbf{\theta }}}$ is the Jacobian Matrix of $f$.

Link to original

Confidence Interval for Parameter of Nonlinear Regression

Definition

Suppose a Nonlinear Regression Model $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ The $100(1-\alpha )\%$ confidence interval for ${\theta }_j$ is obtained as ${\hat{\theta }}_j\pm t_{\alpha /2}(n-p)SE({\hat{\theta }}_j)$ where $SE({\hat{\theta }}_j)$ is the Standard Error for Parameter of Nonlinear Regression ${\theta }_j$, and $p-1$ is the number of explanatory variables.

Link to original

Hypothesis testing of Parameter of Nonlinear Regression

Definition

Suppose a Nonlinear Regression Model $Y_i=f({\mathbf{X}}_i,\mathbf{\theta })+{ϵ}_i,i=1,2,\dots ,n$ The null hypothesis $H_0:{\theta }_j={\theta }_{j,0}$ can be tested with the test statistic $t_j=\frac{{\hat{\theta }}_j-{\theta }_{j,0}}{SE({\hat{\theta }}_j)}\sim t(n-p)$ where $SE({\hat{\theta }}_j)$ is the Standard Error for Parameter of Nonlinear Regression ${\theta }_j$, and $p-1$ is the number of explanatory variables.

Link to original

Non-Parametric Regression

Kernel Estimator

K-Nearest Neighbors Algorithm

Definition

Suppose $d_i(x):=d(x,X_i)$ is the distance between the training sample $X_i$ and the given point $x$. Let $d_{(1)}(x)\le d_{(2)}(x)\le \cdots \le d_{(n)}(x)$ $\underset{1\le k\le K}{\mathrm{argmax}}{\hat{f}}_k=\underset{1\le j\le K}{\mathrm{argmax}}{p}_j(x)=\underset{1\le j\le K}{\mathrm{argmax}}\frac{1}{k}\sum _{i=1}^{n}I(d_i(x)\le d_{(k)}(x))I(Y_i=j)$

The k-NN determines a sample’s class using the $k$-nearest training datum.

Facts

$k$ and $d$ are hyperparameters.

Link to original

Kernel Density Estimation

Definition

${\hat{f}}_h(x)=\frac{1}{n}\sum _{i=1}^n{K}_h(x-X_i)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-X_i}{h}\right)$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(t)dt=1,{\int }_{-\infty }^{\infty }{K}^2(t)dt<\infty$ is the kernel, $K_h(x):=\frac{1}{h}K(\frac{x}{h})$ is the scaled kernel, and $h$ is a smoothing parameter (bandwidth, or window width)

KDE is a non-parametric method to estimate PDF of a random variable based on kernel and wrights.

Multidimensional KDE

${\hat{f}}_H(x)=\frac{1}{n}\sum _{i=1}^n{K}_H(\mathbf{x}-{\mathbf{x}}_i)=\frac{1}{n}|H{|}^{-1}\sum _{i=1}^{n}K(H^{-1}(\mathbf{x}-{\mathbf{x}}_i))$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(\mathbf{u})d\mathbf{u}=1,K\ge 0$ is the kernel, $K_H(x):=|H{|}^{-1}K(H^{-1}x)$ is the scaled kernel, and $H$ is a positive-definite bandwidth matrix.

Facts

Widely used kernels

Link to original

Kernel Regression

Definition

Kernel regression is an extension of Kernel Density Estimation to estimate the conditional expectation of a random variable $Y$ given $X$. The kernel regression estimator is defined as ${\hat{f}}_{\lambda }(x)=\frac{1}{n}\sum _{i=1}^n{K}_{\lambda }(x-X_i)Y_i=\frac{1}{n\lambda }\sum _{i=1}^{n}K(\frac{x-X_i}{\lambda })Y_i$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(t)dt=1,{\int }_{-\infty }^{\infty }{K}^2(t)dt<\infty$ is the kernel, $K_{\lambda }(x):=\frac{1}{h}K(\frac{x}{h})$ is the scaled kernel, and $\lambda$ is a smoothing parameter

Link to original

Nadaraya–Watson Kernel Regression

Definition

Nadaraya–Watson kernel regression is a Kernel Regression with the sum of the weights is $1$. The Nadaraya–Watson kernel regression estimator is defined as ${\hat{f}}_{\lambda }(x)=\frac{{\sum }_{j=1}^{n}K(\frac{x-X_i}{\lambda })Y_i}{{\sum }_{j=1}^{n}K(\frac{x-X_i}{\lambda })}$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(t)dt=1,{\int }_{-\infty }^{\infty }{K}^2(t)dt<\infty$ is the kernel, and $\lambda$ is a smoothing parameter

Facts

Nadaraya–Watson kernel regression estimator is the same as Local Polynomial Regression with $p=0$

Link to original

Local Polynomial Regression

Definition

Local polynomial regression fits the polynomial function with the data near the given point. The polynomial function has a form similar to Taylor Series to exploit the relation between the coefficients and the differentiation at the given point.

The $p$-th degree polynomial function at a point $x_0$ is defined as $f(x_0)=\sum _{i=0}^p{\beta }_i(x_0-x{)}^i={\beta }_0+{\beta }_1(x_0-x)+\cdots +{\beta }_p(x_0-x{)}^p$ where $p\in {\mathbb{N}}_0$

Local polynomial regression estimator is obtained by the relation $\hat{f}(x)={\hat{\beta }}_0$. The coefficients are estimated by an optimization problem $\hat{\mathbf{\beta }}=\underset{\mathbf{\beta }}{\mathrm{argmin}}\sum _{i=1}^n\left([Y_i-\sum _{j=0}^p{\beta }_i(X_i-x{)}^j{]}^{2}K\left(\frac{X_i-x}{\lambda }\right)\right)$ where $\mathbf{\beta }=({\beta }_0,\dots ,{\beta }_p)$, $K$ is a kernel function, and $\lambda$ is a smoothing parameter.

In a matrix notation, the optimization problem can be written as $\hat{\mathbf{\beta }}=\underset{\mathbf{\beta }}{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^{⊺}\mathbf{W}(\mathbf{y}-\mathbf{X}\mathbf{\beta })$ where $\mathbf{y}=(Y_1,\dots ,Y_n{)}^{⊺}$, $\mathbf{X}=\left(\begin{array}{ccc}(X_1-x{)}^0& \dots & (X_1-x{)}^p\\ ⋮& \ddots & ⋮\\ (X_n-x{)}^0& \dots & (X_n-x{)}^p\end{array}\right)$ is an $n\times (p+1)$ matrix, and $\mathbf{W}=\mathrm{diag}\left[K\left(\frac{X_1-x}{\lambda }\right),\dots ,K\left(\frac{X_n-x}{\lambda }\right)\right]$

A solution of the problem is obtained as ${\hat{f}}^{(j)}(x)=j!{\hat{\beta }}_j=j!1_j^{⊺}({\mathbf{X}}^{⊺}\mathbf{WX}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{Wy}$ where $1_j$ is a $(p+1)$-vector whose $j$-th element is $1$ and the others are $0$.

Link to original

Series Estimator

Series Estimator

Definition

Series Estimator

Suppose a non-parametric regression function $y_i=f(x_i)+{ϵ}_i$ where $f(x)\in L_2[a,b]$ and let a sequence of functions $\{g_j(x)\}$ be a complete orthonormal basis. Then, the function can be written as a series. $f(x)=\sum _{j=1}^{\infty }{\theta }_j{g}_j(x)$ where ${\theta }_j=\int g_j(x)f(x)dx$ is a parameter subject to estimation.

The regression function can be approximated with the estimated parameter ${\hat{\theta }}_j$ and finite sum. It is called a series estimator ${\hat{f}}_{\lambda }(x)=\sum _{j=1}^{\lambda }{\hat{\theta }}_j{g}_j(x)$ where $\lambda \in \mathbb{N}$ is a smoothing parameter.

The series estimator can be written as a form of a kernel estimator with the estimation of the parameter ${\hat{\theta }}_j=\frac{1}{n}\sum _{i=1}^n{g}_j(x_i)Y_i$. ${\hat{f}}_{\lambda }(x)=\sum _{j=1}^{\lambda }{\hat{\theta }}_j{g}_j(x)=\frac{1}{n}\sum _{i=1}^n\left(\sum _{j=1}^{\lambda }{g}_j(x_i)g_j(x)\right)Y_i=\frac{1}{n}\sum _{i=1}^n{K}_{\lambda }(x,x_i)Y_i$ where $K_{\lambda }$ is a kernel.

Thresholding

In general, as the $j$ is larger the $|{\theta }_j|$ smaller. So the terms of the series can be selected with a certain algorithm.

The hard-thresholding estimator is defined as ${\hat{f}}_H(x)=\sum _{j=1}^{\infty }{\hat{\theta }}_j{I}_{\{|{\hat{\theta }}_j|>\delta \}}\cdot g_j(x)$ where $\delta$ is a thresholding parameter

And the soft-thresholding estimator is defined as ${\hat{f}}_S(x)=\sum _{j=1}^{\infty }\mathrm{sign}({\hat{\theta }}_j)(|{\hat{\theta }}_j|-\delta {)}_{+}\cdot g_j(x)$

Link to original

Multiresolution Analysis

Definition

Multiresolution analysis decomposes a signal into a coarse approximation and a series of details at different scales. $f(t)=\sum _{k=-\infty }^{\infty }{c}_k{\varphi }_k(t)+\sum _{j=-\infty }^{\infty }\sum _{k=-\infty }^{\infty }{d}_{j,k}{\psi }_{j,k}(t)$ where the first term $\sum _{k=-\infty }^{\infty }{c}_k{\varphi }_k(t)$ represents the coarse approximation and the second term $\sum _{j=-\infty }^{\infty }\sum _{k=-\infty }^{\infty }{d}_{j,k}{\psi }_{j,k}(t)$ represents the sum of details. Here, ${\varphi }_k(t)$ is a scaling function (father wavelet), ${\psi }_{j,k}(t)$ is a wavelet function (mother wavelet), $c_k$ is approximation coefficient, $d_{j,k}$ is a detail coefficient.

The set $\{{\varphi }_k(x),{\psi }_{j,k}(x),j,k\in \mathbb{Z}\}$ is an orthonormal basis of the $L^2$ function space.

Link to original

Wavelet Estimator

Definition

Wavelet estimator is a kind of Series Estimator using wavelet basis. $\hat{f}(t)=\sum _{k=-\infty }^{\infty }{\hat{c}}_k{\varphi }_k(t)+\sum _{j=-\infty }^{\infty }\sum _{k=-\infty }^{\infty }{\hat{d}}_{j,k}{\psi }_{j,k}(t)$ where ${\hat{c}}_k=\frac{1}{n}\sum _{i=1}^n{\varphi }_k(x_i)Y_i$ and ${\hat{d}}_{j,k}=\frac{1}{n}\sum _{i=1}^n{\psi }_{j,k}(x_i)Y_i$

To find a optimal resolution Thresholding is used.

Link to original

Spline Estimator

Piecewise Polynomials

Definition

Suppose that knots are ${\xi }_j,j=1,\dots ,K$. Then, the basis functions for order-M spline are defined as: $\begin{array}{rl}{h}_j(X)& =X^{j-1},j=1,\dots ,M\\ h_{M+l}(X)& =(X-{\xi }_l{)}_{+}^{M-1},l=1,\dots ,K\end{array}$ where $(x{)}_{+}=\{\begin{array}{ll}x& x\ge 0\\ 0& x<0\end{array}$

Examples

For cubic Polynomials with two knots ${\xi }_1,{\xi }_2$, the basis functions are: $h_1(X)=1,h_2(X)=X,h_3(X)=X^2,h_4(X)=X^3,h_5(X)=(X-{\xi }_1{)}_{+}^3,h_6(X)=(X-{\xi }_2{)}_{+}^3$

Facts

Order of piecewise polynomials

• constant: 1
• linear: 2
• quadratic: 3
• cubic: 4

Link to original

Natural Cubic Splines

Definition

Suppose that knots are ${\xi }_j,j=1,\dots ,K$. Then, the basis functions for natural cubic Spline are defined as: $N_1(X)=1,N_2(X)=X,N_{k+2}(X)=d_k(X)-d_{K-1}(X)$ where $d_k(X)=\frac{(X-{\xi }_k{)}_{+}^3-(X-{\xi }_K{)}_{+}^3}{{\xi }_K-{\xi }_k}$

Link to original

B-Spline

Definition

B-Splines are basis functions for spline function in which order is the same, meaning that all possible spline function can be expressed as a linear combination of B-splines.

Consider the knots $t_0,t_1,\dots t_m$ sorted into non-decreasing order. For the given sequence of knots, the B-splines of degree $p=0$ is defined as $B_{j,0}(x)=\{\begin{array}{ll}1& t_j\le x\le t_{j+1}\\ 0& \text{otherwise}\end{array}$ And the higher $(p+1)$-degree B-splines are defined by recursion $B_{j,p}(x)=\frac{x-t_j}{t_{j+p}-t_j}{B}_{j,p-1}(x)+\frac{t_{j+p+1}-x}{t_{j+p+1}-t_{j+1}}{B}_{j+1,p-1}(x)$

The order $n$ spline function on a given set of knots can be expressed as $s(x)=\sum _{j=1}^m{\gamma }_j{B}_{j,p}(x)$

Link to original

Smoothing Splines

Definition

$RSS(f,\lambda )=\sum _{i=1}^N(y_i-f(x_i){)}^2+\lambda \int f^{\prime \prime }(t{)}^{2}dt$ where $\lambda \ge 0$ is a smoothing parameter. If $\lambda =0$, the smoothing spline is any interpolating spline. If $\lambda =\infty$, the smoothing spline is Least Square line.

Smoothing splines is a Spline basis method that avoids the knot selection problem. Among all functions $f(x)\in C^2$, find one that minimizes the penalized RSS.

Facts

The unique analytic solution of smoothing spline is a Natural Cubic Splines with knots at the unique values of the $x_i$‘s

Link to original

Regression Model for Censored Data

Survival Function and Hazard Function

Survival Function

Definition

Let $T\ge 0$ be a non-negative Random Variable representing survival time with PDF $f(t)$ and CDF $F(t)$. The survival function of the random variable is defined as $S(t)=1-F(t)=P(T>t)$

The survival function is a function that gives the probability that a patient, device, or other object of interest will survive past a certain time.

Link to original

Hazard Function

Definition

Hazard Function

Let $T\ge 0$ be a non-negative Random Variable representing survival time with PDF $f(t)$ and CDF $F(t)$. The hazard function $0<\lambda (t)<\infty$ of the random variable $T$ is defined as $\lambda (t)=\frac{f(t)}{S(t)}={\lim }_{\Delta t\to 0}\frac{P(T<t+\Delta t|T>t)}{\Delta t}$ where $S(t)$ is the Survival Function.

The hazard function refers to the rate of occurring event at a given time $t$.

Cumulative Hazard Function

The hazard function can alternatively be represented in terms of the cumulative hazard function, defined as $\Lambda (t)={\int }_0^t\lambda (u)du$

Facts

The Survival Function $S(t)$, the Cumulative Hazard Function, the density (PDF) $f(t)$, the Hazard Function, and the distribution function (CDF) of survival time $F(t)$ are related through $S(t)=\exp [-\Lambda (t)]=\frac{f(t)}{\lambda (t)}=1-F(t)$

Link to original

Censored Data

Right Censoring

Kinds

Type 1 Censoring

Definition

Type 1 Censoring

when $c_1=\cdots =c_n=t_c$ ($c_i$’s are constant)

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. We observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

In a type 1 censoring setting, the censored times $C_i$ are fixed constants, not random variables.

The PDF of the observation $Y_i$ is derived as $\begin{array}{rl}P(Y_i=C_i,{\delta }_i=0)& =P(T_i>C_i)=S(y_i)\\ P(Y_i=T_i,{\delta }_i=1)& =f(y_i)\end{array}$

Likelihood of Type 1 Censoring Data

The Likelihood Function of type 1 censoring data is defined as $L(\theta )=\prod _{i=1}^{n}f(y_i{)}^{{\delta }_i}S(y_i{)}^{1-{\delta }_i}=\prod _{i=1}^n\lambda (y_i{)}^{{\delta }_i}S(y_i)=\prod _{i=1}^n\lambda (y_i{)}^{{\delta }_i}\exp [-\Lambda (y_i)]$

Link to original

Type 2 Censoring

Definition

Type 2 Censoring

when $r=3$

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. We observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

In a type 2 censoring setting, we observe first $r$ out of $n$ experiment $T_1,\dots ,T_r$. In other words, for the order statistics $T_{(i)}$ of $T_i$, we only observe $T_{(1)}\le T_{(2)}\le \cdots \le T_{(r)}$. Where $C_1=C_2=\cdots =C_n$ are not constants, but random variables.

Likelihood of Type 2 Censoring Data

The likelihood function of type 2 censored data can be computed using the same equation used for type 1 censored data but computing the joint-PDF of the Order Statistic $T_{(1)}\le \cdots \le T_{(r)}$ is easier.

The joint PDF of $T_{(1)}\le \cdots \le T_{(r)}$ is derived as $L(\theta )=\frac{n!}{(n-r)!}\left[\prod _{i=1}^{r}f(t_{(i)})\right]S(t_{(r)}{)}^{n-r}$

Link to original

Random Censoring

Definition

Random Censoring

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times, which can be both random variable or constant. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator.

In a random censoring setting, we only observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$, and the censoring times are $C_1,\dots ,C_n$ are i.i.d. Random Variable follows PDF $g_{\gamma }$ and CDF $G_{\gamma }$.

The PDF of the observation $Y_i$ is derived as $\begin{array}{rl}P(Y_i=t,{\delta }_i=0)& =P(C_i=t,T_i>C_i)=g_{\gamma }(t)S_{\theta }(t)\\ P(Y_i=t,{\delta }_i=1)& =P(T_i=t,T_i\le C_i)=f_{\theta }(t)(1-G_{\gamma }(t))\end{array}$ where $\theta$ is the parameter of interest, and $\gamma$ is the: nuisance parameter

Likelihood of Random Censoring Data

The Likelihood Function of random censored data is defined as $L(\theta )=\prod _{i=1}^n[f_{\theta }(t)(1-G_{\gamma }(t)){]}^{{\delta }_i}[S_{\theta }(t)g_{\gamma }(t){]}^{1-{\delta }_i}=c\cdot \prod _{i=1}^n{f}_{\theta }(t{)}^{{\delta }_i}{S}_{\theta }(y_i{)}^{1-{\delta }_i}$ where $c$ is a constant.

Link to original

Link to original

Left Censoring

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times, and $C_1,C_2,\dots ,C_n$ are the censoring times, which can be both random variable or constant. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\ge C_i\text{: uncensored}\\ 0& T_i<C_i\text{: censored}\end{array}$ be the censoring indicator.

In a left censoring setting, we observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\max (T_i,C_i)$. In other words the event has already occurred before it becomes the subject of observed.

Link to original

Interval-Censored Data

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times. Interval censored data is given as an interval, not an exact point of time. We only observe the interval $(L_i,R_i]$ that includes $T_i$ Interval-censored data is divided into four cases.

Case 1 Interval-Censored Data (Current Status Data)

Data is given the form of $[0,C_i]$ or $(C_i,\infty )$, where $C_i$ is a fixed time point.

Case 2 Interval-Censored Data

Data is given the form of $(L_i,R_i]$ where $0<L_i<R_i<\infty$.

Double Censored Data

Data is given the form of $(L_i,R_i]$ where $0\le L_i\le R_i\le \infty$.

If $R_i=\infty$, then it is right-censored data, if $L_i=0$, then it is left-censored data, and if $L_i=R_i$, then it is uncensored data.

Panel Data

Observations are made at discrete time points. The period between these observations can be viewed as an interval.

Link to original

Mean Imputation Method

Definition

Mean imputation method substitute the given interval of case 2 interval-censored data $(L_i,R_i]$ with the mean of the interval $Y_i=\frac{(L_i+R_i)}{2}$ if the interval is finite. If $R_i=\infty$, the data $(L_i,\infty )$ is substituted with the left value $Y_i:=L_i$ and treated as right-censored data.

Link to original

Estimation of Survival Function

Kaplan-Meier Estimator

Definition

Kaplan-Meier Estimator

Consider a Random Censoring case $Y_i=\min (T_i,C_i)$. Assume that $y_1\le y_2\le \cdots \le y_n$, and the distinct failure times are ${\tau }_1<{\tau }_2<\cdots <{\tau }_k$ where $k\le n$. Let $d_j:=\sum _{i=1}^{n}I(y_i={\tau }_j,{\delta }_i=1)$ be the number of deaths at ${\tau }_j$, and $n_j:=\sum _{i=1}^{n}I(y_i\ge {\tau }_j)$ be the number of alive at ${\tau }_j$ where the set$\{y_i|y_i\ge {\tau }_j\}$ is called the risk set at ${\tau }_j$. We only observe $(Y_i,{\delta }_i),i=1,\dots ,n$

The Kaplan-Meier estimator is derived from the expression $S(t)=\prod _{i=1}^k{P}_j,P_j=P(T>{\tau }_j|T>{\tau }_{j-1})$ where ${\tau }_0=0$

General Case

As an estimator of $P_j$, consider ${\hat{P}}_j=\hat{P}(T>{\tau }_j|T>{\tau }_{j-1})=\frac{\hat{P}(T>{\tau }_j)}{\hat{P}(T>{\tau }_{j-1})}=\frac{n_{j+1}/n_1}{n_j/n_1}=\frac{n_j-d_j}{n_j}=1-\frac{d_j}{n_j}$

The Kaplan-Meier estimator is defined with the estimated ${\hat{P}}_j$‘s $\hat{S}(t)=\prod _{j=1}^k{\hat{P}}_j=\prod _{\{j|{\tau }_j\le t\}}\left(1-\frac{d_j}{n_j}\right)$ The cumulative hazard function is estimated by Nelson-Aalen Estimator in the same logic $\hat{\Lambda }(t)=\sum _{\{j|{\tau }_j\le t\}}\frac{d_j}{n_j}$

No Ties Case

When there’s no tie in the observation, $y_1<y_2<\cdots <y_n$, then the failure times are equal to the observation ${\tau }_j=y_i$, death is equal to the censoring indicator $d_j={\delta }_j$, and $n_j=n-j+1$. Thus, the Kaplan-Meier estimator is defined as $\hat{S}(t)=\prod _{\{j|y_j\le t\}}(1-\frac{{\delta }_j}{n-j+1})=\prod _{\{j|{\tau }_j\le t\}}(1-\frac{1}{n-j+1}{)}^{{\delta }_j}$

Properties

Self-Consistency

An estimator $\hat{S}(t)$ is self-consistent if $\begin{array}{rl}\hat{S}(t)& =\frac{1}{n}\sum _{i=1}^n\left[1I(Y_i>t)+0I(Y_i\le t,{\delta }_i=1)+\frac{\hat{S}(t)}{\hat{S}(Y_i)}I(Y_i<t,{\delta }_i=0)\right]\\ & =\frac{1}{n}\left[N_y(t)+\sum _{i:y_i\le t}(1-{\delta }_i)\frac{\hat{S}(t)}{\hat{S}(Y_i)}\right]\end{array}$ where $N_y(t)=\sum _{i=1}^{n}I(Y_i>t)=\#(Y_i>t)$

The Kaplan-Meier estimator is the unique self-consistent estimator for $t<Y_{(n)}$ where $Y_{(n)}$ is the largest observation.

Generalized MLE

The Kaplan-Meier estimator gives the Generalized Maximum Likelihood Estimation of the Survival Function $S$.

Strong Consistency

The Kaplan-Meier estimator uniformly Almost Surely converges to $S(t)$ $\hat{S}(t)\stackrel{a.s.}{\to }S(t),\forall t\in {\mathbb{R}}^{+}$

Proof

Consider a function $S^{\ast }(t)=P(Y>t)$ and decompose it to the sum of the subsurvival functions $S_u^{\ast }(t)$ and $S_c^{\ast }(t)$. $S^{\ast }(t)=S_u^{\ast }(t)+S_c^{\ast }(t)$ where $S_u^{\ast }(t)=P(Y>t,\delta =1)$ is the uncensored case and $S_c^{\ast }(t)=P(Y>t,\delta =0)$ is the censored case.

Then, the survival function $S(t)=P(T>t)$ can be expressed as a function of the subsurvival functions. $S(t)=\Psi (S_u^{\ast },S_c^{\ast },t)$

Define the empirical subsurvival functions ${\hat{S}}_u^{\ast }(t)=\frac{1}{n}\sum _{i=1}^{n}I(Y_i>t,{\delta }_i=1)$ and ${\hat{S}}_c^{\ast }(t)=\frac{1}{n}\sum _{i=1}^{n}I(Y_i>t,{\delta }_i=0)$. The Kaplan-Meier estimator also can be expressed as a function of the empirical subsurvival functions. $\hat{S}(t)=\Psi ({\hat{S}}_u^{\ast },{\hat{S}}_c^{\ast },t)$

By Glivenko-Cantelli theorem, ${\hat{S}}_u^{\ast }(t)\stackrel{a.s.}{\to }{S}_u^{\ast }(t)$ and ${\hat{S}}_c^{\ast }(t)\stackrel{a.s.}{\to }{S}_c^{\ast }(t)$ for all $t\in {\mathbb{R}}^{+}$. Since $\Psi$ is a continuous function of $S_u^{\ast }(t)$ and $S_c^{\ast }(t)$, $\hat{S}(t)=\Psi ({\hat{S}}_u^{\ast },{\hat{S}}_c^{\ast },t)\stackrel{a.s.}{\to }\Psi (S_u^{\ast },S_c^{\ast },t)=S(t),\forall t\in {\mathbb{R}}^{+}$

Asymptotic Normality

Kaplan-Meier estimator has asymptotic normality. $\sqrt{n}(\hat{S}(t)-S(t))\stackrel{D}{\to }N\left(0,S(t{)}^2{\int }_0^t\frac{d{F}_u(X)}{(1-H(x){)}^2}\right)$ where $F_u(t)=P(T\le t,\delta =1)={\int }_0^t(1-G(x))dF(x)$, $C\sim G$, and $Y\sim H$.

The variance of the estimator is estimated by Greenwood’s formula ${\sigma }_S^2(t)\approx \hat{S}(t{)}^2\sum _{\{j|{\tau }_j\le t\}}\frac{d_j}{n_j(n_j-d_j)}$ For the no ties case, the formula is ${\sigma }_S^2(t)\approx \hat{S}(t{)}^2\sum _{\{j|{\tau }_j\le t\}}\frac{{\delta }_j}{(n-j)(n-j+1)}$

Examples

$n=8$ case where $\{\begin{array}{ll}\times & {\delta }_i=1\\ \circ & {\delta }_i=0\end{array}$

Facts

Kaplan-Meier estimator has Self-Consistency and Asymptotic Normality, and it is generalized MLE

If no censoring, Kaplan-Meier estimator is just the Empirical Survival Function.

Link to original

Kernel Estimation for Survival Analysis

Definition

Kaplan-Meier estimator is a step function. So it is difficult to calculate its quantile function and Density Function. The Kernel Density Estimation is used to make it smooth function.

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time with a distribution $F$, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time with a distribution $G$. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$ where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator.

Without Censoring

For the complete data, the Kernel Density Estimation is defined as $\hat{f}(t)=\frac{1}{n}\sum _{i=1}^n{K}_{\lambda }(t-Y_i)=\frac{1}{n\lambda }\sum _{i=1}^{n}K\left(\frac{t-Y_i}{\lambda }\right)$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(t)dt=1,{\int }_{-\infty }^{\infty }{K}^2(t)dt<\infty$ is the kernel, $K_{\lambda }(x):=\frac{1}{h}K(\frac{x}{h})$ is the scaled kernel, and $\lambda$ is a smoothing parameter.

The kernel estimator for the Distribution Function is defined as $\hat{F}(t)={\int }_{-\infty }^t\hat{f}(u)du=\frac{1}{n}\sum _{i=1}^{n}W\left(\frac{t-Y_i}{\lambda }\right)$ where $W(t)={\int }_{-\infty }^{t}K(u)du$

With Censoring

For the censored data, the weights for each observation is defined as a jump size in Kaplan-Meier Estimator. $\hat{f}(t)=\sum _{i=1}^n{S}_i{K}_{\lambda }(t-Y_i)$ $\hat{F}(t)=\sum _{i=1}^n{S}_{i}W\left(\frac{t-Y_i}{h}\right)$ where $S_i$ is the jump size at $Y_i$ in Kaplan-Meier Estimator.

Thus, the Kernel Density Estimation for the Survival Function is $1-\sum _{i=1}^n{S}_{i}W\left(\frac{t-Y_i}{h}\right)$

Link to original

Cox Regression Model

Cox Proportional Hazards Model

Definition

Cox proportional hazards model assume that covariates affect the Hazard Function.

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator, and have covariates ${\mathbf{x}}_i=(x_{i1},\dots ,x_{ip})$. Then, the Cox proportional hazards model is defined as $\lambda (t:\mathbf{x})={\lambda }_0(t)\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })$ where ${\lambda }_0$ is called the baseline hazard function, i.e. hazard at $\mathbf{x}=0$

Conditional Likelihood

Let $Y_1<Y_2<\cdots <Y_n$ (no ties case), and ${\mathcal{R}}_i$ be the risk set. For each uncensored time $Y_i$, $P\{\text{a death in}[y_i,y_i+\Delta )|{\mathcal{R}}_i\}\approx \sum _{j\in {\mathcal{R}}_i}{\lambda }_0(y_j)\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })\Delta$ Therefore, $P\{\text{a death of }i\text{ at time}{y}_i|\text{one death in }{\mathcal{R}}_i\text{ at time }{y}_i\}=\frac{{\lambda }_0(y_i)\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })\Delta }{\sum _{j\in {\mathcal{R}}_i}{\lambda }_0(y_i)\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })\Delta }\approx \frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}$ Taking the product of these conditional probabilities gives a conditional likelihood $L_C(\mathbf{\beta })={\prod }_{i:I_u}\frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}={\prod }_{i=1}^n{\left[\frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}\right]}^{{\delta }_i}$ where $I_u$ is the indicator set for uncensored samples.

The $L_C(\mathbf{\beta })$ is not a likelihood. However, Cox suggested treating the conditional likelihood as an ordinary likelihood to find the Maximum Likelihood Estimation.

Since there’s no analytic solution for the MLE, iterative methods such as Newton–Raphson method is used to estimate the coefficient $\mathbf{\beta }$.

The hazard ratio $\exp ({\hat{\beta }}_j)$ represents the relative change in Hazard Rate for a one-unit increase in the covariate $x_j$.

Goodness-of-Fit Test

For testing the null hypothesis $H_0:\mathbf{\beta }=0$, Cox suggested the Rao Test. ${\left(\frac{\partial \ln L_C(0)}{\partial \mathbf{\beta }}\right)}^{⊺}{\left(-\frac{{\partial }^2\ln L_C(0)}{\partial {\mathbf{\beta }}^2}\right)}^{-1}\left(\frac{\partial \ln L_C(0)}{\partial \mathbf{\beta }}\right)\stackrel{a}{\sim }{\chi }^2(p)$

Asymptotic Normality of MLE

$\hat{\mathbf{\beta }}\stackrel{D}{\to }{N}_p(\beta ,I^{-1}(\mathbf{\beta }))$ where $I(\mathbf{\beta })$ is the observed Fisher Information

Estimation of Survival Function

Under the Cox proportional hazards model, $S(t;\mathbf{x})=S_0(t{)}^{\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })}$ To estimate $S(t;\mathbf{x})$, we can use $\hat{\mathbf{\beta }}$ for $\mathbf{\beta }$ but we still need to estimate $S_0(t)$, ${\lambda }_0(t)$, or ${\Lambda }_0(t)$.

Breslow suggested the estimators of ${\lambda }_0(t)$ and $S_{0,B}(t)$ as ${\hat{\lambda }}_{0,B}=\frac{1}{(y_{u_i}-y_{u_{i-1}})\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}$ If $Y_{u_{i-1}}<t<Y_{u_i}$

${\hat{S}}_{0,B}(t)={\prod }_{i:y_i\le t}\left(1-\frac{{\delta }_i}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}\right)$ If $\beta =0$, then ${\hat{S}}_{0,B}$ is the Kaplan-Meier Estimator

It has a few drawbacks

• ${\hat{S}}_{0,B}(t)\ne \exp (-{\hat{\Lambda }}_0(t))$
• ${\hat{S}}_{0,B}(t)$ can take negative values.

Tsiatis suggested a non-negative version of ${\hat{S}}_{0,B}(t)$ ${\hat{S}}_{0,T}(t)=\exp (-{\hat{\Lambda }}_{0,T}(t))$ where ${\hat{\Lambda }}_{0,T}(t)=\sum _{y_i\le t}\frac{{\delta }_i}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}$

Link suggested using the linear smooth of ${\hat{\Lambda }}_{0,T}(t)$.

Discrete on Grouped Data

When data is discrete or grouped, there are ties at each failure. Denote the ordered discrete failure time by $Y_1<Y_2<\cdots <Y_r$ and let ${\mathcal{R}}_i$ be the risk set at $Y_i^{-}$, ${\mathcal{D}}_i$ be the death set at $Y_i$, and $d_i=\#(D_i)$.

Cox suggested combining the all possible permutations. However, it is computationally infeasible. $L_C={\prod }_{i=1}^r\left[\frac{\prod _{j\in {\mathcal{D}}_i}{\psi }_j}{\sum _{{\mathcal{D}}_i^{\ast }}\prod _{j\in {\mathcal{D}}_i^{\ast }}{\psi }_j}\right]$ where ${\psi }_j=\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })$, and ${\mathcal{D}}_i^{\ast }$ is the size $d_i$ subset of ${\mathcal{R}}_i$

Peto suggested an alternative likelihood that instead of all possible permutations, use the same contribution. $L_C={\prod }_{i=1}^r\left[\frac{\prod _{j\in {\mathcal{D}}_i}{\psi }_j}{(\sum _{j\in {\mathcal{R}}_i}{\psi }_j{)}^{d_i}}\right]$

Time Dependent Covariates

In the case, the covariate depends on time. We observe ${\mathbf{x}}_i(t)$ and the conditional likelihood defined as $L_C(\mathbf{\beta })={\prod }_{i:I_u}\frac{\exp ({\mathbf{x}}_i(y_i{)}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j(y_i{)}^{⊺}\mathbf{\beta })}$ where $I_u$ is the indicator set for uncensored samples, and ${\mathcal{R}}_i$ is the risk set.

Facts

Any two individuals have hazard functions that are constant multiples of the one another.

The Survival Function of the Cox proportional hazard model is a family of Lehmann alternatives. $\forall S\in \mathcal{S},\exists S_0\in \mathcal{S}s.t.S=S_0^{\gamma }$ where $\gamma \in {\mathbb{R}}^{+}$.

If $p=1$, $x_i=I(i\in I_1)$, where $I_1$ is the indicator set for sample 1, and there are no ties, then Cox test is exactly equal to the Mantel-Haenszel Test.

Link to original

Linear Regression Model

Accelerated Life Model

Definition

Consider the random variable $T_0$ represents survival time with $\mathbf{x}=0$ Hazard Function with $f_0(t)$, ${\lambda }_0(t)$, and $S_0(t)$ And assume that the survival time of individual with covariate $\mathbf{x}$ is defined as $T_{\mathbf{x}}=T_0\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })$ If ${\mathbf{x}}^{⊺}\mathbf{\beta }<0$ then the covariate $\mathbf{x}$ accelerates the time to failure. The model based on this assumption is called accelerated failure time (AFT) model.

Under AFT model $S(t:\mathbf{x})=S_0(t\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta }))$ $\lambda (t:\mathbf{x})={\lambda }_0(t\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta }))\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta })$

Let $Y=\ln T_{\mathbf{x}}$, then $E(Y)=E(\ln T_0)+{\mathbf{x}}^{⊺}\mathbf{\beta }\equiv \alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }$ $Y=\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }+ϵ\equiv E(\ln T_0)+{\mathbf{x}}^{⊺}\mathbf{\beta }+\underset{ϵ}{\underbrace{\ln T_0-E(\ln T_0)}}$ where $ϵ=\ln T_0-E(\ln T_0)$.

Assume that $ϵ=\sigma W$, where $W$ is a Random Variable represents error term. Then, the AFT model becomes $Y\equiv \ln T_{\mathbf{x}}=\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }+\sigma W$

Relationship between $Y=\ln T_{\mathbf{x}}$ and $W$. $F_T(t)=P_W(\frac{\ln t-\mu }{\sigma })$ $f_T(t)=\frac{1}{\sigma t}{f}_W(\frac{\ln t-\mu }{\sigma })$ where $\mu =\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }$

If $W\sim N(0,1)$ then $\ln T\sim N(\mu ,{\sigma }^2)$, and if $W\sim \mathrm{Gumbel}(0,1)$, then $T\sim \mathrm{Exp}(\lambda )$ where $\lambda =\exp [-(\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta })]$.

Link to original

Miller Estimator

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator

Suppose a Simple Linear Regression model $Y_i=\alpha +\beta X_i+{ϵ}_i,i=1,2,\dots ,n$

With no censoring present, the least squares estimators of the parameters $\alpha ,\beta$ are obtained by minimizing $\frac{1}{n}\sum _{i=1}^n{ϵ}_i^2=\frac{1}{n}\sum _{i=1}^n(T_i-\alpha -\beta x_i{)}^2={\int }_{-\infty }^{\infty }{z}^{2}d{F}_n(z)$ where $F_n(z):=\frac{1}{n}\sum _{i=1}^{n}I({ϵ}_i\le z)$ is the Empirical Distribution Function of $z_1,z_2,\dots ,z_n$ where $z_i=y_i-\alpha -\beta x_i$.

With censoring present, Miller proposed to minimize ${\int }_{-\infty }^{\infty }{z}^{2}d{\hat{F}}_n(z)=\sum _{i=1}^n{\hat{w}}_i(\beta )(Y_i-\alpha -\beta x_i{)}^2$ where $\hat{F}$ is the Kaplan-Meier Estimator based on $(z_i,{\delta }_i)$ and the weights ${\hat{w}}_i(\beta )$ is its jump size.

If the last observation is censored, then $\sum _{i=1}^n{\hat{w}}_i(\beta )<1$. Hence, change the last observation to be uncensored, so that $\sum _{i=1}^n{\hat{w}}_i(\beta )=1$.

Link to original

Buckley-James Estimator

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator

If we can observe the true survival time $T_i$, we can make a model $E(T_i)=\alpha +\beta x_i$ However, we can’t observe $T_i$, but only censored $Y_i$, and $E(Y_i)\ne \alpha +\beta x_i$ Buckley and James proposed an Unbiased Estimator for $\alpha +\beta x_i$ $Y_i^{\ast }=Y_i{\delta }_i+E(T_i|T_i>Y_i)(1-{\delta }_i)$ Since we also can not observe $y_i^{\ast }$, we estimate it again. ${\hat{y}}_i^{\ast }=y_i{\delta }_i+\left[\hat{\beta }{x}_i+\frac{\sum _{k:{\hat{z}}_k>{\hat{z}}_i}{\hat{w}}_k(\hat{\beta }){\hat{z}}_k}{1-\hat{F}({\hat{z}}_i)}\right](1-{\delta }_i)$ where ${\hat{z}}_i=y_i-\hat{\beta }{x}_i$, $\hat{F}$ is the Kaplan-Meier Estimator based on $({\hat{z}}_i,{\delta }_i)$ and the weights ${\hat{w}}_i(\hat{\beta })$ is its jump size.

The variance of the estimator is estimated by $\widehat{\mathrm{Var}}(\hat{\beta })=\frac{{\hat{\sigma }}_u^2}{\sum _{i:{\delta }_i=1}(x_i-{\bar{x}}_u{)}^2}$ where ${\hat{\sigma }}_u^2=\frac{1}{n_u-2}\sum _{i:{\delta }_i=1}[y_i-{\bar{y}}_u-\hat{\beta }(x_i-{\bar{x}}_u){]}^2$, $n_u=\sum _{i=1}^n{\delta }_i$, ${\bar{x}}_u=\sum _{i:{\delta }_i=1}{x}_i$, and ${\bar{y}}_u=\sum _{i:{\delta }_i=1}{y}_i$

Link to original