Linear models

• We’ve seen several linear models in the course, most likelihood based…

• Generic setup: independent \(Y_i\)’s parameterized by \(\eta_i\).

• Linear model imposes \(\eta_i=X_i^T\beta\).

• Deviance (-2 log likelihood) provides objective function to estimate \(\beta\).

• These models have nice structure when \(Y_i\)’s are from some exponential family: often convex (negative log-likelihoods) which make it easy to regularize, go to high dimensions, etc.

Likelihoods we’ve looked at

1. Binomial

2. Poisson

3. Multinomial

4. Cox partial likelihood

Other models: multivariate Gaussian

• Setup: \(Y_{n \times q}\), \(X_{n \times p}\)

• Model: \(Y|X \sim N(XB, I_{n \times n} \otimes \Sigma_{q \times q})\) for \(B_{p \times q}\) (covariance is multivariate Kronecker notation for IID draws from \(N(0,\Sigma)\).

• A likelihood based model: MLE (exactly analogous to univariate reponse):

\[ \begin{aligned} \widehat{B} &= (X^TX)^{-1}X^TY \\ \widehat{\Sigma} &= \frac{1}{n}Y^T(I - X(X^TX)^{-1}X^T)Y \end{aligned} \]

• Can regularize in groups for group sparsity

\[ {\cal P}(B) = \sum_j \|B[j,]\|_2 \]

Other models: robust regression rlm

• Setup: \(Y_{n \times 1}\), \(X_{n \times p}\)

• Model: typically not thought of as a likelihood…

• Estimator:

\[ (\widehat{\beta}, \widehat{\sigma}) = \text{argmin}_{\beta} \sum_{i=1}^n \rho((Y_i-X_i^T\beta)/\sigma) + n \log(\sigma) \]

Other models: robust regression rlm

• Examples of \(\rho\): \(\ell_1\), Huber, Tukey biweight or other redescending estimator.

• Regularization? Sure… (Huber and \(\ell_1\) are at least convex…)

Other models: quantile regression

• Consider \(\ell_1\):

\[ \text{minimize}_{\beta} \sum_{i=1}^n |Y_i-X_i'\beta| \]

• Instead of modeling \(E[Y|X]\) this models \(\text{median}(Y|X)\).

Pinball loss: \(\tau \in (0,1)\)

\[ \ell_{\tau}(Y,\eta) = \begin{cases} \tau (Y-\eta) & Y > \eta \\ -(1-\tau)(Y-\eta) & Y\leq \eta \end{cases} \]

• Models the \(\tau\)-quantile of \(Y|X\).

Other models: robust regression rlm

• Score equations

\[ \begin{aligned} X^T\psi((Y - X\widehat{\beta})/\widehat{\sigma}) &= 0 \\ \sum_{i=1}^n\psi((Y_i - X_i^T\widehat{\beta})/\widehat{\sigma}) (Y_i - X_i^T\widehat{\beta})/\widehat{\sigma} &= n \end{aligned} \]

• Inference? Sandwich or bootstrap.

Other models: support vector machine (SVM)

• Not always phrased this way, but equivalent to using hinge loss

• For binary outcome \(Y\) coded as \(\pm 1\) and linear predictor \(\eta\)

\[ \ell_{\text{Hinge}}(Y,\eta) = \max(0, 1 - Y \cdot \eta) \]

• Huberized SVM smooths out this loss as Huber does the \(\ell_1\) (Moreau smoothing)

Other models: support vector machine (SVM)

• Regularization: convex in \(\eta\) so, sure…

• Inference? Bootstrap?

Mixed effects models

• Setup: \(Y_{n \times 1}\), \(X_{n \times p}\)

• Model: \(Y|X \overset{D}{=} X\beta + Z\gamma + \epsilon\)

• Measurement error \(\epsilon \sim N(0, \sigma^2 I)\)

• Random effect (usually independent of \(\epsilon\)) \(\gamma \sim N(0, \Sigma(\theta))\).

Mixed effects models

• Log-likelihood

\[ \begin{aligned} \ell(\beta,\sigma,\theta | Y,X) &= -\frac{1}{2} (Y-X\beta)^T(\sigma^2 I + Z\Sigma(\theta)Z^T)^{-1}(Y-X\beta) \\ & \qquad - \frac{1}{2} \log\left(\det(\sigma^2 I + Z\Sigma(\theta)Z^T)\right) \end{aligned} \]

• Estimation: not convex but smooth…

• Inference: likelihood based (if you believe the model)

Mixed effects models: estimation

• EM algorithm (if we observed \(\gamma\)) likelihood would be nice and simple.

• Gradient descent…

• Messy, but feasible to run Newton-Raphson…

• Since we have a likelihood, Bayesian methods are applicable here.

Mixed effects models: REML

• REML: notes that law of residuals does not depend on \(\beta\)

\[ R = (I - X(X^TX)^{-1}X^T)Y \]

• Decomposes problem into noise parameter estimation and fixed effects estimation.

• This likelihood can be used to estimate \((\theta,\sigma^2)\).

• Given this estimate \((\widehat{\sigma}^2, \widehat{\theta})\) the optimal \(\beta\) is

\[ (X^T(\widehat{\sigma}^2I + \Sigma(\widehat{\theta}))^{-1}X)^{-1} X^T(\widehat{\sigma}^2I + \Sigma(\widehat{\theta}))^{-1}Y \]

• REML is often a better estimate than MLE for noise: in usual Gaussian model \(\widehat{\sigma}^2_{REML} = \|Y-X\widehat{\beta}\|^2_2/(n-p)\)…

Generalized linear mixed models

• Gaussian mixed effects model essentially models linear predictor

\[ \begin{aligned} \eta &= X\beta + Z\gamma \\ &= \text{fixed effect} + \text{noise} \end{aligned} \]

• Lends itself naturally to glm likelihoods as well.

• E.g. for Binomial data with logistic link

\[ Y|X,Z,\gamma \sim \text{Binomial}\left(\frac{e^{\eta}}{1+e^{\eta}}\right) \]

Generalized linear mixed models

• Likelihood marginalizes over \(\gamma\) – no closed form (unlike Gaussian model).

• Many computational issues to deal with here.

• Well suited to Bayesian modeling here (though if optimization is difficult, likely that posteriors are quite multi-modal for non-informative priors…)

Beyond linear models: basis expansions

• We’ve always used linear predictor

\[ \eta = X\beta = \sum_{j=1}^p \beta_j X_j \]

• Could use a basis expansion with \(N_j\) functions for feature \(j\)

\[ \eta = \sum_{j=1}^p \sum_{l=1}^{N_j} \beta_{jl} h_{jl}(X_j) \]

• Any linear method that takes linear predictor can profit from basis expansion.

• Generally, these form additive models

\[ \eta = \sum_{j=1}^p f_j(X_j) \]

Implicit basis expansions

• As basis expansion grows in complexity we should regularize in some fashion…

• For explicit basis, we could choose sparsity or some ridge

\[ {\cal P}_j\left(\sum_{l=1}^{N_j} \beta_{jl} h_{jl}\right) = \sum_{l=1}^{N_j} \beta_{jl}^2 \]

Kernel trick

• A natural class of penalties

• A covariance function \(K:\mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}\) is a function that satisfies:

  1. \(K(s,t)=K(t,s)\) (symmetry)

  2. For any \(\{t_1, \dots, t_k\} \subset\mathbb{R}\) and coefficients \(\{a_1, \dots, a_k\}\) we have \(\sum_{i,j=1}^k a_i a_j K(t_i, t_j)>0\). (non-negative definite)

• Example: radial basis function \(K(t,s) = e^{-(t-s)^2/(2\gamma)}\).

• Example: Brownian motion \(K(t,s) = \min(t,s)\).

Reproducing kernel penalty

• Consider functions of the form:

\[ h(s) = \sum_{j=1}^k a_j K(t_j,s) \]

• The penalty is

\[ \|h\|^2_{K} = \sum_{l,j=1}^k a_i a_j K(t_i,t_j) \]

Penalized regression

\[ \text{minimize}_{f_1, \dots, f_p} -\log L(f_1,\dots,f_p|X,Y) + \sum_{j=1}^p \lambda_j \|f_j\|^2_{K_j}. \]

• Kernels can take multivariate arguments \(K:\mathbb{R}^k \times \mathbb{R}^k \rightarrow \mathbb{R}\).

• Greatly expands flexibility of regression modelling…. (at cost of having more choices to make!)

Combining ideas

• Kernel trick tells us how to make lots of cool feature mappings…

• Can be plugged into other regression problems, e.g. “non-parametric” quantile regression

\[ \text{minimize}_f \sum_{i=1}^n \ell_{\tau}(Y_i,f(X_i)) + \lambda \|f\|^2_{K}. \]

• “Non-parametric” survival analysis

\[ \text{minimize}_f \sum_{i:\delta_i=1} \left[ \log \left(\sum_{j:O_j \geq O_i} \exp(f(X_j)) \right) - f(X_i)\right] + \lambda \|f\|^2_K. \]

Review

Key topics from each chapter

Not exhaustive

Chapter 1

• Key distributions for categorical data: Poisson, Binomial, Multinomial.

• Examples of exponential families.

• Inference for sampling from such families.

Chapter 2

• Contingency tables.

• Sampling models: independent rows / columns; Poisson; case-control.

• Odds ratios.

• Conditional independence in three-way tables.

Chapter 3

• Inference for two-way tables.

• Tests of independence (Pearson / Fisher)

• Residuals: Pearson and deviance.

• Ordinal variables.

• Bayesian models

Chapter 4

• Generalized Linear Models

• Basic definitions and examples: binary, counts

• Likelihood / Score / Hessian (Information)

• Deviance: likelihood ratio tests

• Residuals

• Algorithm: Fisher scoring / Newton-Raphson

Chapter 5

• Odds ratios

• Retrospective studies

• Inference

• A little bit of inference under misspecification: e.g. bootstrap

• Multiple logistic regression

Chapter 6

• Model selection and diagnostics (deviance / Pearson residuals) for logistic regression

• Analogs of \(R^2\)

• Mantel-Haenszel tests for stratified 2x2 tables

• Power in logistic regression

Bayesian GLMs (Touched on in Chapter 7 of Agresti)

• Brief description of Metropolis-Hastings / Gibbs sampling

• Simple implementations in STAN

Chapter 8

• Multinomial regression

• Multicategory logit

• Cumulative logit

• Latent variable representations (data augmentation / Bayesian)

• Ordinal responses

• Conditional independence in 3-way tables

Chapter 9

• Loglinear models

• Homogeneous association (interactions up to 2-way only)

• Graphical models

• Inference (deviance, Wald)

• Connection to logistic regression

• Lindsey’s method: density modeling as GLM

Chapter 10

• Model selection for loglinear models

• Deviance tests

• Residuals for diagnostics

Chapter 11

• Matched pairs and square contingency tables

• McNemar’s test

• Subject specific models / conditonal inference

• Symmetry and quasi-symmetry (Bradley-Terry)

• Quasi-independence (agreement between raters)

Regularization

• For GLMs we started off talking about LASSO: \[ {\cal P}(\beta) = \lambda \|\beta\|_1 = \lambda \sum_{j=1}^p |\beta_j|. \]

• KKT conditions: solving LASSO yields pair \((\hat{\beta}, \widehat{u}=-\nabla L(\hat{\beta}))\)

Regularization

LASSO picture

Regularization

Solving LASSO: coordinate descent

\[ \hat{\beta}^{(t+1)}_j = \text{argmin}_{\beta_j \in \mathbb{R}} \left[\beta_j \mapsto \Lambda (X_{-j}\hat{\beta}^{(t)}_{-j} + X_j \beta_j) - \beta_j \cdot (X_j^TY) + p_j(\beta_j) \right] \]

• Sort of analogous to Gibbs sampling but instead of “drawing from dbn” we find its mode.

Regularization

Solving LASSO: proximal gradient

\[ \hat{\beta}^{(t+1)} = S_{\lambda / L}\left( \hat{\beta}^{(t)} - \frac{1}{L} \nabla \ell(\hat{\beta}^{(t)})\right) \]

• \(S_{\lambda}\) is soft-threshold operator (proximal operator of \(\ell_1\) penalty)

• Updates all coordinates at once

Survival analysis

• Lifetime data, often right-censored.

• Observations are \((O_i, \delta_i)\) with \(O_i=\min(T_i, C_i)\) and \[ \delta_i = \begin{cases} 1 & T_i \leq C_i \\ 0 & \text{otherwise.} \end{cases} \]

• Goal: understand distribution of \(T|X\) based on response \((O,C)\).

Survival analysis

Hazard

• Basic object modeled in survival analysis

• Presuming density of \(T\) \[ 1 - F(t) = S(t) = \exp(-H(t)), \qquad H(t) = \int_{[0,t]} h(s) \; ds \]

• Estimates w/o features: Nelson-Aalen for hazard, Kaplan-Meier for survival

• Log-rank tests to compare populations

Survival analysis

Accelerated failure time model

\[ S_{\beta}(t;x) = S_0\left(\frac{t}{\theta(x^T\beta)}\right). \]

• Common choice of link \(\theta(u)=e^{-u}\).

Survival analysis

Proportional hazards model

• Cox model (fit with partial likelihood) \[ H_{\beta}(t;x) = e^{x^T\beta} H_0(t) \]

• Expressible in terms of linear predictor \(\eta_i=X_i^T\beta\).

• Allows linear modeling / regularization / etc.

• Various residuals for assessing goodness of fits.