Title: Statistical Modelling

Q: Characterising orthogonal projections
A: $\exists V\leq\R^n$ such that $\Pi |_V=\id$ and $\Pi |_{V^\perp}=0$ $\iff$ $\Pi$ is symmetric idempotent.
\begin{proof}
($\Rightarrow$) From the definition $\Pi^T(I-\Pi)=0$. So $\Pi^T=\Pi^T\Pi$ is symmetric idempotent.
($\Leftarrow$) Let $V=\im\Pi$. If $v\in\im\Pi$, say $v=\Pi u$, then by idempotency $\Pi v = \Pi^2 u = \Pi u = v$. Whereas if $v\in (\im\Pi)^\perp$, then by symmetry $0=\Pi^T v = \Pi v$.
\end{proof}

Q: Rank of projection matrix
A: For a real symmetric idempotent matrix $A$
\[
\rank(A) = \tr(A).
\]
\begin{proof}\leavevmode
Simply note that $A$ has orthogonal diagonalisation
\[
Q^TAQ=\diag(\underbrace{1, \cdots, 1}_r, \underbrace{0, \cdots, 0}_{n-r}).
\]
\end{proof}

Q: Correlation of residuals and fitted values
A: \[
\cov(PY, (I-P)Y) = P\sigma^2 I (I-P) =0.
\]

Q: Single OLS coefficient
A: Let $X_j^\perp := (I-P_{-j})X_j$. Then
\[
\hat{\beta}_j = \frac{(X_j^\perp)^T Y}{\|X_j^\perp \|^2}, \quad \var(\hat\beta_j) = \frac{\sigma^2}{\|X_j^\perp \|^2}.
\]
\begin{proof}
Writing $Y=PY+(I-P)Y$ and using $X_j\perp \im(I-P)$:
\begin{align*}
\frac{(X_j^\perp)^T Y}{\|X_j^\perp \|^2} &= \frac{(X_j^\perp)^T PY}{\|X_j^\perp \|^2} \\
&= \frac{(X_j^\perp)^T X \hat\beta }{\|X_j^\perp \|^2} \\
&= \frac{(X_j^\perp)^T X_j \hat\beta_j }{\|X_j^\perp \|^2} \\
&= \frac{\|X_j^\perp \|^2 \hat\beta_j }{\|X_j^\perp \|^2}  \quad\text{(idempotency)}.
\end{align*}
\end{proof}

Q: MSPE of OLS estimator
A: The Mean Squared Prediction Error of the OLS estimator is
\[
\frac{1}{n}\E \| X\beta - X\hat\beta\|^2 = \sigma^2 \frac{p}{n}.
\]
\begin{proof}
We can compute
\begin{align*}
\E \| X\beta - X\hat\beta\|^2 &= \E \|P\varepsilon\|^2 \\
&=\E[\tr(\varepsilon^T P \varepsilon)] \\
&=\E[\tr(\varepsilon\varepsilon^T P )] \\
&=\tr(\sigma^2 P) \\
&=\sigma^2 p.
\end{align*}
\end{proof}

Q: Gauss Markov Theorem
A: The OLS $\hat{\boldsymbol\beta}$ is the BLUE, i.e. for any unbiased $\boldsymbol\beta^*=C Y$
  \[
    \var(\mathbf{t}^T\hat{\boldsymbol\beta}) \leq \var (\mathbf{t}^T\boldsymbol\beta^*)
  \]
  for all $\mathbf{t}\in \R^p$.
\begin{proof}
Since it is unbiased $CX\boldsymbol\beta^*=\boldsymbol\beta^*$. So $CX = I_p$.
  Now let $D= C - (X^TX)^{-1}X^T$, such that $DX=I_p-I_p=0$. Hence the cross-terms vanish leaving
  \begin{align*}
    \cov(\boldsymbol\beta^*) &= \sigma^2CC^T\\
    &= \sigma^2(D + (X^TX)^{-1}X^T)(D + (X^TX)^{-1}X^T)^T\\
    &= \sigma^2(DD^T + (X^TX)^{-1})\\
    &= \sigma^2 DD^T + \cov(\hat{\boldsymbol\beta})
  \end{align*}
where $DD^T$ is postive semi-definite.
\end{proof}

Q: Normal quadric
A: If $\mathbf{Z}\sim N_n(\mathbf{0}, \sigma^2 I)$ and $A$ is $n\times n$, symmetric, idempotent with rank $r$, then $\|A\mathbf{Z}\|^2=\mathbf{Z}^TA\mathbf{Z} \sim \sigma^2 \chi_r^2$.
\begin{proof}\leavevmode
Since $A$ is idempotent, the eigenvalues are either 0 or 1.
      Since $A$ is also symmetric, it is orthogonally diagonalisable as
      \[
        \Lambda = Q^TAQ = \diag(\underbrace{1, \cdots, 1}_r, \underbrace{0, \cdots, 0}_{n-r}).
      \]
Then $ Q^T\mathbf{Z}\sim N_n(\mathbf{0}, \sigma^2 I)$ and
      \begin{align*}
        \mathbf{Z}^TA\mathbf{Z} &=\mathbf{Z}^TQ\Lambda Q^T\mathbf{Z}\\
&= (Q^T\mathbf{Z})^T\Lambda Q^T\mathbf{Z}\\
&= \sum_{i=1}^r (Q^T\mathbf{Z})_i^2 \\
&\sim \sigma^2\chi_r^2.
      \end{align*}
\end{proof}

Q: Distribution of RSS
A: For NLMs
\[
\text{RSS}\sim \sigma^2 \chi_{n - p}^2, \quad \hat{\sigma}^2 \sim \frac{\sigma^2}{n}\chi_{n -p}^2.
\]
\begin{proof}
Since $(I_n - P)X = \mathbf{0}$:
\begin{align*}
\text{RSS} &= \|(I_n - P)\mathbf{Y} \|^2 \\
&= \|(I_n - P)\mathbf{\boldsymbol\varepsilon} \|^2  \\
&\sim \sigma^2 \chi^2_{n-p}.
\end{align*}
\end{proof}

Q: Prediction intervals
A: The expected response $\E[\mathbf{Y}\mid \mathbf{x}^*] = \mathbf{x}^{*T}\boldsymbol\beta$ is estimated by $\mathbf{x}^{*T}\hat{\boldsymbol\beta}$ where
\[
  \mathbf{x}^{*T}(\hat{\boldsymbol\beta} - \boldsymbol\beta)\sim  N(0, \sigma^2\mathbf{x}^{*T}(X^TX)^{-1}\mathbf{x}^*).
\]
Let $\tau^2 = \mathbf{x}^{*T} (X^TX)^{-1}\mathbf{x}^*$. Then
\[
  \frac{\mathbf{x}^{*T}(\hat{\boldsymbol\beta} - \boldsymbol\beta)}{\tilde{\sigma}\tau} = \frac{\mathbf{x}^{*T}(\hat{\boldsymbol\beta} - \boldsymbol\beta)/\sigma\tau}{\sqrt{\mathrm{RSS}/((n - p)\sigma^2)}}\sim t_{n - p} .
\]

Q: Independence of estimators in NLM
A: $\hat{\boldsymbol\beta}$ and $\hat{\sigma}^2$ are independent.
\begin{proof}\leavevmode
We can compute
\begin{align*}
\cov(\hat{\boldsymbol\beta}, \mathbf{R}) &= (X^TX)^{-1}X^T \var(Y) (I_n-P) \\
&= \sigma^2 (X^TX)^{-1}X^T (I_n-P) \\
&= 0.
\end{align*}
So $\hat{\boldsymbol\beta}$ and $\mathbf{R}$ are independent.
\end{proof}

Q: Inference for $\hat{\boldsymbol\beta}$
A: We know that $\hat{\boldsymbol\beta} \sim N_p(\boldsymbol\beta, \sigma^2(X^TX)^{-1})$. So
\[
  \hat{\beta}_j \sim N(\beta_j, \sigma^2(X^TX)^{-1}_{jj}).
\]
The \emph{standard error} of $\hat{\beta}_j$ is defined to be
\[
  \mathrm{SE}(\hat{\beta}_j) = \sqrt{\tilde{\sigma}^2 (X^TX)_{jj}^{-1}},
\]
where $\tilde{\sigma}^2 = \mathrm{RSS}/(n - p)$. Then
\[
  \frac{\hat{\beta}_j - \beta_j}{\mathrm{SE}(\hat{\beta}_j)} = \frac{\hat{\beta}_j - \beta_j}{\sqrt{\tilde{\sigma}^2(X^TX)^{-1}_{jj}}} = \frac{(\hat{\beta}_j - \beta_j)/\sqrt{\sigma^2(X^TX)_{jj}^{-1}}}{\sqrt{\mathrm{RSS}/((n - p)\sigma^2)}} \sim t_{n - p}.
\]

Q: Residual standard error
A: Since RSS $\sim \sigma^2\chi_{n-p}^2$, we know that $\E(\mathrm{RSS}) = \sigma^2(n - p)$.
Hence the RSE on $n - p$ degrees of freedom $\tilde{\sigma}^2 = \frac{\mathrm{RSS}}{n - p}$ is an unbiased estimator of $\sigma^2$.

Q: Confidence ellipsoids of NLM parameter
A: We know that
\[
\frac{\|X(\hat\beta - \beta)\|^2/p}{\tilde\sigma^2} \sim \frac{\sigma^2 \chi^2_p / p}{\sigma^2\chi_{n-p}^2 / (n-p)} \sim F_{p, n-p}.
\]
This gives the $1-\alpha$ confidence set
\[
\left\{ b\in\R^p: \frac{\|X(b-\hat\beta)\|^2/p}{\tilde\sigma^2} \leq F_{p, n-p}(\alpha) \right\}.
\]

Q: GLRT for NLMs
A: Let $\beta = (\beta_0, \beta_1)^T$ for $\beta_0\in\R^{p_0}$. The GLR is
\begin{align*}
  \Lambda_\mathbf{Y}(H_0, H_1) &= \frac{\left(\frac{1}{\sqrt{2\pi \hat{\sigma}^2}}\right)\exp\left(-\frac{1}{2\hat\sigma^2}(\|\mathbf{Y} - X\hat{\boldsymbol\beta}\|^2)\right)}{\left(\frac{1}{\sqrt{2\pi \hat\sigma_0^2}}\right)\exp\left(-\frac{1}{2\hat\sigma_0^2}(\|\mathbf{Y} - X\hat{\boldsymbol\beta}_0\|^2)\right)}\\
  &= \left(\frac{\hat\sigma_0^2}{\hat\sigma^2}\right)^{n/2} \\
  &= \left(\frac{\text{RSS}_0}{\mathrm{RSS}}\right)^{n/2} \\
  &= \left(1 + \frac{\mathrm{RSS}_0 - \mathrm{RSS}}{\mathrm{RSS}}\right)^{n/2}.
\end{align*}
We can apx $2\log \Lambda \sim \chi_{p - p_0}^2$. For an exact test, under $H_0$
\[
  \mathrm{RSS}_0 - \mathrm{RSS} = \boldsymbol\varepsilon^T(I_n - P_0)\boldsymbol\varepsilon - \boldsymbol\varepsilon^T(I_n - P)\boldsymbol\varepsilon = \boldsymbol\varepsilon^T(P - P_0)\boldsymbol\varepsilon
\]
where
\[
  (I_n - P)(P - P_0)  = 0.
\]
So we have the independent rvs
\begin{align*}
  \mathrm{RSS} = \boldsymbol\varepsilon^T(I_n - P)\boldsymbol\varepsilon &\sim \sigma^2\chi_{n - p}^2\\
  \mathrm{RSS}_0 - \mathrm{RSS} = \boldsymbol\varepsilon^T(P - P_0)\boldsymbol\varepsilon &\sim \sigma^2\chi^2_{p - p_0}.
\end{align*}
So under $H_0$,
\[
 \frac{(\mathrm{RSS}_0 - \mathrm{RSS})/(p - p_0)}{\mathrm{RSS}/(n - p)} \sim F_{p - p_0, n - p}.
\]

Q: One-way Analysis of Variance
A: Consider the model
\[
Y_{ij} = \mu + \alpha_j + \varepsilon_{ij}, \quad 1\leq i \leq n_j, \quad 1\leq j \leq J
\]
where each candidate is in exactly one of $J$ categories. Each \emph{categorical predictor} or \emph{factor} $j$ has $n_j$ candidates. If the $n_j$ are equal, we say the ANOVA is \emph{balanced}.
Write $Y = (Y_{11}, Y_{21}, \cdots Y_{n_1 1}, Y_{12} \cdots )^T$, $\beta = (\mu, \alpha_1, \cdots \alpha_n)$.
Then $r(X)=J<J+1$. However observe that we have the symmetry
\[
\mu\mapsto \mu + c, \quad \alpha_i \mapsto \alpha_i - c
\]
which does not affect $\hat Y$, $\hat\sigma^2$, MSPE or any testing procedures. Solutions:
\begin{enumerate}
\item set the \emph{corner point constraint} $\alpha_1=0$ (useful for when there is a natural control or baseline)
\item set $\mu = 0$,
\item \emph{sum-to-zero constraint} $\sum_{j=1}^J n_j\alpha_j=0$.
\end{enumerate}

Q: Two-way Analysis of Variance
A: Let each candidate be assigned one of $J$ factors and one of $K$ factors. The \emph{additive} model is
\[
Y_{ijk} = \mu + \alpha_j + \beta_k + \varepsilon_{ijk}, \quad 1\leq i \leq n_{jk}, 1 \leq j \leq J, 1 \leq k \leq K.
\]
Alternatively, we can add \emph{interaction} terms
\[
Y_{ijk} = \mu + \alpha_j + \beta_k + \gamma_{jk} + \varepsilon_{ijk}, \quad 1\leq i \leq n_{jk}, 1 \leq j \leq J, 1 \leq k \leq K.
\]
For identifiability of $(\mu, \alpha, \beta)$, constraints such as the corner point constraints $\alpha_1=\beta_1=\gamma_{1k}=\gamma_{j1}$ are necessary.

Q: Analysis of Covariance
A: A linear model that contains both categorical and numerical predictors.

Q: Coefficient of determination
A: Comparing to a model with just an intercept, we have
\[
R^2 = \frac{\|Y-\bar{Y}\mathbf{1}\|^2 - \|(I-P)Y\|^2}{\|Y-\bar{Y}\mathbf{1}\|^2}\in [0,1]
\]
with larger values indicating a better fit. This is the fraction of the sample variance explained by the predictors.

Q: Adjusted $R^2$
A: We can write
\[
R^2 = 1 - \frac{RSS}{RSS_0} = 1 - \frac{\hat\sigma^2}{\hat\sigma_0^2}.
\]
Using the unbiased versions of the variance estimators, we define
\begin{align*}
\tilde{R}^2 &=1 - \frac{RSS / (n-p)}{RSS_0 / (n-1)} \\
&= 1 - \frac{n-1}{n-p} (1-R^2) \\
&\in [-\infty, R^2].
\end{align*}
This is useful for comparing models with nested sets of predictors.

Q: AIC
A: Assume data $(Y_1, x_1^T), \cdots (Y_n, x_n^T)$ with $Y_i$ indep. Consider a model
\[
\mathcal{F}=\{(f_{x_i}(\ph, \theta))_{i=1}^n, \theta\in\Theta\subseteq\R^p\}
\]
with MLE $\hat{f}_{x_i}$ of true density $g_{x_i}$. Define the mean Kullback-Leibler divergence as
\[
\bar{K} = \frac{1}{n} \sum_{i=1}^n K(g_{x_i}, \hat{f}_{x_i}).
\]
We seek to minimise $\bar{K}\geq 0$, or equivalently minimise the variable part
\[
\tilde{K} = -\frac{1}{n} \sum_{i=1}^n \int_\R g_{x_i} \log \hat{f}_{x_i}.
\]
It can be shown that $\text{AIC}=-2\ell(\hat\theta)+2p$ satisfies
\[
\E[\text{AIC}] / n\approx 2 \E[\tilde{K}]
\]
for $n>>1$ provided $\mathcal{F}$ is well-specified.
In the case of a NLM
\[
\text{AIC} = n\left\{ 1 + \log(2\pi\hat\sigma^2)\right\} + 2(p+1) = n\log(\hat\sigma^2) + 2p + \text{const}.
\]

Q: Gibbs inequality
A: For continuous distributions $P,Q$, we define the \emph{Kullback-Leibler Divergence}
\[
D_{KL}(P||Q) = \int_\R p\log\frac{p}{q}.
\]
Then $D_{KL}\geq 0$ with equality iff $P=Q$ a.s.
\begin{proof}
Since $-\log$ is strictly convex, Jensen's inequality gives
\[
D_{KL}(P||Q) = \E_P\left[-\log\frac{q}{p}\right] \geq -\log\left(\E_P\left[\frac{q}{p}\right]\right) = - \log\left(\int_\R q\right)=0 .
\]
\end{proof}

Q: Orthogonal parameters
A: Let $X=(X_0, X_1), \beta = (\beta_0^T, \beta_1^T)^T$. We say that $\beta_0,\beta_1$ are \emph{orthogonal sets of parameters} if $X_0^TX_1=0$. Then
\begin{align*}
\hat\beta &= \left(\begin{pmatrix} X_0^T \\ X_1^T \end{pmatrix} (X_0 X_1)  \right)^{-1} \begin{pmatrix} X_0^T \\ X_1^T \end{pmatrix} Y \\
&= \begin{pmatrix} (X_0^TX_0)^{-1} & 0 \\
0 & (X_1^TX_1)^{-1}\end{pmatrix} \begin{pmatrix} X_0^T \\ X_1^T \end{pmatrix} Y \\
&= \begin{pmatrix} (X_0^TX_0)^{-1}X_0^TY \\
(X_1^TX_1)^{-1}X_1^TY\end{pmatrix} \begin{pmatrix} X_0^T \\ X_1^T \end{pmatrix} Y \\
&= \begin{pmatrix} \hat\beta_0 \\ \hat\beta_1 \end{pmatrix}.
\end{align*}
Suppose further the columns of $X$ are orthogonal. Write $P_S$ for the projection onto $\im (X_S)$ where $S\subseteq \{1, \cdots p\}$. By Pythagoras
\[
\| (I-P_S)Y\|^2 = \|Y\|^2 - \|P_{S}Y\|^2 =  \|Y\|^2 -\sum_{j\in S} \|\hat\beta_j  X_j\|^2.
\]
Hence to find the best model with $p_0\leq p$ variables, we list predictors in decreasing order of
\[
\|\hat\beta_j  X_j\|^2 = \frac{(X_j^TY)^2 \|X_j\|^2}{\|X_j\|^4} = \frac{(X_j^TY)^2}{\|X_j\|^2}.
\]

Q: Forward selection
A: \begin{enumerate}
\item Let $S_0$ be an intercept only model.
\item Add the predictor which reduces the RSS the most.
\item Terminate at large model size.
\item Pick one of the models $S_0\subseteq S_1\subseteq \cdots$ using AIC, $R^2$ os similar.
\end{enumerate}

Q: Backward selection
A: \begin{enumerate}
\item Let $S_0$ be the largest available model.
\item Remove the predictor which least increases the RSS.
\item Continue until a large number of predictors have been removed.
\item Pick a model from $S_0\supseteq S_1 \cdots$
\end{enumerate}

Q: Studentised residuals
A: The residuals $\hat\varepsilon=(I-P)Y\sim N(0, \sigma^2(I-P))$ can be standardised to
\[
\eta_i = \frac{\varepsilon_i}{\tilde{\sigma}\sqrt{1-p_i}}
\]
where $p_i=P_{ii}<1$ a.s. We hope these are apx iid $N(0,1)$. We can verify the small covariance by computing the average squared entry of $P$:
\[
\frac{1}{n^2} \sum_{i,j}P_{ij}^2 = \frac{1}{n^2} \tr(P^TP) = \frac{1}{n^2} \tr(P) = \frac{p}{n^2} \to 0
\]
as $n/p\to\infty$.
The name follows since
\[
 \frac{\varepsilon_i}{\tilde{\sigma}_{(-i)}\sqrt{1-p_i}} \sim t_{n-p-1}.
\]

Q: Assumptions of normal linear models
A: \begin{enumerate}
\item Centredness: $\E[\varepsilon]=0$.
\item Homoscedasticity: $\var(\epsilon_i)=\sigma^2$.
\item Uncorrelatedness: $\cov(\varepsilon_i, \varepsilon_j)=0$.
\item Normality: $\varepsilon \sim N(0, \sigma^2 I)$.
\end{enumerate}

Q: Approximating a Q-Q plot
A: If $U_1, \cdots U_n \sim U[0,1]$ are iid, then for any cts cdf $F$, $F^{-1}(U_i)\sim F$. So the $i$th order statistic can be simulated by $F^{-1}(U_{(i)})$.
We approximate
\[
\E[F^{-1}(U_{(i)})] \approx F^{-1} \E[U_{(i)}] = F^{-1}\left( \frac{i}{n+1}\right).
\]
Conclusion: plot $X_{(i)}$ against $F^{-1}\left( \frac{i}{n+1}\right)$. If $X_i\sim F$, then we expect a straight line through the origin with gradient 1.

Q: Order statistics of uniform
A: Let $U_1, \cdots U_n \sim U[0,1]$ be iid. Then $U_{(k)}\sim \text{Beta}(k, n+1-k)$.
\begin{proof}
$U_{(k)}$ has pdf
\[
f_{(k)}(x)\propto x^{k-1}(1-x)^{n-k}.
\]
\end{proof}

Q: Box-Cox transformations
A: If the dependent variables $y_i$ are not normally distributed, we can apply one of the transformations
\[
y \mapsto y^{(\lambda)} = \begin{cases}
\frac{y^\lambda - 1}{\lambda} \quad & \lambda\neq 0 \\
\log(y) \quad & \lambda = 0.
\end{cases}
\]
We pick $\lambda$ to maximise the log-likelihood.

Q: Leverage
A: The fitted values satisfy
\[
\hat{Y}_i = (PY)_i = P_{i1}Y_1 + \cdots + P_{in}Y_n.
\]
The value $p_i= P_{ii}$ is called the \emph{leverage} of the $i$th observation.
Since $\var(\hat\varepsilon_i) = \sigma^2(1-p_i)$, we see that $p_i\approx 1 \implies Y_i\approx \hat{Y}_i$.
Their sum is $\sum_{i=1}^n p_i = \tr(P) = p$. Rule of thumb: $p_i>3p/n$ is bad.

Q: Cook's distance: formula and algorithm
A: The Cooks distance $D_i$ of observation $(Y_i, x_i)$ is
\[
D_i = \frac{\tfrac{1}{p} \| X(\hat\beta_{(-i)} - \hat\beta) \|^2}{\tilde{\sigma}^2}.
\]
We do not need to fit $n+1$ models as we have from ES
\[
D_i = \frac{1}{p} \frac{p_i}{1-p_i} \hat\eta_i^2.
\]
Rule of thumb: $D_i>F_{p,n-p}(0.5)$.

Q: Model function
A: A probability density function or probability mass function

Q: Exponential tilting
A: Given non-degenerate model function $f_0$, cgf $K$, mgf $M$, \emph{natural parameter space} $\Theta=\{\theta: K(\theta)<\infty\}$, exponential tilting generates the \emph{natural exponential family}
\[
f(y, \theta) = \frac{e^{y\theta}f_0(y)}{\int e^{u\theta}f_0(u) \; \d u} = \frac{e^{y\theta}f_0(y)}{M(\theta)} = e^{\theta y - K(\theta)} f_0(y)
\]
for \emph{natural parameter} $\theta\in\Theta$.

Q: Cumulant generating function
A: The \emph{cumulant generating function} is logarithm of the MGF, i.e.
\[
K(t) = \log (M(t)).
\]
We write
\[
K(t) = \sum_{r=0}^\infty \kappa_r \frac{t^r}{r!}
\]
where $\kappa_r$ is the $r$th \emph{cumulant}.
In particular $\kappa_1=\E[Y]$, $\kappa_2=\var(Y)$.

Q: Cumulants for natural exponential family
A: If $Y\sim f(y,\theta)=e^{\theta y-K(\theta)} f_0(y)$ for $\theta\in\Theta$ natural, then
\[
\E_\theta[Y] = \frac{\d}{\d t} K(t,\theta)|_{t=0} = K'(\theta), \quad \var_\theta(Y) = \frac{\d^2}{\d t^2} K(t,\theta)|_{t=0} = K''(\theta).
\]
\begin{proof}
Since $\Theta$ is open, $\exists t$ such that $ \theta + t \in\Theta$ and thus
\begin{align*}
M(t,\theta) &= \int e^{ty} e^{\theta y-K(\theta)} f_0(y) \; \d y \\
&= e^{K(\theta + t) - K(\theta)}\int e^{(\theta + t) y-K(\theta + t)} f_0(y) \; \d y \\
&= e^{K(\theta + t) - K(\theta)} \\
&<\infty.
\end{align*}
Since $0\in\Theta$, the result follows.
\end{proof}

Q: Mean value parameterisation
A: If $Y\sim f(y,\theta)$, then we define the \emph{mean function} $\mu:\Theta\to \mathcal{M}$, where $\mathcal{M}=\mu(\Theta)$ is the \emph{mean space}, given by
\[
\mu(\theta):= \E_\theta[Y] = K'(\theta)
\]
satisfies $\mu'=K''=\var(Y)>0$. So $\mu$ is smooth and strictly increasing. Write $\theta(\mu)$ for the inverse.
Then we can reparameterise
\[
f(y, \mu) = e^{\theta(\mu)y - K(\theta(\mu))} f_0(y), \quad \mu\in\mathcal{M}.
\]
We also define the \emph{variance function} $V: \mathcal{M}\to (0,\infty) $ by
\[
V(\mu) = \var_{\theta(\mu)} Y = K''(\theta(\mu)).
\]

Q: Exponential dispersion family
A: A family of non-degenerate model functions
\[
f(y,\theta, \sigma^2) = a(\sigma^2, y) \exp\left[ \frac{\theta y - K(\theta)}{\sigma^2} \right]
\]
where
\begin{itemize}
\item $a(\ph, \ph)>0$ is known
\item the \emph{dispersion parameter} $\sigma^2\in\Phi\subseteq (0,\infty)$
\item the \emph{natural parameter} $\theta\in\Theta=(a,b)$.
\end{itemize}
Can show the CGF is
\[
K(t,\theta,\sigma^2) = \frac{K(\sigma^2 t + \theta) - K(\theta)}{\sigma^2}.
\]
Differentiating wrt $t$, we have the first two cumulants
\[
\mu(\theta)=K'(\theta), \qquad \sigma^2 V(\mu)= \sigma^2 K''(\theta(\mu)).
\]

Q: Generalised linear model
A: \emph{Model functions}. $Y_i$ independent model functions from same family
\[
f(y, \theta_i, \sigma_i^2) = a(\sigma_i^2, y) \exp\left[ \frac{1}{\sigma_i^2}\{ \theta_i y - K(\theta_i)\}\right] \quad \theta_i\in\Theta, \sigma_i^2\in\Phi\in(0,\infty).
\]
with $\sigma_i^2=a_i\sigma^2$ where $a_i>0$ are known.
\emph{Link function}. $\mu_i=\E[Y_i]$ and $\eta_i=x_i^T\beta$ are linked by
\[
\eta_i = g(\mu_i)
\]
where $g$ is twice differentiable strictly increasing. The \emph{canonical link function} is $g(\mu)=\theta(\mu)$, the natural parameter.

Q: Likelihood equations
A: The log-likelihood for a GLM is
\begin{align*}
\ell(\beta, \sigma^2) &= \sum_{i=1}^n \left[\frac{\theta(\mu_i) - K(\theta(\mu_i))}{\sigma^2a_i} + \log(a(\sigma_i^2,y_i))  \right]. \\
\intertext{Canonically $\theta(\mu_i)=x_i^T \beta$ so}
\ell(\beta, \sigma^2) &= \sum_{i=1}^n \left[\frac{y_ix_i^T\beta - K(x_i^T\beta)}{\sigma_i^2} + \log(a(\sigma_i^2,y_i)\right]. \\
\frac{\partial\ell(\beta,\sigma^2)}{\partial\beta} &= \sum_{i=1}^n \frac{x_i}{\sigma_i^2}\{y_i - \mu_i\} \\
\frac{\partial^2\ell(\beta,\sigma^2)}{\partial\beta\partial\beta^T} &= -\sum_{i=1}^n \frac{x_ix_i^T}{\sigma_i^2}V_i \prec 0.
\end{align*}

Q: Parametric inference in GLM
A: Let $i(\beta, \sigma^2)$ be the FIM of a GLM. Then $\exists i(\beta)\in\R^{p\times p}$ such that
\[
i(\beta, \sigma^2) = \begin{pmatrix}
i(\beta) & 0 \\
0 & i(\sigma^2)
\end{pmatrix},
\qquad i(\beta, \sigma^2)^{-1}=\begin{pmatrix}
i(\beta)^{-1} & 0 \\
0 & i(\sigma^2)^{-1}
\end{pmatrix}.
\]
From consistency of the MLE
\[
\hat\beta - \beta \sim AN(0, i^{-1}(\beta)).
\]
N.B. $i(\beta)$ actually still depends on $\sigma^2$, so may need to replace with $\tilde{\sigma}^2$.

Q: Estimating dispersion parameter
A: We can compute that
\[
\E[(Y_i-\mu_i)^2] = \sigma^2a_iV(\mu_i).
\]
so we can estimate with
\[
\tilde{\sigma}^2 = \frac{1}{n-p} \sum_{i=1}^n \frac{(y_i-\hat\mu_i)^2}{a_i V(\hat\mu_i)}.
\]

Q: GLRT for GLM
A: Let $\beta = (\beta_0^T, \beta_1^T)^T$, $H_0: \beta_1=0$, $\check\beta_0=\beta_\text{MLE}(H_0)$. Define
\[
\tilde\ell(\mu, \sigma^2) = \sum_{i=1}^n \left[ \frac{1}{a_i\sigma^2} \left\{(y_i\theta(\mu_i) - K(\theta(\mu_i))\right\} + \log(a(a_i\sigma^2, y_i))\right]
\]
and
\[
\hat{\mu}_i = g^{-1}(x_i^T\hat\beta), \qquad \check{\mu}_i = g^{-1}(x_i^T\check\beta).
\]
Taking maxima for e.g. the \emph{saturated model} and $H_1$:
\[
\tilde\ell(y, \sigma^2) = \max_{\mu\in\R^n} \tilde\ell(\mu, \sigma^2) , \qquad \ell(\hat\mu, \sigma^2) = \max_{\mu_i\in g^{-1}(x_i^T\R^p)} \tilde\ell(\mu, \sigma^2).
\]
Analogously to RSS, the \emph{deviances} of $H_0,H_1$ are
\begin{align*}
D(y, \check\mu) &= 2\sigma^2\left[ \tilde{\ell}(y,\sigma^2) - \ell(\check\mu, \sigma^2) \right] \\
D(y, \hat\mu) &= 2\sigma^2\left[ \tilde{\ell}(y,\sigma^2) - \ell(\hat\mu, \sigma^2) \right].
\end{align*}
By Wilk's Theorem
\[
w_{LR}(H_0) = \frac{D(y, \check\mu) - D(y, \hat\mu)}{\sigma^2} \sim \chi^2_{p-p_0}.
\]
When $\sigma^2$ is unknown
\[
 \frac{\tfrac{1}{p-p_0} (D(y, \check\mu) - D(y, \hat\mu))}{\tilde\sigma^2} \sim F_{p-p_0, n-p}.
\]

Q: Iterative weighted least squares
A: Applying Newton-Raphson to the score function $U=\frac{\partial\ell}{\partial\beta}$ gives
\[
\hat\beta^{(m)} = \hat\beta^{(m-1)} + j^{-1}(\hat\beta^{(m-1)})U(\hat\beta^{(m-1)}).
\]
To avoid singularity, Fisher replaces $j\mapsto i\succ 0$.
Define $W(\mu)=\diag(W_{ii})\in\R^{n\times n}, \tilde{Y}(\mu)\in\R^n$ where
\begin{align*}
 W_{ii}(\mu)&= \frac{1}{a_iV(\mu_i) \{g'(\mu_i)\}^2} \\
\tilde{Y}_i(\mu) &= g'(\mu_i)(y_i-\mu_i).
\end{align*}
Then
\begin{align*}
U(\beta) &= \sigma^{-2}X^TW(\mu)\tilde{Y}(\mu) \\
i(\beta)&= \sigma^{-2}X^TW(\mu)X.
\end{align*}
Further define $\hat\mu^{(m)}_i = g^{-1}(x_i^T\hat{\beta}^{(m)})$ and
\begin{align*}
W^{(m)} &= W(\hat\mu^{(m)}) \\
\tilde{Y}^{(m)} &= \tilde{Y}(\hat\mu^{(m)}).
\end{align*}
Then \emph{Fisher scoring} becomes
\[
\hat\beta^{(m)} = \hat\beta^{(m-1)} + ( X^T W^{(m-1)} X)^{-1} X^TW^{(m-1)}\tilde{Y}^{(m-1)}.
\]
Define the \emph{adjusted dependent variable}
\[
Z^{(m)}= \tilde{Y}^{(m)} +X\hat\beta^{(m)}.
\]
This gives the IWLS formula
\[
\hat\beta^{(m)} = ( X^T W^{(m-1)} X)^{-1}X^TW^{(m-1)} Z^{(m-1)} = \arg\min_{b\in\R^p} \left\{ \sum_{i=1}^n W_{ii}^{(m-1)}(Z_i^{(m-1)} - x_i^Tb)\right\}.
\]
Often use $\hat\mu^{(0)}=y$ or $\hat\mu^{(0)}=\max\{y, \varepsilon\}$.

Q: GLM residuals
A: The \emph{raw residuals} $Y_i-\hat\mu_i$ are not homoscedastic. The \emph{Pearson residuals} correct for this:
\[
e_i = \frac{y_i-\hat\mu_i}{\sqrt{a_iV(\hat\mu_i)}}.
\]
The \emph{deviance residuals} are
\[
d_i = \sgn (y_i-\hat\mu_i) \sqrt{D(y_i,\hat\mu_i)}.
\]

Q: Binomial regression is a GLM
A: Let $Y_i\sim\frac{1}{n_i}\text{Bin}(n_i,\mu_i)$. Then
\begin{align*}
f(y_i, \mu_i) &= {n_i \choose n_i y_i} \mu_i^{n_iy_i} (1-\mu_i)^{n_i(1-y_i)} \\
&= \underbrace{ {n_i \choose n_iy_i} }_{a(a_i,y_i)} \exp\bigg[ \frac{1}{n_i^{-1}} \bigg\{ y_i \underbrace{\log\left( \frac{\mu_i}{1-\mu_i} \right)}_{\theta_i=\theta(\mu_i)} + \underbrace{\log(1-\mu_i)}_{K(\theta_i)} \bigg\} \bigg].
\end{align*}
We take dispersion parameter 1 and $a_i=n_i^{-1}$.

Q: Binomial link functions
A: \begin{enumerate}
\item \emph{Complementary log-log link} $\log(-\log(1-\mu))$
\item \emph{Probit link} $\Phi^{-1}(\mu)$
\item \emph{Logit link} (canonical) $\log(\frac{\mu}{1-\mu})$
\end{enumerate}

Q: Logistic regression coefficients
A: The value $e^{\beta_j}$ from a GLM with logit link gives the multiplicative change in the odds $\frac{\mu_i}{1-\mu_i}$ for a unit increase in the $j$th predictor.
\begin{proof}
\[
\frac{\mu_i}{1-\mu_i} = \exp\left(\sum_{j=1}^p X_{ij}\beta_j\right) = \prod_{j=1}^p(e^{\beta_j})^{X_{ij}}.
\]
\end{proof}

Q: Poisson regression link
A: Log link $\log(\mu)$, aka \emph{log-linear model}.
Then $\log(0,\infty)=\R$ and
\[
\mu_i = \exp\left( \sum_{j=1}^pX_{ij}\beta_j \right) = \prod_{j=1}^p (e^{\beta_j})^{X_{ij}}
\]
so $e^{\beta_j}$ is the multiplicative change in the expected response for a unit increase in the $j$th predictor.

Q: Poisson regression likelihood equations
A: With $\log(\mu_i)=x_i^T\beta$
\begin{align*}
L(\beta) &= \prod_{i=1}^n e^{-\mu_i(\beta)} \frac{(\mu_i\beta)^y_i}{y_i!} \\
\ell(\beta) &= \sum_{i=1}^n \left[-\exp(x_i^T\beta) + y_ix_i^T\beta\right] + \text{const}. \\
\intertext{Letting $x_i^T\beta\mapsto \alpha + x_i^T\beta$ and differentiating}
\ell(\alpha, \beta) &= \sum_{i=1}^n \left[-\exp(\alpha+x_i^T\beta) + y_i(\alpha+x_i^T\beta)\right] + \text{const}. \\
\frac{\partial\ell}{\partial\alpha} &=\sum_{i=1}^n \left[-\exp(\alpha+x_i^T\beta) +  y_i \right]. \\
\intertext{Writing $\hat\mu_i = \exp(\hat\alpha+x_i^T\hat\beta)$:}
\sum_{i=1}^n\hat\mu_i &= \sum_{i=1}^ny_i.
\end{align*}

Q: Poisson regression deviance
A: Assuming an intercept is included
\begin{align*}
\tilde\ell(\mu, \sigma^2) &= -\sum_{i=1}^n \mu_i +\sum_{i=1}^n y_i \log (\mu_i) \\
D(y, \hat\mu) &= 2\sum_{i=1}^n y_i\log(\frac{y_i}{\hat\mu_i}) - 2\sum_{i=1}^n (y_i-\hat\mu_i) \\
&= 2\sum_{i=1}^n y_i\log(\frac{y_i}{\hat\mu_i}) \\
&=: 2\sum_{i=1}^n (\hat\mu_i+\delta_i)\log(1+\frac{\delta_i}{\hat\mu_i}) \\
&= 2\sum_{i=1}^n \left[ \delta_i+\frac{\delta_i^2}{\hat\mu_i} - \frac{\delta_i^2}{2\mu_i}\right] + O(\delta^3)  \\
&= \sum_{i=1}^n \frac{(y_i-\hat\mu_i)^2}{\hat\mu_i} + O(\delta^3) \\
&\sim \chi^2_{n-p}.
\end{align*}

Q: Poisson regression with an offset
A: If $Y_i\propto t_i$ and $Y_i\sim \Po(\mu_i)$, then
\[
\log(\frac{\mu_i}{t_i}) = \log(\mu_i) - \log(t_i) = x_i^T\beta
\]
so
\[
\mu_i = t_i\exp(x_i^T\beta).
\]
This is the usual Poisson regression with an \emph{offset} of $log(t_i)$.

Q: Poisson with fixed total
A: Let $(Z_i)_1^m \sim \Po(\mu_i)$. Conditional on $\sum Z_i=n$, we have that
\[
Z\sim \mathrm{Multinomial}(n;p_1, \cdots p_m)
\]
where $p_i = \mu_i / \sum_j \mu_j$.
\begin{proof}
If $\sum z_j=n$, since $\sum Z_j \sim \Po(\sum \mu_j)$, using Bayes
\begin{align*}
\P(\forall i: Z_i=z_i | \sum Z_j = n) &= \frac{\exp(-\sum \mu_j) \prod (\mu^{z_i} / z_i!)}{\exp(-\sum \mu_j)(\sum \mu_j)^n / n!} \\
&= \frac{n!}{z_1! \cdots z_m!} p_1^{z_1} \cdots p_m^{z_m}.
\end{align*}
\end{proof}

Q: Poisson multinomial likelihood
A: Let $(Y_{ij})_{i\leq I, j\leq J}\sim \mathrm{Multinomial} (n; (p_{ij}))$ where $p_{ij}=\mu_{ij}/\mu_{++}$ and $\log(\mu_{ij})=\alpha + x_{ij}^T\beta$. Then
\[
p_{ij} = \frac{\exp(x_{ij}^T\beta)}{\sum \exp(x_{k\ell}^T\beta)}.
\]
Parameterise as $x_{ij} = \mathbbm{1}_{\{i\} \cup \{I+j\}}$. Then the log-likelihood is
\begin{align*}
\ell_m(\beta | n) &= \sum_{i,j} y_{ij} \log p_{ij}(\beta) + \text{constant} \\
&= \sum_{i,j} y_{ij}x_{ij}^T\beta -n\log\left( \sum_{k,\ell} \exp(x_{k\ell}^T\beta) \right) + \text{constant}.
\end{align*}

Q: Poisson surrogate likelihood
A: Let $(Y_{ij})\sim\Po(\mu_{ij})$ where $\log(\mu_{ij})=\alpha + x_{ij}^T\beta$. If $\sum_{ij}y_{ij}=n$, then the log-likelihood is
\begin{align*}
\ell_P(\alpha, \beta) &= - \sum_{i,j} \mu_{ij}(\alpha,\beta)+ \sum_{i,j} y_{ij} \log \mu_{ij}(\alpha, \beta) + \text{constant} \\
&= - \sum_{i,j}\exp(\alpha + x_{ij}^T\beta)+ \sum_{i,j} y_{ij} (\alpha + x_{ij}^T\beta) + \text{constant} \\
&= - \exp(\alpha) \sum_{i,j}\exp(x_{ij}^T\beta)+ n\alpha + \sum_{i,j} y_{ij}  x_{ij}^T\beta + \text{constant} \\
\ell(\tau, \beta) &= \sum_{i,j} y_{ij}  x_{ij}^T\beta - n\log\left( \sum_{ij} \exp(x_{ij}^T\beta)\right) + n\log\tau - \tau \\
&= \ell_m(\beta | n) + \ell_P(\tau)
\end{align*}
where we reparemeterise using
\[
\tau = \sum_{ij}\mu_{ij} = \exp(\alpha) \sum_{i,j}\exp(x_{ij}^T\beta).
\]
In particular $\hat\beta_m=\hat\beta_P$.

Q: Surrogate Poisson models
A: The Poisson and Multinomial models have the same deviances, FIM and fitted values.
Hence we can fit Multinomial models using Poisson log-linear models with an intercept.

Q: Three-way contigency table
A: Eight hypotheses for independence:
\begin{enumerate}
\item $p_{ijk}=\alpha_i\beta_j\gamma_k$, i.e. $(A,B,C)$ indep [1]
\item $p_{ijk}=\alpha_i\beta_{jk}$, i.e. $A$ indep of $(B,C)$ [3]
\item $p_{ijk}=\beta_{ij}\gamma_{ik}$, i.e. given $A$; $B,C$ are indep [3]
\item $p_{ijk}=\alpha_{jk}\beta_{ik}\gamma_{ij}$: no three-way interactions. [1]
\end{enumerate}

Q: Observed information matrix
A: Define
\[
j(\theta) = -\frac{\partial^2}{\partial\theta\partial\theta^T} \ell(\theta, Y).
\]
such that $i(\theta)=\E_\theta [ j(\theta)]$.